##' @importFrom ggplot2 fortify
 ##' @method fortify treedata
 ##' @export

 fortify.treedata <- function(model, data,
                              layout        = "rectangular",
                              yscale        = "none",
                              ladderize     = TRUE,
                              right         = FALSE,
                              branch.length = "branch.length",
                              mrsd          = NULL,
                              as.Date       = FALSE, ...) {

     model <- set_branch_length(model, branch.length)

     fortify.phylo(model, data,
                   layout        = layout,
                   yscale        = yscale,
                   ladderize     = ladderize,
                   right         = right,
                   branch.length = branch.length,
                   mrsd          = mrsd,
                   as.Date       = as.Date, ...)
 }
 
 ##' @importFrom ape ladderize
 ##' @method fortify phylo
 ##' @export
 fortify.phylo <- function(model, data,
                           layout        = "rectangular",
                           ladderize     = TRUE,
                           right         = FALSE,
                           branch.length = "branch.length",
                           mrsd          = NULL,
                           as.Date       = FALSE,
                           yscale        = "none",
                           ...) {
 
     x <- as.phylo(model) ## reorder.phylo(get.tree(model), "postorder")

     if (ladderize == TRUE) {
         x <- ladderize(x, right=right)
     }

 
     if (! is.null(x$edge.length)) {
         if (anyNA(x$edge.length)) {
             warning("'edge.length' contains NA values...\n## setting 'edge.length' to NULL automatically when plotting the tree...")
             x$edge.length <- NULL
         }
     }
 

     if (is.null(x$edge.length) || branch.length == "none") {
         xpos <- getXcoord_no_length(x)
     } else {
         xpos <- getXcoord(x)
     }

     ypos <- getYcoord(x)
     N <- Nnode(x, internal.only=FALSE)

     xypos <- data_frame(node=1:N, x=xpos, y=ypos)

     df <- as_data_frame(model)

     res <- full_join(df, xypos, by = "node")

 
     ## add branch mid position
     res <- calculate_branch_mid(res)
 

     if (!is.null(mrsd)) {
         res <- scaleX_by_time_from_mrsd(res, mrsd, as.Date)
     }
 

     if (layout == "slanted") {
         res <- add_angle_slanted(res)
     } else {
         ## angle for all layout, if 'rectangular', user use coord_polar, can still use angle
         res <- calculate_angle(res)
     }

     scaleY(as.phylo(model), res, yscale, layout, ...)

 }
 

 ##' @method as_data_frame treedata
 ##' @importFrom tibble as_data_frame

 ##' @export

 ##' @importFrom treeio Nnode
 ##' @importFrom treeio Ntip
 as_data_frame.treedata <- function(x, ...) {
     res <- as_data_frame(x@phylo)

     tree_anno <- as_data_frame(get_tree_data(x))

     if (nrow(tree_anno) > 0) {

         by <- "node"
         tree_anno$node <- as.integer(tree_anno$node)
         if ("parent" %in% colnames(tree_anno)) {
             by <- c(by, "parent")
             tree_anno$parent <- as.integer(tree_anno$parent)
         }
 
         res <- full_join(res, tree_anno, by=by)

     }
     return(res)
 }
 

 ##' @method as_data_frame phylo
 ##' @export

 ##' @importFrom tibble data_frame
 ##' @importFrom dplyr full_join

 as_data_frame.phylo <- function(x, ...) {

     phylo <- x
     ntip <- Ntip(phylo)
     N <- Nnode(phylo, internal.only=FALSE)
 
     tip.label <- phylo[["tip.label"]]

     res <- as_data_frame(phylo[["edge"]])

     colnames(res) <- c("parent", "node")
     if (!is.null(phylo$edge.length))
         res$branch.length <- phylo$edge.length
 
     label <- rep(NA, N)
     label[1:ntip] <- tip.label
     if ( !is.null(phylo$node.label) ) {
         label[(ntip+1):N] <- phylo$node.label
     }
     isTip <- rep(FALSE, N)
     isTip[1:ntip] <- TRUE

 
     label.df <- data_frame(node=1:N, label=label, isTip = isTip)
     res <- full_join(res, label.df, by='node')
 
     idx <- is.na(res$parent)
     res$parent[idx] <- res$node[idx]

     res <- res[order(res$node),]
     aa <- names(attributes(phylo))
     group <- aa[ ! aa %in% c("names", "class", "order", "reroot", "node_map")]
     if (length(group) > 0) {
         for (group_ in group) {
             ## groupOTU & groupClade
             group_info <- attr(phylo, group_)
             if (length(group_info) == nrow(res)) {
                 res[[group_]] <- group_info
             }
         }
     }

     return(res)
 }
 
 get_tree_data <- function(tree_object) {

     tree_anno <- tree_object@data
     extraInfo <- tree_object@extraInfo

     if (nrow(tree_anno) == 0) {
         return(extraInfo)

     }

     if (nrow(extraInfo) == 0) {
         return(tree_anno)
     }
 
     full_join(tree_anno, extraInfo, by = "node")

 }
 
 

 ##' convert tip or node label(s) to internal node number
 ##'
 ##'
 ##' @title nodeid
 ##' @param x tree object or graphic object return by ggtree
 ##' @param label tip or node label(s)
 ##' @return internal node number
 ##' @export
 ##' @author Guangchuang Yu
 nodeid <- function(x, label) {
     if (is(x, "gg"))
         return(nodeid.gg(x, label))
 
     nodeid.tree(x, label)
 }
 
 nodeid.tree <- function(tree, label) {
     tr <- get.tree(tree)
     lab <- c(tr$tip.label, tr$node.label)
     match(label, lab)
 }
 
 nodeid.gg <- function(p, label) {
     p$data$node[match(label, p$data$label)]
 }
 
 
 reroot_node_mapping <- function(tree, tree2) {
     root <- getRoot(tree)
 
     node_map <- data.frame(from=1:getNodeNum(tree), to=NA, visited=FALSE)
     node_map[1:Ntip(tree), 2] <- match(tree$tip.label, tree2$tip.label)
     node_map[1:Ntip(tree), 3] <- TRUE
 
     node_map[root, 2] <- root
     node_map[root, 3] <- TRUE
 
     node <- rev(tree$edge[,2])
     for (k in node) {
         ip <- getParent(tree, k)
         if (node_map[ip, "visited"])
             next
 
         cc <- getChild(tree, ip)
         node2 <- node_map[cc,2]
         if (anyNA(node2)) {
             node <- c(node, k)
             next
         }
 
         to <- unique(sapply(node2, getParent, tr=tree2))
         to <- to[! to %in% node_map[,2]]
         node_map[ip, 2] <- to
         node_map[ip, 3] <- TRUE
     }
     node_map <- node_map[, -3]
     return(node_map)
 }
 
 
 
 ##' @importFrom ape reorder.phylo

 layout.unrooted <- function(tree, branch.length="branch.length", layout.method="equal_angle", ...) {

     df <- switch(layout.method,
                  equal_angle = layoutEqualAngle(tree, branch.length),
                  daylight = layoutDaylight(tree, branch.length)
                  )

     return(df)
 }

 ##' 'Equal-angle layout algorithm for unrooted trees'
 ##'

 ##' @references

 ##' "Inferring Phylogenies" by Joseph Felsenstein.

 ##'

 ##' @title layoutEqualAngle
 ##' @param tree phylo object
 ##' @param branch.length set to 'none' for edge length of 1. Otherwise the phylogenetic tree edge length is used.
 ##' @return tree as data.frame with equal angle layout.
 layoutEqualAngle <- function(tree, branch.length ){

     root <- getRoot(tree)
     ## Convert Phylo tree to data.frame.
     df <- as.data.frame.phylo_(tree)
 
     ## NOTE: Angles (start, end, angle) are in half-rotation units (radians/pi or degrees/180)
 
     ## create and assign NA to the following fields.
     df$x <- NA
     df$y <- NA
     df$start <- NA # Start angle of segment of subtree.
     df$end   <- NA # End angle of segment of subtree

     df$angle <- NA # Orthogonal angle to beta for tip labels.

     ## Initialize root node position and angles.
     df[root, "x"] <- 0
     df[root, "y"] <- 0
     df[root, "start"] <- 0 # 0-degrees
     df[root, "end"]   <- 2 # 360-degrees
     df[root, "angle"] <- 0 # Angle label.
 
     N <- getNodeNum(tree)
 
     ## Get number of tips for each node in tree.
     nb.sp <- sapply(1:N, function(i) length(get.offspring.tip(tree, i)))
     ## Get list of node id's.
     nodes <- getNodes_by_postorder(tree)
 
     for(curNode in nodes) {
         ## Get number of tips for current node.
         curNtip <- nb.sp[curNode]
         ## Get array of child node indexes of current node.
         children <- getChild(tree, curNode)
 
         ## Get "start" and "end" angles of a segment for current node in the data.frame.
         start <- df[curNode, "start"]
         end <- df[curNode, "end"]
 
         if (length(children) == 0) {
             ## is a tip
             next
         }

         for (i in seq_along(children)) {
             child <- children[i]
             ## Get the number of tips for child node.
             ntip.child <- nb.sp[child]
 
             ## Calculated in half radians.
             ## alpha: angle of segment for i-th child with ntips_ij tips.
             ## alpha = (left_angle - right_angle) * (ntips_ij)/(ntips_current)
             alpha <- (end - start) * ntip.child / curNtip
             ## beta = angle of line from parent node to i-th child.
             beta <- start + alpha / 2
 
             if (branch.length == "none") {
                 length.child <- 1
             } else {
                 length.child <- df[child, "length"]
             }
 
             ## update geometry of data.frame.
             ## Calculate (x,y) position of the i-th child node from current node.
             df[child, "x"] <- df[curNode, "x"] + cospi(beta) * length.child
             df[child, "y"] <- df[curNode, "y"] + sinpi(beta) * length.child

             ## Calculate orthogonal angle to beta for tip label.

             df[child, "angle"] <- -90 - 180 * beta * sign(beta - 1)
             ## Update the start and end angles of the childs segment.
             df[child, "start"] <- start
             df[child, "end"] <- start + alpha
             start <- start + alpha
         }

 
     }

     return(df)

 }

 ##' Equal daylight layout method for unrooted trees.

 ##'

 ##' #' @title
 ##' @param tree phylo object
 ##' @param branch.length set to 'none' for edge length of 1. Otherwise the phylogenetic tree edge length is used.
 ##' @return tree as data.frame with equal angle layout.

 ##' @references

 ##' The following aglorithm aims to implement the vague description of the "Equal-daylight Algorithm"
 ##' in "Inferring Phylogenies" pp 582-584 by Joseph Felsenstein.

 ##'

 ##' ```
 ##' Leafs are subtrees with no children
 ##' Initialise tree using equal angle algorithm
 ##' tree_df = equal_angle(tree)
 ##'
 ##' nodes = get list of nodes in tree_df breadth-first
 ##' nodes = remove tip nodes.

 ##'

 ##' ```
 layoutDaylight <- function( tree, branch.length ){
 

     ## How to set optimal

     MAX_COUNT <- 5
     MINIMUM_AVERAGE_ANGLE_CHANGE <- 0.05

     ## Initialize tree.
     tree_df <- layoutEqualAngle(tree, branch.length)

     ## nodes = get list of nodes in tree_df
     ## Get list of node id's.
     ## nodes <- getNodes_by_postorder(tree)

     ## nodes <- getSubtree.df(tree_df, root)

     ## Get list of internal nodes
     ## nodes <- tree_df[tree_df$IsTip != TRUE]$nodes

     nodes <- getNodesBreadthFirst.df(tree_df)
     ## select only internal nodes
     internal_nodes <- tree_df[!tree_df$isTip,]$node
     ## Remove tips from nodes list, but keeping order.
     nodes <- intersect(nodes, internal_nodes)

     i <- 1
     ave_change <- 1.0
     while( i <= MAX_COUNT & ave_change > MINIMUM_AVERAGE_ANGLE_CHANGE ){

         message('Iteration: ', i)

         ## Reset max_change after iterating over tree.
         total_max <- 0.0

         ## for node in nodes {
         for( j in seq_along(nodes)){
             currentNode_id <- nodes[j]

             result <- applyLayoutDaylight(tree_df, currentNode_id)
             tree_df <- result$tree
             total_max <- total_max + result$max_change

         }

         # Calculate the running average of angle changes.
         ave_change <- total_max / length(nodes) * length(i)

         cat('Average angle change [',i,']', ave_change,'\n')

         i <- i + 1
     }

     return(tree_df)

 }

 ##' Apply the daylight alorithm to adjust the spacing between the subtrees and tips of the
 ##' specified node.

 ##'

 ##' @title applyLayoutDaylight
 ##' @param df tree data.frame
 ##' @param node_id is id of the node from which daylight is measured to the other subtrees.

 ##' @return list with tree data.frame with updated layout using daylight algorithm and max_change angle.
 ##
 ##
 ## ```
 ## for node in nodes {
 ##   if node is a leaf {
 ##     next
 ##   }
 ##
 ##   subtrees = get subtrees of node
 ##
 ##   for i-th subtree in subtrees {
 ##     [end, start] = get left and right angles of tree from node id.
 ##     angle_list[i, 'left'] = end
 ##     angle_list[i, 'beta'] = start - end  # subtree arc angle
 ##     angle_list[i, 'index'] = i-th # index of subtree/leaf
 ##   }
 ##
 ##   sort angle_list by 'left' column in ascending order.
 ##
 ##   D = 360 - sum( angle_list['beta'] ) # total daylight angle
 ##   d = D / |subtrees| # equal daylight angle.
 ##
 ##   new_L = left angle of first subtree.
 ##
 ##   for n-th row in angle_list{
 ##     # Calculate angle to rotate subtree/leaf to create correct daylight angle.
 ##     new_left_angle = new_left_angle + d + angle_list[n, 'beta']
 ##     Calculate the difference between the old and new left angles.
 ##     adjust_angle = new_left_angle - angle_list[n, 'left']
 ##
 ##     index = angle_list['index']
 ##     rotate subtree[index] wrt n-th node by adjust_angle
 ##     }
 ##   }
 ## }
 ## ```

 applyLayoutDaylight <- function(df, node_id){

   max_change <- 0.0

   # Get lists of node ids for each subtree, including  rest of unrooted tree.
   subtrees <- getSubtreeUnrooted.df(df, node_id)
   angle_list <- data.frame(left=numeric(0), beta=numeric(0), subtree_id=integer(0) )
 
   # Return tree if only 2 or less subtrees to adjust.
   if(length(subtrees) <= 2){
     return( list(tree = df, max_change = max_change) )
   }

   # Find start and end angles for each subtree.
   #   subtrees = get subtrees of node
   #   for i-th subtree in subtrees {
   for (i in seq_along(subtrees) ) {
     subtree <- subtrees[[i]]
     # [end, start] = get start and end angles of tree.

     angles <- getTreeArcAngles(df, node_id, subtree)
     angle_list[ i, 'subtree_id'] <- i
     angle_list[ i, 'left'] <- angles['left']
     angle_list[ i, 'beta'] <- angles['left'] - angles['right'] # subtree arc angle
     # If subtree arc angle is -ve, then + 2 (360).
     if(angle_list[ i, 'beta'] < 0 ){
       angle_list[ i, 'beta'] <- angle_list[ i, 'beta'] + 2

     }

   }
   #   sort angle_list by 'left angle' column in ascending order.
   angle_list <- angle_list[with(angle_list, order(left)), ]
   #   D = 360 - sum( angle_list['beta'] ) # total day
   #   d = D / |subtrees| # equal daylight angle.

   total_daylight <- 2 - colSums(angle_list['beta'])

   d <- total_daylight / length(subtrees)

   # Initialise new left-angle as first subtree left-angle.
   new_left_angle <- angle_list[1, 'left']
 
   # Adjust angles of subtrees and tips connected to current node.
   # for n-th row in angle_list{
   # Skip the first subtree as it is not adjusted.
   for (i in 2:nrow(angle_list) ) {
     # Calculate angle to rotate subtree/leaf to create correct daylight angle.
     new_left_angle <- new_left_angle + d + angle_list[i, 'beta']
     # Calculate the difference between the old and new left angles.
     adjust_angle <- new_left_angle - angle_list[i, 'left']

     max_change <- max(max_change, abs(adjust_angle))
     #cat('Adjust angle:', abs(adjust_angle), ' Max change:', max_change ,'\n')

     # rotate subtree[index] wrt current node
     subtree_id <- angle_list[i, 'subtree_id']
     subtree_nodes <- subtrees[[subtree_id]]$subtree
     # update tree_df for all subtrees with rotated points.
     df <- rotateTreePoints.df(df, node_id, subtree_nodes, adjust_angle)
   }

 
   return( list(tree = df, max_change = max_change) )
 

 }
 
 

 ##' Find the right (clockwise rotation, angle from +ve x-axis to furthest subtree nodes) and

 ##' left (anti-clockwise angle from +ve x-axis to subtree) Returning arc angle in [0, 2] (0 to 360) domain.

 ##'

 ##' @title getTreeArcAngles
 ##' @param df tree data.frame
 ##' @param origin_id node id from which to calculate left and right hand angles of subtree.

 ##' @param subtree named list of root id of subtree (node) and list of node ids for given subtree (subtree).

 ##' @return named list with right and left angles in range [0,2] i.e 1 = 180 degrees, 1.5 = 270 degrees.
 getTreeArcAngles <- function(df, origin_id, subtree) {
   # Initialise variables
   theta_child <- 0.0
   subtree_root_id <- subtree$node
   subtree_node_ids <- subtree$subtree

   # Initialise angle from origin node to parent node.
   # If subtree_root_id is child of origin_id
   if( any(subtree_root_id == getChild.df(df, origin_id)) ){

     # get angle from original node to parent of subtree.
     theta_left <- getNodeAngle.df(df, origin_id, subtree_root_id)
     theta_right <- theta_left
   }else if( subtree_root_id == origin_id){
     # Special case.
     # get angle from parent of subtree to children
     children_ids <- getChild.df(df, subtree_root_id)

     if(length(children_ids) == 2){
       # get angles from parent to it's two children.
       theta1 <- getNodeAngle.df(df, origin_id, children_ids[1])
       theta2 <- getNodeAngle.df(df, origin_id, children_ids[2])

       delta <- theta1 - theta2

       # correct delta for points crossing 180/-180 quadrant.
       if(delta > 1){
         delta_adj = delta - 2
       }else if(delta < -1){
         delta_adj = delta + 2
       }else{
         delta_adj <- delta
       }

       if(delta_adj >= 0){
         theta_left = theta1
         theta_right = theta2
       }else if(delta_adj < 0){
         theta_left = theta2
         theta_right = theta1
       }
     }else{
       # subtree only has one child node.
       theta_left <- getNodeAngle.df(df, origin_id, children_ids[1])
       theta_right <- theta_left
     }

   }else{
     # get the real root of df tree to initialise left and right angles.

     tree_root <- getRoot.df(df)
     if( !is.na(tree_root) & is.numeric(tree_root) ){
       theta_left <- getNodeAngle.df(df, origin_id, tree_root)
       theta_right <- theta_left
     }else{
       print('ERROR: no root found!')
       theta_left <- NA
     }
 

   }

   # no parent angle found.
   if (is.na(theta_left) ){
     return(0)
   }

   # create vector with named columns
   # left-hand and right-hand angles between origin node and the extremities of the tree nodes.

   arc <- c('left' = theta_left, 'right' = theta_right)

   # Subtree has to have 1 or more nodes to compare.
   if (length(subtree_node_ids) == 0 ){
     return(0)

   }
 

 
   # Remove tips from nodes list, but keeping order.
   # internal_nodes <- df[!df$isTip,]$node
   # subtree_node_ids <- intersect(subtree_node_ids, internal_nodes)

   # Calculate the angle from the origin node to each child node.
   # Moving from parent to children in depth-first traversal.
   for( i in seq_along(subtree_node_ids) ){
     parent_id <- subtree_node_ids[i]
     # Get angle from origin node to parent node.

     # Skip if parent_id is a tip or parent and child node are the same.
     if(origin_id == parent_id | isTip.df(df, parent_id) ){

       next

     }

 
     theta_parent <- getNodeAngle.df(df, origin_id, parent_id)

     children_ids <- getChild.df(df, parent_id)

     for( j in seq_along(children_ids)){
       #delta_x <- df[subtree_node_id, 'x'] - df[origin_id, 'x']
       #delta_y <- df[subtree_node_id, 'y'] - df[origin_id, 'y']
       #angles[i] <- atan2(delta_y, delta_x) / pi
       child_id <- children_ids[j]
       # Skip if child is parent node of subtree.
       if( child_id == origin_id ){
         next
       }

       theta_child <- getNodeAngle.df(df, origin_id, child_id)
 
       # Skip if child node is already inside arc.
       # if left < right angle (arc crosses 180/-180 quadrant) and child node is not inside arc of tree.
       # OR if left > right angle (arc crosses 0/360 quadrant) and child node is inside gap
       if ( (arc['left'] < arc['right'] & !( theta_child > arc['left'] & theta_child < arc['right'])) |
         (arc['left'] > arc['right'] & ( theta_child < arc['left'] & theta_child > arc['right'])) ){
         # child node inside arc.
         next
       }

       delta <- theta_child - theta_parent
       delta_adj <- delta
       # Correct the delta if parent and child angles cross the 180/-180 half of circle.
       # If delta > 180
       if( delta > 1){ # Edge between parent and child cross upper and lower quadrants of cirlce on 180/-180 side.
         delta_adj <- delta - 2 # delta' = delta - 360
       # If delta < -180
       }else if( delta < -1){ # Edge between parent and child cross upper and lower quadrants of cirlce
         delta_adj <- delta + 2 # delta' = delta - 360
       }

       theta_child_adj <- theta_child
 
       # If angle change from parent to node is positive (anti-clockwise), check left angle
       if(delta_adj > 0){

         # If child/parent edges cross the -180/180 quadrant (angle between them is > 180),

         # check if right angle and child angle are different signs and adjust if needed.
         if( abs(delta) > 1){
           if( arc['left'] > 0 & theta_child < 0){
             theta_child_adj <- theta_child + 2
           }else if (arc['left'] < 0 & theta_child > 0){
             theta_child_adj <- theta_child - 2
           }
         }

           # check if left angle of arc is less than angle of child. Update if true.
         if( arc['left'] < theta_child_adj ){
           arc['left'] <- theta_child

         }

       # If angle change from parent to node is negative (clockwise), check right angle

       }else if(delta_adj < 0){

         # If child/parent edges cross the -180/180 quadrant (angle between them is > 180),

         # check if right angle and child angle are different signs and adjust if needed.
         if( abs(delta) > 1){
           # Else change in angle from parent to child is negative, then adjust child angle if right angle is a different sign.
           if( arc['right'] > 0 & theta_child < 0){
             theta_child_adj <- theta_child + 2
           }else if (arc['right'] < 0 & theta_child > 0){
             theta_child_adj <- theta_child - 2
           }
         }
         # check if right angle of arc is greater than angle of child. Update if true.
         if( arc['right'] > theta_child_adj  ){
           arc['right'] <- theta_child

         }
 

       }
     }

   }

   # Convert arc angles of [1, -1] to [2,0] domain.
   arc[arc<0] <- arc[arc<0] + 2
   return(arc)

 }

 ##' Rotate the points in a tree data.frame around a pivot node by the angle specified.

 ##'

 ##' @title rotateTreePoints.data.fram
 ##' @param df tree data.frame
 ##' @param pivot_node is the id of the pivot node.
 ##' @param nodes list of node numbers that are to be rotated by angle around the pivot_node
 ##' @param angle in range [0,2], ie degrees/180, radians/pi
 ##' @return updated tree data.frame with points rotated by angle
 rotateTreePoints.df <- function(df, pivot_node, nodes, angle){
   # Rotate nodes around pivot_node.
   # x' = cos(angle)*delta_x - sin(angle)*delta_y + delta_x
   # y' = sin(angle)*delta_x + cos(angle)*delta_y + delta_y

   cospitheta <- cospi(angle)
   sinpitheta <- sinpi(angle)
   for(node in nodes){
     # Update (x,y) of node
     delta_x <- df[node, 'x'] - df[pivot_node, 'x']
     delta_y <- df[node, 'y'] - df[pivot_node, 'y']
     df[node, 'x'] <- cospitheta * delta_x - sinpitheta * delta_y + df[pivot_node, 'x']
     df[node, 'y'] <- sinpitheta * delta_x + cospitheta * delta_y + df[pivot_node, 'y']

   }

   # Now update tip labels of rotated tree.
   # angle is in range [0, 360]

   for(node in nodes){

     # Update label angle of tipnode if not root node.
     if( isTip.df(df, node) ){
       # get parent
       parent_id <- getParent.df(df, node)
       # if 'node' is not root, then update label angle.
       if( parent_id != 0 ){
         theta_parent_child <- getNodeAngle.df(df, parent_id, node)
         if(!is.na(theta_parent_child)){
           # Update tip label angle, that is parallel to edge.
           #df[node, 'angle'] <- -90 - 180 * theta_parent_child * sign(theta_parent_child - 1)
           if(theta_parent_child > 0 ){

             df[node, 'angle'] <- 180 * theta_parent_child

           }else if(theta_parent_child < 0 ){
             df[node, 'angle'] <- 180 * ( theta_parent_child + 2 )
           }

         }

       }
     }
   }

   return(df)
 }

 ##' Get the angle between the two nodes specified.

 ##'

 ##' @title getNodeAngle.df
 ##' @param df tree data.frame
 ##' @param origin_node_id origin node id number
 ##' @param node_id end node id number
 ##' @return angle in range [-1, 1], i.e. degrees/180, radians/pi
 getNodeAngle.df <- function(df, origin_node_id, node_id){

   if( (origin_node_id != node_id) & any(origin_node_id %in% df$node) & any(node_id %in% df$node) ){

     delta_x <- df[node_id, 'x'] - df[origin_node_id, 'x']
     delta_y <- df[node_id, 'y'] - df[origin_node_id, 'y']
     angle <- atan2(delta_y, delta_x) / pi
     return( angle )
   }else{
     return(NA)
   }
 }

 euc.dist <- function(x1, x2) sqrt(sum((x1 - x2) ^ 2))
 

 ## Get the distances from the node to all other nodes in data.frame (including itself if in df)

 getNodeEuclDistances <- function(df, node){
   # https://stackoverflow.com/questions/24746892/how-to-calculate-euclidian-distance-between-two-points-defined-by-matrix-contain#24747155
   dist <- NULL
   for(i in 1:nrow(df)) dist[i] <- euc.dist(df[df$node==node, c('x', 'y')], df[i, c('x', 'y')])
   return(dist)
 }

 ##' Get all children of node from tree, including start_node.

 ##'

 ##' @title getSubtree
 ##' @param tree ape phylo tree object
 ##' @param node is the tree node id from which the tree is derived.
 ##' @return list of all child node id's from starting node.
 getSubtree <- function(tree, node){

   subtree <- c(node)
   i <- 1
   while( i <= length(subtree)){
     subtree <- c(subtree, getChild(tree, subtree[i]))
     # remove any '0' root nodes
     subtree <- subtree[subtree != 0]
     i <- i + 1
   }
   return(subtree)

 }
 

 ##' Get all children of node from df tree using breath-first.
 ##'

 ##' @title getSubtree.df

 ##' @param df tree data.frame
 ##' @param node id of starting node.
 ##' @return list of all child node id's from starting node.

 getSubtree.df <- function(df, node){

   subtree <- c(node)

   subtree <- subtree[subtree != 0]

   i <- 1
   while( i <= length(subtree)){
     subtree <- c(subtree, getChild.df(df, subtree[i]))
     # remove any '0' root nodes
     subtree <- subtree[subtree != 0]
     i <- i + 1
   }
   return(subtree)
 }
 

 ##' Get all subtrees of specified node. This includes all ancestors and relatives of node and

 ##' return named list of subtrees.

 ##'

 ##' @title getSubtreeUnrooted

 ##' @param tree ape phylo tree object
 ##' @param node is the tree node id from which the subtrees are derived.
 ##' @return named list of subtrees with the root id of subtree and list of node id's making up subtree.

 getSubtreeUnrooted <- function(tree, node){

   # if node leaf, return nothing.
   if( isTip(tree, node) ){
     # return NA
     return(NA)
   }

   subtrees <- list()

   # get subtree for each child node.
   children_ids <- getChild(tree, node)

   remaining_nodes <- getNodes_by_postorder(tree)
   # Remove current node from remaining_nodes list.
   remaining_nodes <- setdiff(remaining_nodes, node)

   for( child in children_ids ){
     # Append subtree nodes to list if not 0 (root).
     subtree <- getSubtree(tree, child)
     subtrees[[length(subtrees)+1]] <- list( node = child, subtree = subtree)
     # remove subtree nodes from remaining nodes.
     remaining_nodes <- setdiff(remaining_nodes, as.integer(unlist(subtrees[[length(subtrees)]]['subtree']) ))
   }

   # The remaining nodes that are not found in the child subtrees are the remaining subtree nodes.
   # ie, parent node and all other nodes. We don't care how they are connect, just their ids.
   parent_id <- getParent(tree, node)
   # If node is not root, add remainder of tree nodes as subtree.
   if( parent_id != 0 & length(remaining_nodes) >= 1){
     subtrees[[length(subtrees)+1]] <- list( node = parent_id, subtree = remaining_nodes)
   }

   return(subtrees)
 }
 
 
 ##' Get all subtrees of node, as well as remaining branches of parent (ie, rest of tree structure as subtree)
 ##' return named list of subtrees with list name as starting node id.

 ##' @title getSubtreeUnrooted

 ##' @param df tree data.frame
 ##' @param node is the tree node id from which the subtrees are derived.
 ##' @return named list of subtrees with the root id of subtree and list of node id's making up subtree.
 getSubtreeUnrooted.df <- function(df, node){
   # if node leaf, return nothing.
   if( isTip.df(df, node) ){
     return(NA)
   }

   subtrees <- list()

   # get subtree for each child node.
   children_ids <- getChild.df(df, node)

   # remaining_nodes <- getNodes_by_postorder(tree)
   remaining_nodes <- df$node

   # Remove current node from remaining_nodes list.
   remaining_nodes <- setdiff(remaining_nodes, node)

   for( child in children_ids ){

     subtree <- getSubtree.df(df, child)

     # Append subtree nodes to list if more than 1 node in subtree (i.e. not a tip)
     #if(length(subtree) >= 2){
       subtrees[[length(subtrees)+1]] <- list( node = child, subtree = subtree)
       # remove subtree nodes from remaining nodes.
       remaining_nodes <- setdiff(remaining_nodes, as.integer(unlist(subtrees[[length(subtrees)]]['subtree']) ))
     #}else{
       # remove remaining nodes
     #  remaining_nodes <- setdiff(remaining_nodes, subtree)
     #}
   }

   # The remaining nodes that are not found in the child subtrees are the remaining subtree nodes.
   # ie, parent node and all other nodes. We don't care how they are connected, just their id.
   parent_id <- getParent.df(df, node)
   # If node is not root.
   if( parent_id != 0 & length(remaining_nodes) >= 1){
     subtrees[[length(subtrees)+1]] <- list( node = parent_id, subtree = remaining_nodes)
   }

   return(subtrees)
 }
 
 
 getRoot.df <- function(df, node){

   root <- which(is.na(df$parent))

   # Check if root was found.
   if(length(root) == 0){
     # Alternatively, root can self reference, eg node = 10, parent = 10
     root <- unlist(apply(df, 1, function(x){ if(x['node'] == x['parent']){ x['node'] } }))
   }

   return(root)
 }
 
 ##' Get parent node id of child node.
 ##'
 ##' @title getParent.df
 ##' @param df tree data.frame
 ##' @param node is the node id of child in tree.
 ##' @return integer node id of parent

 getParent.df <- function(df, node) {
     i <- which(df$node == node)

     parent_id <- df$parent[i]
     if (parent_id == node | is.na(parent_id)) {

         ## root node
         return(0)
     }

     return(parent_id)

 }
 
 getAncestor.df <- function(df, node) {
     anc <- getParent.df(df, node)
     anc <- anc[anc != 0]
     if (length(anc) == 0) {

         # stop("selected node is root...")
       return(0)

     }
     i <- 1
     while(i<= length(anc)) {
         anc <- c(anc, getParent.df(df, anc[i]))
         anc <- anc[anc != 0]
         i <- i+1
     }
     return(anc)
 }
 
 

 ##' Get list of child node id numbers of parent node
 ##'
 ##' @title getChild.df
 ##' @param df tree data.frame
 ##' @param node is the node id of child in tree.
 ##' @return list of child node ids of parent

 getChild.df <- function(df, node) {
     i <- which(df$parent == node)
     if (length(i) == 0) {

         return(0) # it has no children, hence tip node.

     }

     res <- df$node[i]

     res <- res[res != node] ## node may root
     return(res)
 }
 
 get.offspring.df <- function(df, node) {
     sp <- getChild.df(df, node)

     sp <- sp[sp != 0] # Remove root node.

     if (length(sp) == 0) {

         #stop("input node is a tip...")
       return(0)

     }
 
     i <- 1
     while(i <= length(sp)) {
         sp <- c(sp, getChild.df(df, sp[i]))
         sp <- sp[sp != 0]
         i <- i + 1
     }
     return(sp)
 }
 
 
 ##' extract offspring tips
 ##'
 ##'
 ##' @title get.offspring.tip
 ##' @param tr tree
 ##' @param node node
 ##' @return tip label
 ##' @author ygc
 ##' @importFrom ape extract.clade
 ##' @export
 get.offspring.tip <- function(tr, node) {
     if ( ! node %in% tr$edge[,1]) {
         ## return itself
         return(tr$tip.label[node])
     }
     clade <- extract.clade(tr, node)
     clade$tip.label
 }
 
 

 ## ##' calculate total number of nodes
 ## ##'
 ## ##'
 ## ##' @title getNodeNum
 ## ##' @param tr phylo object
 ## ##' @return number
 ## ##' @author Guangchuang Yu
 ## ##' @export
 ## getNodeNum <- function(tr) {
 ##     Ntip <- length(tr[["tip.label"]])
 ##     Nnode <- tr[["Nnode"]]
 ##     ## total nodes
 ##     N <- Ntip + Nnode
 ##     return(N)
 ## }

 getParent <- function(tr, node) {
     if ( node == getRoot(tr) )
         return(0)
     edge <- tr[["edge"]]
     parent <- edge[,1]
     child <- edge[,2]
     res <- parent[child == node]
     if (length(res) == 0) {
         stop("cannot found parent node...")
     }
     if (length(res) > 1) {
         stop("multiple parent found...")
     }
     return(res)
 }
 
 getChild <- function(tr, node) {

     # Get edge matrix from phylo object.

     edge <- tr[["edge"]]

     # Select all rows that match "node".

     res <- edge[edge[,1] == node, 2]
     ## if (length(res) == 0) {
     ##     ## is a tip
     ##     return(NA)
     ## }
     return(res)
 }
 
 getSibling <- function(tr, node) {
     root <- getRoot(tr)
     if (node == root) {
         return(NA)
     }
 
     parent <- getParent(tr, node)
     child <- getChild(tr, parent)
     sib <- child[child != node]
     return(sib)
 }
 
 
 getAncestor <- function(tr, node) {
     root <- getRoot(tr)
     if (node == root) {
         return(NA)
     }
     parent <- getParent(tr, node)
     res <- parent
     while(parent != root) {
         parent <- getParent(tr, parent)
         res <- c(res, parent)
     }
     return(res)
 }
 
 isRoot <- function(tr, node) {
     getRoot(tr) == node
 }
 

 isTip <- function(tr, node) {
   children_ids <- getChild(tr, node)

   #length(children_ids) == 0 ## getChild returns 0 if nothing found.
   if( length(children_ids) == 0 | any(children_ids == 0) ){
     return(TRUE)
   }
   return(FALSE)

 }
 
 isTip.df <- function(df, node) {

   # df may not have the isTip structure.
   # return(df[node, 'isTip'])
   # Tip has no children.
   children_ids <- getChild.df(df, node)
   if( length(children_ids) == 0 | any(children_ids == 0) ){
     return(TRUE)
   }
   return(FALSE)

 }
 
 

 getNodeName <- function(tr) {
     if (is.null(tr$node.label)) {
         n <- length(tr$tip.label)
         nl <- (n + 1):(2 * n - 2)
         nl <- as.character(nl)
     }
     else {
         nl <- tr$node.label
     }
     nodeName <- c(tr$tip.label, nl)
     return(nodeName)
 }
 

 ## ##' get the root number
 ## ##'
 ## ##'
 ## ##' @title getRoot
 ## ##' @param tr phylo object
 ## ##' @return root number
 ## ##' @export
 ## ##' @author Guangchuang Yu
 ## getRoot <- function(tr) {
 ##     edge <- tr[["edge"]]
 ##     ## 1st col is parent,
 ##     ## 2nd col is child,
 ##     if (!is.null(attr(tr, "order")) && attr(tr, "order") == "postorder")
 ##         return(edge[nrow(edge), 1])
 
 ##     parent <- unique(edge[,1])
 ##     child <- unique(edge[,2])
 ##     ## the node that has no parent should be the root
 ##     root <- parent[ ! parent %in% child ]
 ##     if (length(root) > 1) {
 ##         stop("multiple roots founded...")
 ##     }
 ##     return(root)
 ## }

 get.trunk <- function(tr) {
     root <- getRoot(tr)
     path_length <- sapply(1:(root-1), function(x) get.path_length(tr, root, x))
     i <- which.max(path_length)
     return(get.path(tr, root, i))
 }
 
 ##' path from start node to end node
 ##'
 ##'
 ##' @title get.path
 ##' @param phylo phylo object
 ##' @param from start node
 ##' @param to end node
 ##' @return node vectot
 ##' @export
 ##' @author Guangchuang Yu
 get.path <- function(phylo, from, to) {
     anc_from <- getAncestor(phylo, from)
     anc_from <- c(from, anc_from)
     anc_to <- getAncestor(phylo, to)
     anc_to <- c(to, anc_to)
     mrca <- intersect(anc_from, anc_to)[1]
 
     i <- which(anc_from == mrca)
     j <- which(anc_to == mrca)
 
     path <- c(anc_from[1:i], rev(anc_to[1:(j-1)]))
     return(path)
 }
 
 get.path_length <- function(phylo, from, to, weight=NULL) {
     path <- get.path(phylo, from, to)
     if (is.null(weight)) {
         return(length(path)-1)
     }
 
     df <- fortify(phylo)
     if ( ! (weight %in% colnames(df))) {
         stop("weight should be one of numerical attributes of the tree...")
     }
 
     res <- 0
 
     get_edge_index <- function(df, from, to) {
         which((df[,1] == from | df[,2] == from) &
                   (df[,1] == to | df[,2] == to))
     }
 
     for(i in 1:(length(path)-1)) {
         ee <- get_edge_index(df, path[i], path[i+1])
         res <- res + df[ee, weight]
     }
 
     return(res)
 }
 
 getNodes_by_postorder <- function(tree) {

   tree <- reorder.phylo(tree, "postorder")

     unique(rev(as.vector(t(tree$edge[,c(2,1)]))))
 }
 

 #getNodes_by_postorder.df <- function(df) {
     #tree <- reorder.phylo(tree, "postorder")
     #unique(rev(as.vector(t(tree$edge[,c(2,1)]))))
 #}
 
 ##' Get the nodes of tree from root in breadth-first order.

 ##'

 ##' @title getNodesBreadthFirst.df
 ##' @param df tree data.frame
 ##' @return list of node id's in breadth-first order.
 getNodesBreadthFirst.df <- function(df){
 
   root <- getRoot.df(df)
   if(isTip.df(df, root)){
     return(root)
   }
 
   tree_size <- nrow(df)
   # initialise list of nodes
   res <- root

   i <- 1
   while(length(res) < tree_size){
     parent <- res[i]
     i <- i + 1

     # Skip if parent is a tip.
     if(isTip.df(df, parent)){
       next
     }

     # get children of current parent.
     children <- getChild.df(df,parent)
 
     # add children to result
     res <- c(res, children)

   }

   return(res)

 }
 
 

 getXcoord2 <- function(x, root, parent, child, len, start=0, rev=FALSE) {
     x[root] <- start
     x[-root] <- NA  ## only root is set to start, by default 0
 
     currentNode <- root
     direction <- 1
     if (rev == TRUE) {
         direction <- -1
     }

     while(anyNA(x)) {
         idx <- which(parent %in% currentNode)
         newNode <- child[idx]
         x[newNode] <- x[parent[idx]]+len[idx] * direction
         currentNode <- newNode
     }
 
     return(x)
 }
 
 getXcoord_no_length <- function(tr) {
     edge <- tr$edge
     parent <- edge[,1]
     child <- edge[,2]
     root <- getRoot(tr)
 
     len <- tr$edge.length
 
     N <- getNodeNum(tr)
     x <- numeric(N)
     ntip <- Ntip(tr)
     currentNode <- 1:ntip
     x[-currentNode] <- NA
 
     cl <- split(child, parent)
     child_list <- list()
     child_list[as.numeric(names(cl))] <- cl
 
     while(anyNA(x)) {
         idx <- match(currentNode, child)
         pNode <- parent[idx]
         ## child number table
         p1 <- table(parent[parent %in% pNode])
         p2 <- table(pNode)
         np <- names(p2)
         i <- p1[np] == p2
         newNode <- as.numeric(np[i])
 
         exclude <- rep(NA, max(child))
         for (j in newNode) {
             x[j] <- min(x[child_list[[j]]]) - 1
             exclude[child_list[[j]]] <- child_list[[j]]
         }
         exclude <- exclude[!is.na(exclude)]
 
         ## currentNode %<>% `[`(!(. %in% exclude))
         ## currentNode %<>% c(., newNode) %>% unique
         currentNode <- currentNode[!currentNode %in% exclude]
         currentNode <- unique(c(currentNode, newNode))
 
     }
     x <- x - min(x)
     return(x)
 }
 
 
 getXcoord <- function(tr) {
     edge <- tr$edge
     parent <- edge[,1]
     child <- edge[,2]
     root <- getRoot(tr)
 
     len <- tr$edge.length
 
     N <- getNodeNum(tr)
     x <- numeric(N)
     x <- getXcoord2(x, root, parent, child, len)
     return(x)
 }
 
 getXYcoord_slanted <- function(tr) {
 
     edge <- tr$edge
     parent <- edge[,1]
     child <- edge[,2]
     root <- getRoot(tr)
 
     N <- getNodeNum(tr)
     len <- tr$edge.length
     y <- getYcoord(tr, step=min(len)/2)
 
     len <- sqrt(len^2 - (y[parent]-y[child])^2)
     x <- numeric(N)
     x <- getXcoord2(x, root, parent, child, len)
     res <- data.frame(x=x, y=y)
     return(res)
 }
 
 
 ## @importFrom magrittr %>%
 ##' @importFrom magrittr equals
 getYcoord <- function(tr, step=1) {
     Ntip <- length(tr[["tip.label"]])
     N <- getNodeNum(tr)
 
     edge <- tr[["edge"]]
     parent <- edge[,1]
     child <- edge[,2]
 
     cl <- split(child, parent)
     child_list <- list()
     child_list[as.numeric(names(cl))] <- cl
 
     y <- numeric(N)
     tip.idx <- child[child <= Ntip]
     y[tip.idx] <- 1:Ntip * step
     y[-tip.idx] <- NA
 

     ## use lookup table
     pvec <- integer(max(tr$edge))
     pvec[child] = parent
 

     currentNode <- 1:Ntip
     while(anyNA(y)) {

         ## pNode <- unique(parent[child %in% currentNode])
         pNode <- unique(pvec[currentNode])
 

         ## piping of magrittr is slower than nested function call.
         ## pipeR is fastest, may consider to use pipeR
         ##
         ## child %in% currentNode %>% which %>% parent[.] %>% unique
         ## idx <- sapply(pNode, function(i) all(child[parent == i] %in% currentNode))
         idx <- sapply(pNode, function(i) all(child_list[[i]] %in% currentNode))
         newNode <- pNode[idx]
 
         y[newNode] <- sapply(newNode, function(i) {
             mean(y[child_list[[i]]], na.rm=TRUE)
             ##child[parent == i] %>% y[.] %>% mean(na.rm=TRUE)
         })
 
         currentNode <- c(currentNode[!currentNode %in% unlist(child_list[newNode])], newNode)
         ## currentNode <- c(currentNode[!currentNode %in% child[parent %in% newNode]], newNode)
         ## parent %in% newNode %>% child[.] %>%
         ##     `%in%`(currentNode, .) %>% `!` %>%
         ##         currentNode[.] %>% c(., newNode)
     }
 
     return(y)
 }
 
 
 getYcoord_scale <- function(tr, df, yscale) {
 
     N <- getNodeNum(tr)
     y <- numeric(N)
 
     root <- getRoot(tr)
     y[root] <- 0
     y[-root] <- NA
 
     edge <- tr$edge
     parent <- edge[,1]
     child <- edge[,2]
 
     currentNodes <- root
     while(anyNA(y)) {
         newNodes <- c()
         for (currentNode in currentNodes) {
             idx <- which(parent %in% currentNode)
             newNode <- child[idx]
             direction <- -1
             for (i in seq_along(newNode)) {
                 y[newNode[i]] <- y[currentNode] + df[newNode[i], yscale] * direction
                 direction <- -1 * direction
             }
             newNodes <- c(newNodes, newNode)
         }
         currentNodes <- unique(newNodes)
     }
     if (min(y) < 0) {
         y <- y + abs(min(y))
     }
     return(y)
 }
 
 getYcoord_scale2 <- function(tr, df, yscale) {
     root <- getRoot(tr)
 
     pathLength <- sapply(1:length(tr$tip.label), function(i) {
         get.path_length(tr, i, root, yscale)
     })
 
     ordered_tip <- order(pathLength, decreasing = TRUE)
     ii <- 1
     ntip <- length(ordered_tip)
     while(ii < ntip) {
         sib <- getSibling(tr, ordered_tip[ii])
         if (length(sib) == 0) {
             ii <- ii + 1
             next
         }
         jj <- which(ordered_tip %in% sib)
         if (length(jj) == 0) {
             ii <- ii + 1
             next
         }
         sib <- ordered_tip[jj]
         ordered_tip <- ordered_tip[-jj]
         nn <- length(sib)
         if (ii < length(ordered_tip)) {
             ordered_tip <- c(ordered_tip[1:ii],sib, ordered_tip[(ii+1):length(ordered_tip)])
         } else {
             ordered_tip <- c(ordered_tip[1:ii],sib)
         }
 
         ii <- ii + nn + 1
     }
 
 
     long_branch <- getAncestor(tr, ordered_tip[1]) %>% rev
     long_branch <- c(long_branch, ordered_tip[1])
 
     N <- getNodeNum(tr)
     y <- numeric(N)
 
     y[root] <- 0
     y[-root] <- NA
 
     ## yy <- df[, yscale]
     ## yy[is.na(yy)] <- 0
 
     for (i in 2:length(long_branch)) {
         y[long_branch[i]] <- y[long_branch[i-1]] + df[long_branch[i], yscale]
     }
 
     parent <- df[, "parent"]
     child <- df[, "node"]
 
     currentNodes <- root
     while(anyNA(y)) {
         newNodes <- c()
         for (currentNode in currentNodes) {
             idx <- which(parent %in% currentNode)
             newNode <- child[idx]
             newNode <- c(newNode[! newNode %in% ordered_tip],
                          rev(ordered_tip[ordered_tip %in% newNode]))
             direction <- -1
             for (i in seq_along(newNode)) {
                 if (is.na(y[newNode[i]])) {
                     y[newNode[i]] <- y[currentNode] + df[newNode[i], yscale] * direction
                     direction <- -1 * direction
                 }
             }
             newNodes <- c(newNodes, newNode)
         }
         currentNodes <- unique(newNodes)
     }
     if (min(y) < 0) {
         y <- y + abs(min(y))
     }
     return(y)
 }
 
 getYcoord_scale_numeric <- function(tr, df, yscale, ...) {
     df <- .assign_parent_status(tr, df, yscale)
     df <- .assign_child_status(tr, df, yscale)
 
     y <- df[, yscale]
 
     if (anyNA(y)) {
         warning("NA found in y scale mapping, all were setting to 0")
         y[is.na(y)] <- 0
     }
 
     return(y)
 }
 
 .assign_parent_status <- function(tr, df, variable) {

     yy <- df[[variable]]

     na.idx <- which(is.na(yy))
     if (length(na.idx) > 0) {
         tree <- get.tree(tr)
         nodes <- getNodes_by_postorder(tree)
         for (curNode in nodes) {
             children <- getChild(tree, curNode)
             if (length(children) == 0) {
                 next
             }
             idx <- which(is.na(yy[children]))
             if (length(idx) > 0) {
                 yy[children[idx]] <- yy[curNode]
             }
         }
     }
     df[, variable] <- yy
     return(df)
 }
 
 .assign_child_status <- function(tr, df, variable, yscale_mapping=NULL) {

     yy <- df[[variable]]

     if (!is.null(yscale_mapping)) {
         yy <- yscale_mapping[yy]
     }
 
     na.idx <- which(is.na(yy))
     if (length(na.idx) > 0) {
         tree <- get.tree(tr)
         nodes <- rev(getNodes_by_postorder(tree))
         for (curNode in nodes) {
             parent <- getParent(tree, curNode)
             if (parent == 0) { ## already reach root
                 next
             }
             idx <- which(is.na(yy[parent]))
             if (length(idx) > 0) {
                 child <- getChild(tree, parent)
                 yy[parent[idx]] <- mean(yy[child], na.rm=TRUE)
             }
         }
     }
     df[, variable] <- yy
     return(df)
 }
 
 
 getYcoord_scale_category <- function(tr, df, yscale, yscale_mapping=NULL, ...) {
     if (is.null(yscale_mapping)) {
         stop("yscale is category variable, user should provide yscale_mapping,
              which is a named vector, to convert yscale to numberical values...")
     }
     if (! is(yscale_mapping, "numeric") ||
         is.null(names(yscale_mapping))) {
         stop("yscale_mapping should be a named numeric vector...")
     }
 
     if (yscale == "label") {

         yy <- df[[yscale]]

         ii <- which(is.na(yy))
         if (length(ii)) {
             df[ii, yscale] <- df[ii, "node"]
         }
     }
 
     ## assign to parent status is more prefer...
     df <- .assign_parent_status(tr, df, yscale)
     df <- .assign_child_status(tr, df, yscale, yscale_mapping)
 

     y <- df[[yscale]]

 
     if (anyNA(y)) {
         warning("NA found in y scale mapping, all were setting to 0")
         y[is.na(y)] <- 0
     }
     return(y)
 }
 
 
 add_angle_slanted <- function(res) {

     x <- res[["x"]]
     y <- res[["y"]]
     dy <- (y - y[match(res$parent, res$node)]) / diff(range(y))
     dx <- (x - x[match(res$parent, res$node)]) / diff(range(x))

     theta <- atan(dy/dx)
     theta[is.na(theta)] <- 0 ## root node
     res$angle <- theta/pi * 180

     branch.y <- (y[match(res$parent, res$node)] + y)/2

     idx <- is.na(branch.y)

     branch.y[idx] <- y[idx]

     res[, "branch.y"] <- branch.y
     return(res)
 }
 
 calculate_branch_mid <- function(res) {

     res$branch <- with(res, (x[match(parent, node)] + x)/2)

     if (!is.null(res$branch.length)) {
         res$branch.length[is.na(res$branch.length)] <- 0

     }
     res$branch[is.na(res$branch)] <- 0
     return(res)
 }
 
 
 set_branch_length <- function(tree_object, branch.length) {

     if (branch.length == "branch.length") {
         return(tree_object)
     } else if (branch.length == "none") {

         tree_object@phylo$edge.length <- NULL

         return(tree_object)

     }
 

     if (is(tree_object, "phylo")) {
         return(tree_object)

     }
 

     tree_anno <- get_tree_data(tree_object)

     if (is(tree_anno, "matrix"))
         tree_anno <- as.data.frame(tree_anno)
 

     phylo <- get.tree(tree_object)
 

     cn <- colnames(tree_anno)
     cn <- cn[!cn %in% c('node', 'parent')]
 
     length <- match.arg(branch.length, cn)
 

     if (all(is.na(as.numeric(tree_anno[[length]])))) {

         stop("branch.length should be numerical attributes...")
     }
 
     edge <- as.data.frame(phylo$edge)
     colnames(edge) <- c("parent", "node")
 
     dd <- merge(edge, tree_anno,

                 by  = "node",

                 all.x = TRUE)
     dd <- dd[match(edge$node, dd$node),]

     len <- unlist(dd[[length]])

     len <- as.numeric(len)
     len[is.na(len)] <- 0
 
     phylo$edge.length <- len
 

     tree_object@phylo <- phylo

     return(tree_object)

 }
 

 ## set_branch_length <- function(tree_object, branch.length) {
 ##     if (is(tree_object, "phylo4d")) {
 ##         phylo <- as.phylo.phylo4(tree_object)
 ##         d <- tree_object@data
 ##         tree_anno <- data.frame(node=rownames(d), d)
 ##     } else {
 ##         phylo <- get.tree(tree_object)
 ##     }
 
 ##     if (branch.length %in%  c("branch.length", "none")) {
 ##         return(phylo)
 ##     }
 
 ##     ## if (is(tree_object, "codeml")) {
 ##     ##     tree_anno <- tree_object@mlc@dNdS
 ##     ## } else
 
 ##     if (is(tree_object, "codeml_mlc")) {
 ##         tree_anno <- tree_object@dNdS
 ##     } else if (is(tree_object, "beast")) {
 ##         tree_anno <- tree_object@stats
 ##     }
 
 ##     if (has.extraInfo(tree_object)) {
 ##         tree_anno <- merge(tree_anno, tree_object@extraInfo, by.x="node", by.y="node")
 ##     }
 ##     cn <- colnames(tree_anno)
 ##     cn <- cn[!cn %in% c('node', 'parent')]
 
 ##     length <- match.arg(branch.length, cn)
 
 ##     if (all(is.na(as.numeric(tree_anno[, length])))) {
 ##         stop("branch.length should be numerical attributes...")
 ##     }
 
 ##     edge <- as.data.frame(phylo$edge)
 ##     colnames(edge) <- c("parent", "node")
 
 ##     dd <- merge(edge, tree_anno,
 ##                 by.x  = "node",
 ##                 by.y  = "node",
 ##                 all.x = TRUE)
 ##     dd <- dd[match(edge$node, dd$node),]
 ##     len <- unlist(dd[, length])
 ##     len <- as.numeric(len)
 ##     len[is.na(len)] <- 0
 
 ##     phylo$edge.length <- len
 
 ##     return(phylo)
 ## }

 re_assign_ycoord_df <- function(df, currentNode) {
     while(anyNA(df$y)) {
         pNode <- with(df, parent[match(currentNode, node)]) %>% unique
         idx <- sapply(pNode, function(i) with(df, all(node[parent == i & parent != node] %in% currentNode)))
         newNode <- pNode[idx]
         ## newNode <- newNode[is.na(df[match(newNode, df$node), "y"])]
 
         df[match(newNode, df$node), "y"] <- sapply(newNode, function(i) {
             with(df, mean(y[parent == i], na.rm = TRUE))
         })
         traced_node <- as.vector(sapply(newNode, function(i) with(df, node[parent == i])))
         currentNode <- c(currentNode[! currentNode %in% traced_node], newNode)
     }
     return(df)
 }
 
 

 ## ##' test whether input object is produced by ggtree function
 ## ##'
 ## ##'
 ## ##' @title is.ggtree
 ## ##' @param x object
 ## ##' @return TRUE or FALSE
 ## ##' @export
 ## ##' @author guangchuang yu
 ## is.ggtree <- function(x) inherits(x, 'ggtree')

 
 
 calculate_angle <- function(data) {
     data$angle <- 360/(diff(range(data$y)) + 1) * data$y
     return(data)
 }