##' convert polytomy to binary tree
 ##'
 ##' as.binary method for \code{phylo} object
 ##' @rdname as.binary
 ##' @return binary tree
 ##' @method as.binary phylo
 ##' @importFrom ape is.binary.tree
 ##' @export
 ##' @author Guangchuang Yu \url{http://ygc.name}
 ##' @examples
 ##' require(ape)
 ##' tr <- read.tree(text="((A, B, C), D);")
 ##' is.binary.tree(tr)
 ##' tr2 <- as.binary(tr)
 ##' is.binary.tree(tr2)
 as.binary.phylo <- function(tree, ...) {
     if(is.binary.tree(tree)) {
         cat("The input tree is already binary...\n")
         invisible(tree)
     }
     
     polyNode <- tree$edge[,1] %>% table %>% '>'(2) %>%
         which %>% names %>% as.numeric
 
     N <- getNodeNum(tree)
     ii <- 0
     for (pn in polyNode) {
         idx <- which(tree$edge[,1] == pn)
         while(length(idx) >2) {
             ii <- ii + 1
             newNode <- N+ii
             tree$edge[idx[-1],1] <- newNode
             newEdge <- matrix(c(tree$edge[idx[1],1], newNode), ncol=2)
             tree$edge <- rbind(tree$edge, newEdge)
             idx <- idx[-1]
         }
     }
         
     tree$Nnode <- tree$Nnode+ii
     tree$edge.length <- c(tree$edge.length, rep(0, ii))
     return(tree)
 }
 
 
 ##' remove singleton
 ##'
 ##' 
 ##' @title rm.singleton.newick
 ##' @param nwk newick file
 ##' @param outfile output newick file 
 ##' @return tree text
 ##' @importFrom magrittr %<>%
 ##' @importFrom magrittr add
 ##' @importFrom ape write.tree
 ##' @author Guangchuang Yu \url{http://ygc.name}
 rm.singleton.newick <- function(nwk, outfile = NULL) {    
     tree <- readLines(nwk)
 
     ## remove singleton of tips
     nodePattern <- "\\w+:[\\.0-9Ee\\+\\-]+"
     singletonPattern.with.nodename <- paste0(".*(\\(", nodePattern, "\\)\\w+:[\\.0-9Ee\\+\\-]+).*")

     singletonPattern.wo.nodename <- paste0(".*(\\(", nodePattern, "\\):[\\.0-9Ee\\+\\-]+).*")

     
     while(length(grep("\\([^,]+\\)", tree)) > 0) {
         singleton <- gsub(singletonPattern.with.nodename, "\\1", tree)
         if (singleton == tree) {
             singleton <- gsub(singletonPattern.wo.nodename, "\\1", tree)
         }
         if (singleton == tree) {
             stop("can't parse singleton node...")
         }
 
         tip <- gsub("\\((\\w+).*", "\\1", singleton)
         
         len1 <- gsub(".*[^\\.0-9Ee\\+\\-]+([\\.0-9Ee\\+\\-]+)", "\\1", singleton)
         len2 <- gsub(".*:([\\.0-9Ee\\+\\-]+)\\).*", "\\1", singleton)
         len <- as.numeric(len1) + as.numeric(len2)
         
         tree <- gsub(singleton, paste0(tip, ":", len), tree, fixed = TRUE)
     }
 
     tree <- read.tree(text=tree)
 
     ### remove singleton of internal nodes
     p.singleton <- which(table(tree$edge[,1]) == 1)
     if (length(p.singleton) > 0) {
         p.singleton %<>% names %>% as.numeric
         edge <- tree$edge
         idx <- which(edge[,1] == p.singleton)
         singleton <- edge[idx, 2]
         sidx <- which(edge[,1] == singleton)
         edge[sidx,1] <- p.singleton
         edge <- edge[-idx,]
         tree$edge <- edge
         tree$edge.length[sidx] %<>% add(., tree$edge.length[idx])
         tree$edge.length <- tree$edge.length[-idx]
     }
     
     if (!is.null(outfile)) {
         write.tree(tree, file=outfile)
     }
     invisible(tree)
 }
 
 ##' @method fortify beast
 ##' @export
 fortify.beast <- function(model, data,

                           layout        = "rectangular",
                           yscale        = "none",
                           ladderize     = TRUE,
                           right         = FALSE,
                           branch.length = "branch.length",
                           ndigits       = NULL,

                           mrsd = NULL, ...) {
 

     phylo <- set_branch_length(model, branch.length)
 

     df    <- fortify(phylo, layout=layout, branch.length=branch.length,

                      ladderize=ladderize, right=right, mrsd = mrsd, ...)
     
     stats <- model@stats
 
     scn <- colnames(stats)
     scn <- scn[scn != 'node']
     
     for (cn in scn) {
         if (cn %in% colnames(df)) {
             colnames(stats)[colnames(stats) == cn] <- paste0(cn, "_")
             msg <- paste("feature", cn, "was renamed to", paste0(cn, "_"), "due to name conflict...")
             warning(msg)
         }
     }
 
     idx <- which(colnames(stats) != "node")
     for (ii in idx) {
         if (is.character_beast(stats, ii)) {
             len <- sapply(stats[,ii], length)
             if (any(len > 1)) {
                 stats[,ii] %<>% sapply(., function(x) {
                     y <- unlist(x) %>% as.character %>%
                         gsub("\"", "", .) %>% gsub("'", "", .)
                     if (length(y) == 1) {
                         return(y)
                     } else {
                         return(paste0('{', paste0(y, collapse = ','), '}'))
                     }
                 })
             } else {
                 stats[,ii] %<>% unlist %>% as.character %>%
                     gsub("\"", "", .) %>% gsub("'", "", .)
             }
             next
         }
         
         len <- sapply(stats[,ii], length)
         if ( all(len == 1) ) {
             stats[, ii] %<>% unlist %>% as.character %>% as.numeric
             if (!is.null(ndigits)) {
                 stats[, ii] %<>% round(., ndigits)
             }
         } else if (all(len <= 2)) {
             stats[, ii] %<>% sapply(., function(x) {
                 y <- unlist(x) %>% as.character %>% as.numeric
                 if (!is.null(ndigits)) {
                     y %<>% round(., ndigits)
                 }
                 if (length(y) == 1) {
                     return(y)
                 } else {
                     return(paste0('[', paste0(y, collapse = ','), ']'))
                 }
             })
         } else {
             stats[,ii] %<>% sapply(., function(x) {
                 y <- unlist(x) %>% as.character %>% as.numeric
                 if (!is.null(ndigits)) {
                     y %<>% round(., ndigits)
                 }
                 if (length(y) == 1) {
                     return(y)
                 } else {
                     return(paste0('{', paste0(y, collapse = ','), '}'))
                 }
             })  
         }
     }
             
       
     cn <- colnames(stats)
     lo <- cn[grep("_lower", cn)]
     hi <- gsub("lower$", "upper", lo)
     rid <- gsub("_lower$", "", lo)
     
     for (i in seq_along(rid)) {
         stats[, rid[i]] <- paste0("[", stats[, lo[i]], ",", stats[, hi[i]], "]")
         stats[is.na(stats[, lo[i]]), rid[i]] <- NA
     }
     
     idx   <- match(df$node, stats$node)
     stats <- stats[idx,]
     cn_stats <- colnames(stats)
     stats <- stats[, cn_stats != "node"]
     
     df <- cbind(df, stats)
     if (is(stats, "data.frame") == FALSE) {
         colnames(df)[colnames(df) == "stats"] <- cn_stats[cn_stats != "node"]
     }
     
     df <- scaleY(phylo, df, yscale, layout, ...)
 
     append_extraInfo(df, model)
 }
 
 
 ##' @method fortify codeml
 ##' @export
 fortify.codeml <- function(model, data,
                            layout        = "rectangular",
                            yscale        = "none",
                            ladderize     = TRUE,
                            right         = FALSE,
                            branch.length = "mlc.branch.length",
                            ndigits       = NULL,
                            mrsd          = NULL,
                            ...) {
 
     dNdS <- model@mlc@dNdS
     if (branch.length == "branch.length") {
         message("branch.length setting to mlc.branch.length by default...")
         branch.length <- "mlc.branch.length"
     }
     length <- match.arg(branch.length,
                         c("none",
                           "mlc.branch.length",
                           "rst.branch.length",
                           colnames(dNdS)[-c(1,2)])
                         )
     
     if (length == "rst.branch.length") {
         phylo <- get.tree(model@rst)
     } else {
         if (length == "mlc.branch.length") {

             length <- "branch.length"

         }

         phylo <- set_branch_length(model@mlc, length)

     }
     
     df <- fortify(phylo, data, layout, ladderize, right,
                   branch.length=length, mrsd=mrsd, ...)
     
     res <- merge_phylo_anno.codeml_mlc(df, dNdS, ndigits)
     df <- merge_phylo_anno.paml_rst(res, model@rst)
     df <- scaleY(phylo, df, yscale, layout, ...)

     append_extraInfo(df, model)
 }
 
 
 ##' @method fortify codeml_mlc
 ##' @export
 fortify.codeml_mlc <- function(model, data,
                                layout        = "rectangular",
                                yscale        = "none",
                                ladderize     = TRUE,
                                right         = FALSE,
                                branch.length = "branch.length",
                                ndigits       = NULL,
                                mrsd          = NULL,
                                ...) {

 
     phylo <- set_branch_length(model, branch.length)

     df <- fortify(phylo, data, layout, ladderize, right,
                   branch.length=branch.length, mrsd=mrsd, ...)

     
     dNdS <- model@dNdS
 
     df <- merge_phylo_anno.codeml_mlc(df, dNdS, ndigits)
     df <- scaleY(phylo, df, yscale, layout, ...)
 
     append_extraInfo(df, model)
 }
 
 merge_phylo_anno.codeml_mlc <- function(df, dNdS, ndigits = NULL) {
     if (!is.null(ndigits)) {
         idx <- which(! colnames(dNdS) %in% c("node", "parent"))
         for (ii in idx) {
             if (is.numeric(dNdS[, ii])) {
                 dNdS[, ii] <- round(dNdS[,ii], ndigits)
             }
         }
     }
     
     res <- merge(df, dNdS,
                  by.x  = c("node", "parent"),
                  by.y  = c("node", "parent"),
                  all.x = TRUE)
     
     res[match(df$node, res$node),]
 }
 
 fortify.codeml_mlc_ <- function(model, data,
                                 layout        = "rectangular",
                                 ladderize     = TRUE,
                                 right         = FALSE,
                                 branch.length = "branch.length",
                                 ...) {
 
 }
 

     
 ##' @method fortify paml_rst
 ##' @export
 fortify.paml_rst <- function(model, data, layout = "rectangular", yscale="none",
                              ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {
     df <- fortify.phylo(model@phylo, data, layout, ladderize, right, mrsd=mrsd, ...)
     df <- merge_phylo_anno.paml_rst(df, model)
     df <- scaleY(model@phylo, df, yscale, layout, ...)
 
     append_extraInfo(df, model)
 }
 
 merge_phylo_anno.paml_rst <- function(df, model) {
     for (type in get.fields(model)) {
         anno <- get.subs(model, type=type)
         colnames(anno)[2] <- type
         df <- df %add2% anno
     }
     return(df)
 }
 
 
 ##' @method fortify phangorn
 ##' @export
 fortify.phangorn <- fortify.paml_rst
 
 
 ##' @method fortify hyphy
 ##' @export
 fortify.hyphy <- fortify.paml_rst
 
 
 ##' @method fortify jplace
 ##' @importFrom ape read.tree
 ##' @export
 fortify.jplace <- function(model, data,
                            layout="rectangular", yscale="none",
                            ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {
     df <- get.treeinfo(model, layout, ladderize, right, mrsd=mrsd, ...)
     place <- get.placements(model, by="best")
 
     df <- df %add2% place
 
     df <- scaleY(model@phylo, df, yscale, layout, ...)
 
     append_extraInfo(df, model)    
 }
 
 scaleY <- function(phylo, df, yscale, layout, ...) {
     if (yscale == "none") {
         return(df)
     } 
     if (! yscale %in% colnames(df)) {
         warning("yscale is not available...\n")
         return(df)
     }
     if (is.numeric(df[, yscale])) {
         y <- getYcoord_scale_numeric(phylo, df, yscale, ...)
         ## if (order.y) {
         ##     y <- getYcoord_scale2(phylo, df, yscale)
         ## } else {
         ##     y <- getYcoord_scale(phylo, df, yscale)
         ## }
     } else {
         y <- getYcoord_scale_category(phylo, df, yscale, ...)
     }
     
     df[, "y"] <- y
 
     return(df)
 }
 
 
 ##' @method fortify phylo4
 ##' @export
 fortify.phylo4 <- function(model, data, layout="rectangular", yscale="none",
                            ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {
     phylo <- as.phylo.phylo4(model)
     df <- fortify.phylo(phylo, data,
                         layout, ladderize, right, mrsd=mrsd, ...)
     scaleY(phylo, df, yscale, layout, ...)
 }
 

 ##' @method fortify phylo4d
 ##' @export
 fortify.phylo4d <- function(model, data, layout="rectangular", yscale="none",
                             ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {
     res <- fortify.phylo4(model, data, layout, yscale, ladderize, right, mrsd, ...)
     tdata <- model@data[match(res$node, rownames(model@data)), , drop=FALSE]
     cbind(res, tdata)
 }
 

 as.phylo.phylo4 <- function(phylo4) {
     edge <- phylo4@edge
     edge <- edge[edge[,1] != 0, ]
     edge.length <- phylo4@edge.length
     edge.length <- edge.length[!is.na(edge.length)]
     tip.id <- sort(setdiff(edge[,2], edge[,1]))
     tip.label <- phylo4@label[tip.id]
     phylo <- list(edge = edge,
                   edge.length = edge.length,
                   tip.label = tip.label)
     
     node.id <- sort(unique(edge[,1]))
     node.id <- node.id[node.id != 0]
     node.label <- phylo4@label[node.id]
     if (!all(is.na(node.label))) {
         phylo$node.label <- node.label
     }
     phylo$Nnode <- length(node.id)
     class(phylo) <- "phylo"
     return(phylo)
 }
 
 ##' fortify a phylo to data.frame
 ##'
 ##' 
 ##' @rdname fortify
 ##' @title fortify
 ##' @param model phylo object
 ##' @param data not use here
 ##' @param layout layout
 ##' @param ladderize ladderize, logical
 ##' @param right logical
 ##' @param mrsd most recent sampling date
 ##' @param as.Date logical whether using Date class in time tree
 ##' @param ... additional parameter
 ##' @return data.frame
 ##' @importFrom ape ladderize
 ##' @importFrom ggplot2 fortify
 ##' @method fortify phylo
 ##' @export
 ##' @author Yu Guangchuang
 fortify.phylo <- function(model, data, layout="rectangular", 
                           ladderize=TRUE, right=FALSE, mrsd=NULL, as.Date=FALSE, ...) {
     if (ladderize == TRUE) {
         tree <- ladderize(model, right=right)
     } else {
         tree <- model
     }

 
     if (! is.null(tree$edge.length)) {

         if (anyNA(tree$edge.length)) {

             warning("'edge.length' contains NA values...\n## setting 'edge.length' to NULL automatically when plotting the tree...")
             tree$edge.length <- NULL
         }
     }

     
     df <- as.data.frame(tree, layout=layout, ...)
     idx <- is.na(df$parent)
     df$parent[idx] <- df$node[idx]
     rownames(df) <- df$node
     cn <- colnames(df)
     colnames(df)[grep("length", cn)] <- "branch.length"
     if(layout == "slanted") {
         df <- add_angle_slanted(df)
     }
     aa <- names(attributes(tree))
     group <- aa[ ! aa %in% c("names", "class", "order", "reroot", "node_map")]
     if (length(group) > 0) {
         for (group_ in group) {
             ## groupOTU & groupClade
             group_info <- attr(tree, group_)
             if (length(group_info) == nrow(df)) {
                 df[, group_] <- group_info
             }
         }
     }
     
     if (!is.null(mrsd)) {

         df <- scaleX_by_time_from_mrsd(df, mrsd, as.Date)

     }
     return(df)
 }
 
 ##' convert phylo to data.frame
 ##'
 ##' 
 ##' @title as.data.frame
 ##' @param x phylo object
 ##' @param row.names omitted here
 ##' @param optional omitted here
 ##' @param layout layout
 ##' @param ... additional parameter
 ##' @return data.frame
 ##' @method as.data.frame phylo
 ##' @export
 ##' @author Yu Guangchuang
 as.data.frame.phylo <- function(x, row.names, optional,
                                 layout="rectangular", ...) {
     if (layout == "unrooted") {
         return(layout.unrooted(x))
     } 
     as.data.frame.phylo_(x, layout, ...)
 }
 
 as.data.frame.phylo_ <- function(x, layout="rectangular",
                                  branch.length="branch.length", ...) {

     if (branch.length != 'none') {
         branch.length = "branch.length"
     }
     

     tip.label <- x[["tip.label"]]
     Ntip <- length(tip.label)
     N <- getNodeNum(x)
     
     edge <- as.data.frame(x[["edge"]])
     colnames(edge) <- c("parent", "node")

     if (! is.null(x$edge.length)) {
         edge$length <- x$edge.length
         if (branch.length == "none") {
             xpos <- getXcoord_no_length(x)
             ypos <- getYcoord(x)
         } else {
             xpos <- getXcoord(x)
             ypos <- getYcoord(x)
         }
         ## } else  if (layout != "cladogram") {
         ##     xpos <- getXcoord(x)
         ##     ypos <- getYcoord(x)
         ## } else {
         ##     ## layout == "cladogram" && branch.length != "none"
         ##     xy <- getXYcoord_cladogram(x)
         ##     xpos <- xy$x
         ##     ypos <- xy$y
         ## }
     } else {
         xpos <- getXcoord_no_length(x)
         ypos <- getYcoord(x)
     }
     
     xypos <- data.frame(node=1:N, x=xpos, y=ypos)
 
     res <- merge(edge, xypos, by.x="node", by.y="node", all.y=TRUE)
     label <- rep(NA, N)
     label[1:Ntip] <- tip.label
     if ( !is.null(x$node.label) ) {
         label[(Ntip+1):N] <- x$node.label
     }
     res$label <- label
     isTip <- rep(FALSE, N)
     isTip[1:Ntip] <- TRUE
     res$isTip <- isTip
 
     ## add branch mid position
     res <- calculate_branch_mid(res)
 

     ## ## angle for all layout, if 'rectangular', user use coord_polar, can still use angle

     res <- calculate_angle(res)

     return(res)
 }
 
 ##' @method fortify nhx
 ##' @export
 fortify.nhx <- function(model, data, layout= "rectangular",
                         ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {
     df <- fortify(get.tree(model), layout=layout, ladderize=ladderize, right=right, mrsd=mrsd, ...)
     df <- merge(df, model@nhx_tags, by.x="node", by.y="node", all.x=TRUE)
     append_extraInfo(df, model)
 }
 
 
 ##' @method fortify raxml
 ##' @export
 fortify.raxml <- function(model, data, layout= "rectangular",
                           ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {
     df <- fortify(get.tree(model), layout=layout, ladderize=ladderize, right=right, mrsd=mrsd, ...)
     df <- merge(df, model@bootstrap, by.x="node", by.y="node", all.x=TRUE)
     append_extraInfo(df, model)
 }
 
 ##' @method fortify apeBootstrap
 ##' @export
 fortify.apeBootstrap <- fortify.raxml
 
 
 ##' @method fortify multiPhylo
 ##' @export
 fortify.multiPhylo <-  function(model, data, layout="rectangular", 
                                 ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {
 
     df.list <- lapply(model, function(x) fortify(x, layout=layout, ladderize=ladderize, right=right, mrsd=mrsd, ...))
     if (is.null(names(model))) {
         names(df.list) <- paste0("Tree ", "#", seq_along(model))
     } else {
         names(df.list) <- names(model)
     }
     df <- do.call("rbind", df.list)
     df$.id <- rep(names(df.list), times=sapply(df.list, nrow))
     df$.id <- factor(df$.id, levels=names(df.list))
     

     ## nNode <- sapply(df.list, nrow)
     ## nNode2 <- cumsum(c(0, nNode[-length(nNode)])) 
     ## df$parent <- df$parent + rep(nNode2, times=nNode)

     return(df)
 }
 

 ##' @method fortify phylip
 ##' @export
 fortify.phylip <- function(model, data, layout="rectangular",
                         ladderize=TRUE, right=FALSE,
                         branch.length = "TREE", mrsd=NULL, ...) {
     trees <- get.tree(model)
     fortify(trees, layout=layout, ladderize = ladderize, right=right, mrsd=mrsd, ...)
 }
     

 ##' @method fortify r8s
 ##' @export
 fortify.r8s <- function(model, data, layout="rectangular",
                         ladderize=TRUE, right=FALSE,
                         branch.length = "TREE", mrsd=NULL, ...) {
     trees <- get.tree(model)
     branch.length %<>% match.arg(names(trees))
     phylo <- trees[[branch.length]]
     fortify(phylo, layout=layout, ladderize = ladderize, right=right, mrsd=mrsd, ...)
 }

 
 ##' @method fortify obkData
 ##' @export
 fortify.obkData <- function(model, data, layout="rectangular",
                             ladderize=TRUE, right=FALSE, mrsd = NULL, ...) {
 
     df <- fortify(model@trees[[1]], layout=layout, ladderize=ladderize, right=right, mrsd=mrsd, ...)
 
     meta.df <- model@dna@meta
     meta.df <- data.frame(taxa=rownames(meta.df), meta.df)
     loc <- model@individuals
     loc <- data.frame(individualID=rownames(loc), loc)
     meta_loc <- merge(meta.df, loc, by="individualID")
     meta_loc <- meta_loc[,-1]
 
     df <- merge(df, meta_loc, by.x="label", by.y="taxa", all.x=TRUE)
     df <- df[order(df$node, decreasing = FALSE),]
     return(df)
 }

 
 ##' @method fortify phyloseq
 ##' @export
 fortify.phyloseq <- function(model, data, layout="rectangular",
                              ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {
 
     df <- fortify(model@phy_tree, layout=layout, ladderize=ladderize, right=right, mrsd=mrsd, ...)
     phyloseq <- "phyloseq"
     require(phyloseq, character.only=TRUE)
     psmelt <- eval(parse(text="psmelt"))
     dd <- psmelt(model)
     if ('Abundance' %in% colnames(dd)) {
         dd <- dd[dd$Abundance > 0, ]
     }
     
     data <- merge(df, dd, by.x="label", by.y="OTU", all.x=TRUE)
     spacing <- 0.02
     idx <- with(data, sapply(table(node)[unique(node)], function(i) 1:i)) %>% unlist
     data$hjust <- spacing * idx * max(data$x)
     ## data$hjust <- data$x + hjust
 
     data[order(data$node, decreasing = FALSE), ]
 }
 

                          
 ## fortify.cophylo <- function(model, data, layout="rectangular",
 ##                             ladderize=TRUE, right=FALSE, mrsd = NULL, ...) {
 ##     trees <- model$trees
 ##     df.list <- lapply(trees, function(x) fortify(x, layout=layout, ladderize=ladderize, right=right, mrsd=mrsd, ...))
 ##     df1 <- df.list[[1]]
 ##     df2 <- df.list[[2]]
 
 ##     df2$x <- max(df2$x) + df2$x * -1 + max(df1$x) * 1.1
 ##     df2$parent <- df2$parent+nrow(df1)
 ##     df <- rbind(df1, df2)
 
 ##     ggplot(df) + geom_tree()
 
 ## }

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