## ##' site mask
## ##'
## ##'
## ##' @title mask
## ##' @param tree_object tree object
## ##' @param field selected field
## ##' @param site site
## ##' @param mask_site if TRUE, site will be masked.
## ##'                  if FALSE, selected site will not be masked, while other sites will be masked.
## ##' @return updated tree object
## ##' @export
## ##' @author Guangchuang Yu
## mask <- function(tree_object, field, site, mask_site=FALSE) {
##     has_field <- has.field(tree_object, field)
##     if (has_field == FALSE) {
##         stop("'field' is not available in 'tree_object'...")
##     }

##     has_slot <- attr(has_field, "has_slot")
##     is_codeml <- attr(has_field, "is_codeml")

##     if (has_slot) {
##         if (is_codeml) {
##             field_info <- slot(tree_object@rst, field)
##         } else {
##             field_info <- slot(tree_object, field)
##         }
##         field_data <- field_info[,2]
##     } else {
##         field_data <- tree_object@extraInfo[, field]
##     }

##     field_data <- sapply(field_data, gsub, pattern="\n", replacement="/")

##     x <- field_data[field_data != ""]
##     x <- x[!is.na(x)]
##     pos <- strsplit(x, " / ") %>% unlist %>%
##         gsub("^[a-zA-Z]+", "", . ) %>%
##         gsub("[a-zA-Z]\\s*$", "", .) %>%
##         as.numeric

##     if (mask_site == FALSE) {
##         pos2 <- 1:max(pos)
##         pos2 <- pos2[-site]
##         site <- pos2
##     }

##     site <- site[site %in% pos]

##     for (i in seq_along(field_data)) {
##         if (is.na(field_data[i]))
##             next
##         for (j in seq_along(site)) {
##             pattern <- paste0("/*\\s*[a-zA-Z]", site[j], "[a-zA-Z]\\s*")
##             field_data[i] <- gsub(pattern, "",  field_data[i])
##         }
##         field_data[i] <- gsub("^/\\s", "", field_data[i]) %>% .add_new_line
##     }

##     if (has_slot) {
##         field_info[,2] <- field_data
##         if (is_codeml) {
##             slot(tree_object@rst, field) <- field_info
##         } else {
##             slot(tree_object, field) <- field_info
##         }
##     } else {
##         tree_object@extraInfo[, field] <- field_data
##     }
##     tree_object
## }


## read.tip_seq_mlc <- function(mlcfile) {
##     info <- getPhyInfo(mlcfile)
##     mlc <- readLines(mlcfile)
##     ## remove blank lines
##     blk <- grep("^\\s*$", mlc)
##     if (length(blk) > 0) {
##         mlc <- mlc[-blk]
##     }

##     seqs <- mlc[2:(info$num+1)]
##     seqs <- gsub("\\s+", "", seqs)
##     wd <- info$width
##     res <- sapply(seqs, function(x) substring(x, (base::nchar(x) - wd + 1), base::nchar(x)))
##     nn <- sapply(seqs, function(x) substring(x, 1, (base::nchar(x) - wd)))
##     names(res) <- nn
##     return(res)
## }

## read.tip_seq_mlb <- read.tip_seq_mlc

## read.dnds_mlc <- function(mlcfile) {
##     mlc <- readLines(mlcfile)
##     i <- grep("dN & dS for each branch", mlc)
##     j <- grep("tree length for dN", mlc)

##     mlc <- mlc[i:j]
##     hi <- grep("dN/dS", mlc)
##     cn <- strsplit(mlc[hi], " ") %>% unlist %>% `[`(nzchar(.))

##     ii <- grep("\\d+\\.\\.\\d+", mlc)
##     info <- mlc[ii]
##     info %<>% sub("^\\s+", "", .)
##     info %<>% sub("\\s+$", "", .)
##     res <- t(sapply(info, function(x) {
##         y <- unlist(strsplit(x, "\\s+"))
##         edge <- unlist(strsplit(y[1], "\\.\\."))
##         yy <- c(edge, y[-1])
##         as.numeric(yy)
##     }))

##     row.names(res) <- NULL
##     colnames(res) <- c("parent", "node", cn[-1])
##     colnames(res) <- gsub("\\*", "_x_", colnames(res))
##     colnames(res) <- gsub("\\/", "_vs_", colnames(res))
##     return(res)
## }

## read.treetext_paml_mlc <- function(mlcfile) {
##     read.treetext_paml(mlcfile, "mlc")
## }

## read.treetext_paml_rst<- function(rstfile) {
##     read.treetext_paml(rstfile, "rst")
## }

## read.treetext_paml <- function(file, by) {
##     ## works fine with baseml and codeml
##     x <- readLines(file)
##     tr.idx <- get_tree_index_paml(x)

##     if ( by == "rst" ) {
##         ii <- 1
##     } else if (by == "mlc") {
##         ii <- 3
##     } else {
##         stop("_by_ should be one of 'rst' or 'mlc'")
##     }

##     return(x[tr.idx][ii])
## }

## read.phylo_paml_mlc <- function(mlcfile) {
##     parent <- node <- label <- NULL

##     mlc <- readLines(mlcfile)
##     edge <- get_tree_edge_paml(mlc)

##     tr.idx <- get_tree_index_paml(mlc)
##     tr2 <- read.tree(text=mlc[tr.idx[2]])
##     tr3 <- read.tree(text=mlc[tr.idx[3]])

##     treeinfo <- as.data.frame(edge)
##     colnames(treeinfo) <- c("parent", "node")
##     treeinfo$length <- NA
##     treeinfo$label <- NA
##     treeinfo$isTip <- FALSE
##     ntip <- Ntip(tr3)
##     ii <- match(tr2$tip.label, treeinfo[, "node"])
##     treeinfo[ii, "label"] <- tr3$tip.label
##     treeinfo[ii, "isTip"] <- TRUE
##     ## jj <- match(1:ntip, tr3$edge[,2])
##     treeinfo[ii, "length"] <- tr2$edge.length[ii] ##tr3$edge.length[jj]

##     root <- getRoot(tr3) ## always == (Ntip(tr3) + 1)
##     currentNode <- treeinfo$label[ii]
##     treeinfo.tr3 <- fortify(tr3)

##     treeinfo$visited <- FALSE
##     while(any(treeinfo$visited == FALSE)) {
##         pNode <- c()
##         for( kk in currentNode ) {
##             i <- which(treeinfo$label == kk)
##             treeinfo[i, "visited"] <- TRUE
##             j <- which(treeinfo.tr3$label == kk)
##             ip <- treeinfo[i, "parent"]
##             if (ip != root) {
##                 ii <- which(treeinfo[, "node"] == ip)
##                 if (treeinfo$visited[ii] == FALSE) {
##                     jp <- treeinfo.tr3[j, "parent"]
##                     jj <- which(treeinfo.tr3[, "node"] == jp)
##                     treeinfo[ii, "label"] <- as.character(ip)
##                     treeinfo.tr3[jj, "label"] <- as.character(ip)
##                     treeinfo[ii, "length"] <- treeinfo.tr3[jj, "branch.length"]
##                     pNode <- c(pNode, ip)
##                 }
##                 treeinfo[ii, "visited"] <- TRUE
##             }

##         }
##         currentNode <- unique(pNode)
##     }

##     phylo <- with(treeinfo,
##                   list(edge= cbind(as.numeric(parent),
##                            as.numeric(node)),
##                        edge.length = length,
##                        tip.label = label[order(node)][1:ntip],
##                        Nnode = tr3$Nnode,
##                        node.label = c(root, label[order(node)][-c(1:ntip)])
##                        )
##                   )
##     class(phylo) <- "phylo"
##     phylo <- reorder.phylo(phylo, "cladewise")
##     return(phylo)
## }

## ##' @importFrom ape reorder.phylo
## read.phylo_paml_rst <- function(rstfile) {
##     parent <- node <- label <- NULL

##     ## works fine with baseml and codeml
##     rst <- readLines(rstfile)
##     tr.idx <- get_tree_index_paml(rst)

##     tr1 <- read.tree(text=rst[tr.idx][1])
##     tr3 <- read.tree(text=rst[tr.idx][3])

##     edge <- get_tree_edge_paml(rst)

##     label=c(tr3$tip.label, tr3$node.label)
##     root <- getRoot(tr3)
##     label %<>% `[`(. != root)

##     node.length <- data.frame(label=label,
##                               length=tr1$edge.length)

##     ## node.length$node <- sub("_\\w+", "", node.length$label
##     node.length$node <- gsub("^(\\d+)_.*", "\\1", node.length$label)
##     node.length$label %<>% sub("\\d+_", "", .)

##     edge <- as.data.frame(edge)
##     colnames(edge) <- c("parent", "node")

##     treeinfo <- merge(edge, node.length, by.x="node", by.y="node")
##     edge2 <- treeinfo[, c("parent", "node")]
##     edge2 %<>% as.matrix

##     ntip <- Ntip(tr3)

##     phylo <- with(treeinfo,
##                   list(edge= cbind(as.numeric(parent),
##                       as.numeric(node)),
##                        edge.length = length,
##                        tip.label = label[order(node)][1:ntip],
##                        Nnode = tr3$Nnode,
##                        node.label = c(root, label[order(node)][-c(1:ntip)])
##                        )
##                   )

##     class(phylo) <- "phylo"
##     phylo <- reorder.phylo(phylo, "cladewise")

##     return(phylo)
## }

## read.ancseq_paml_rst <- function(rstfile, by="Marginal") {
##     ## works fine with baseml and codeml
##     rst <- readLines(rstfile)

##     by <- match.arg(by, c("Marginal", "Joint"))
##     query <- paste(by, "reconstruction of ancestral sequences")
##     idx <- grep(query, rst)
##     if(length(idx) == 0) {
##         ## in some paml setting, joint_ancseq are not available.
##         return("")
##     }
##     si <- grep("reconstructed sequences", rst)
##     idx <- si[which.min(abs(si-idx))]

##     nl <- strsplit(rst[idx+2], split=" ") %>% unlist %<>% `[`(nzchar(.))
##     N <- as.numeric(nl[1])
##     seq.leng <- as.numeric(nl[2])

##     seqs <- rst[(idx+4):(idx+3+N)]

##     seq.name <- character(N)
##     res <- character(N)
##     for (i in 1:N) {
##         ss <- gsub(" ", "", seqs[i])
##         nn <- base::nchar(ss)
##         res[i] <- substring(ss, nn-seq.leng+1,nn)
##         seq.name[i] <- substring(ss, 1, nn-seq.leng)
##     }
##     seq.name <- sub("\\w+#", "", seq.name)
##     names(res) <- seq.name

##     return(res)
## }


## get_tree_index_paml <- function(paml) {
##     grep("\\)[ \\.0-9]*;", paml)
## }

## get_tree_edge_index_paml <- function(paml) {
##     grep("\\d+\\.\\.\\d+", paml)
## }

## get_tree_edge_paml <- function(paml) {
##     tr.idx <- get_tree_index_paml(paml)

##     edge.idx <- get_tree_edge_index_paml(paml)
##     edge.idx <- edge.idx[edge.idx < tr.idx[3]]

##     nodeNum <- strsplit(paml[edge.idx], split="\\.\\.") %>%
##         unlist %>% strsplit(split="[[:space:]]") %>% unlist

##     nodeNum %<>% `[`(nzchar(.))

##     edge <- matrix(as.numeric(nodeNum), ncol=2, byrow = TRUE)

##     return(edge)
## }

## set.paml_rst_ <- function(object) {
##     if (!is(object, "paml_rst")) {
##         stop("object should be an instance of 'paml_rst'")
##     }
##     if (length(object@tip_seq) == 0) {
##         return(object)
##     }

##     types <- get.fields(object)
##     for (type in types) {
##         value <- subs_paml_rst(object, type)
##         if (all(is.na(value)))
##             next

##         if (type == "marginal_subs") {
##             object@marginal_subs <- value
##         } else if (type == "marginal_AA_subs") {
##             object@marginal_AA_subs <- value
##         } else if (type == "joint_subs") {
##             object@joint_subs <- value
##         } else if (type == "joint_AA_subs") {
##             object@joint_AA_subs <- value
##         }
##     }
##     return(object)
## }


## get.subs_paml_rst <- function(object, type) {
##     if (!is(object, "paml_rst")) {
##         stop("object should be an instance of 'paml_rst'")
##     }
##     if (type == "marginal_subs") {
##         res <- object@marginal_subs
##     } else if (type == "marginal_AA_subs") {
##         res <- object@marginal_AA_subs
##     } else if (type == "joint_subs") {
##         res <- object@joint_subs
##     } else if (type == "joint_AA_subs") {
##         res <- object@joint_AA_subs
##     } else {
##         stop("type should be one of 'marginal_subs',
##                              'marginal_AA_subs', 'joint_subs' or 'joint_AA_subs'. ")
##     }
##     return(res)
## }

## subs_paml_rst <- function(x, type, ...) {
##     if (class(x) != "paml_rst") {
##         stop("x should be an object of paml_rst...")
##     }
##     seqs <- x@tip_seq
##     if (length(seqs) == 0) {
##         stop("tip sequences is not available...")
##     }
##     if (type %in% c("marginal_subs", "marginal_AA_subs")) {
##         ancseq <- x@marginal_ancseq
##         ## seqs <- c(seqs, x@marginal_ancseq)
##     } else if (type %in% c("joint_subs", "joint_AA_subs")){
##         ancseq <- x@joint_ancseq
##         ## seqs <- c(seqs, x@joint_ancseq)
##     } else {
##         stop("type should be one of 'marginal_subs',
##                              'marginal_AA_subs', 'joint_subs' or 'joint_AA_subs'. ")
##     }
##     if( type %in% c("marginal_subs", "joint_subs")) {
##         translate <- FALSE
##     } else {
##         translate <- TRUE
##     }

##     if (all(ancseq == "")) {
##         return(NA)
##     }
##     seqs <- c(seqs, ancseq)
##     get.subs_(x@phylo, seqs, translate=translate, ...)
## }