R/barPhylo.R

In fkeck/phylosignal: Exploring the Phylogenetic Signal in Continuous Traits

#' Plots of Traits Values along a Phylogeny
#'
#' This function provides a general interface to plot \code{phylo4d} object
#' (i.e. phylogenetic tree and data).
#' 
#' @param p4d a \code{phylo4d} object.
#' @param trait the traits in the \code{phylo4d} object to include in the plot.
#' Can be a character vector giving the name of the traits or numbers giving the column index
#' in the table of the data slot of the p4d object. Can be used to reorder the traits in the plot.
#' @param center a logical indicating whether traits values should be centered.
#' @param scale a logical indicating whether traits values should be scaled.
#' @param plot.type a character string specifying the type of plot for traits data.
#' Can be "\code{barplot}", "\code{dotplot}" or "\code{gridplot}".
#' 
#' @param tree.ladderize a logical indicating whether the tree should be (right) ladderized.
#' @param tree.type a character string specifying the type of phylogeny to be drawn.
#' Can be "\code{phylogram}", "\code{cladogram}" or "\code{fan}".
#' @param tree.ratio a numeric value in [0, 1] giving the proportion of width of the figure for the tree.
#' @param tree.xlim a numeric vector of length 2 giving the limits of the x-axis for the tree. If \code{NULL},
#' it is determined automatically.
#' @param tree.open.angle a numeric value giving the angle in degrees left blank if \code{tree.type = "fan"}.
#' @param tree.open.crown a logical indicating whether the crowns should be drawn following the value
#' of \code{tree.open.angle} (default \code{TRUE}).
#' 
#' @param show.tip logical indicating whether tips labels should be drawn.
#' @param tip.labels character vector to label the tips.
#' If \code{NULL} the tips labels of the \code{phylo4d} object are used
#' @param tip.col a vector of R colors to use for the tips labels.
#' Recycled if necessary.
#' @param tip.cex a numeric vector to control character size of the tips labels.
#' Recycled if necessary.
#' @param tip.font an integer vector specifying the type of font for the tips labels:
#' 1 (plain text), 2 (bold), 3 (italic), or 4 (bold italic). Recycled if necessary.
#' @param tip.adj a vector of numeric in [0, 1] to control tips labels justification:
#' 0 (left-justification), 0.5 (centering), or 1 (right-justification). Recycled if necessary.
#' 
#' @param data.xlim numeric vector of length 2 or matrix giving the x coordinates range for
#' the barplots/dotplots (see Details).
#' @param bar.lwd a vector of numeric giving bar widths of the barplot(s).
#' Recycled along the tips, reapeated for each trait.
#' @param bar.col a vector of R colors to use for the bars.
#' Recycled along the tips, reapeated for each trait.
#' The user can also provide a matrix for a finer tuning (see Details)
#' @param show.data.axis logical indicating whether barplots/dotplots axes should be drawn.
#' 
#' @param dot.col a vector of R colors to use for the points.
#' Recycled along the tips, reapeated for each trait.
#' The user can also provide a matrix for a finer tuning (see Details)
#' @param dot.pch a numerical vector of symbol to use for the points.
#' Recycled along the tips, reapeated for each trait.
#' The user can also provide a matrix for a finer tuning (see Details)
#' @param dot.cex a numerical vector. Character (or symbol) expansion for the points.
#' Recycled along the tips, reapeated for each trait.
#' The user can also provide a matrix for a finer tuning (see Details)
#' 
#' @param cell.col a vector of colors for \code{gridplot} cells.
#' Easily generated by \code{\link[grDevices]{heat.colors}},
#' \code{\link[grDevices]{topo.colors}}, \code{\link[grDevices]{terrain.colors}}
#' or other functions created with \code{\link[grDevices]{colorRampPalette}}.
#' @param show.color.scale logical indicating whether color scale should be drawn.
#' 
#' @param show.trait logical indicating whether traits labels should be drawn.
#' @param trait.labels character vector to label the traits.
#' If \code{NULL} the traits labels of the \code{phylo4d} object are used.
#' @param trait.col a vector of R colors to use for the traits labels.
#' Recycled if necessary.
#' @param trait.cex a numeric vector to control character size of the trait labels.
#' Recycled if necessary.
#' @param trait.font an integer vector specifying the type of font for the traits labels:
#' 1 (plain text), 2 (bold), 3 (italic), or 4 (bold italic). Recycled if necessary.
#' @param trait.bg.col a vector of R colors to use for the background of the barplots.
#' Recycled if necessary.
#' 
#' @param error.bar.sup a matrix giving the superior limit for error bars.
#' Columns and rows names must match with traits and tips labels, respectively.
#' @param error.bar.inf a matrix giving the inferior limit for error bars.
#' Columns and rows names must match with traits and tips labels, respectively.
#' @param error.bar.col a vector of R colors to use for the bars.
#' Recycled along the tips, reapeated for each trait.
#' The user can also provide a matrix for a finer tuning (see Details)
#' 
#' @param show.box a logical indicating whether a box should be drawn around the plots.
#' @param grid.vertical a logical incating whether vertical lines of the grid should be drawn.
#' @param grid.horizontal a logical incating whether horizontal lines of the grid should be drawn.
#' @param grid.col a vector of R colors to use for the lines of the grid.
#' @param grid.lty the lines type of the grid. Possibly a vector.
#' @param ... further arguments to be passed to \code{plot.phylo}.
#'  
#' @export
multiplot.phylo4d <- function(p4d, trait = names(tdata(p4d)), center = TRUE, scale = TRUE, plot.type = "barplot",
                              tree.ladderize = FALSE, tree.type = "phylogram", tree.ratio = NULL, tree.xlim = NULL,
                              tree.open.angle = 0, tree.open.crown = TRUE,
                              show.tip = TRUE, tip.labels = NULL, tip.col = "black", tip.cex = 1, tip.font = 3, tip.adj = 0,
                              data.xlim = NULL, bar.lwd = 10, bar.col = "grey35", show.data.axis = TRUE,
                              dot.col = "black", dot.pch = 20, dot.cex = 2,
                              cell.col = white2red(100), show.color.scale = TRUE,
                              show.trait = TRUE, trait.labels = NULL, trait.col = "black", trait.cex = 1, trait.font = 1,
                              trait.bg.col = "grey90", error.bar.sup = NULL, error.bar.inf = NULL, error.bar.col = 1,
                              show.box = FALSE, grid.vertical = TRUE, grid.horizontal = FALSE, grid.col = "grey25",
                              grid.lty = "dashed", ...){
  
  
  p4 <- phylobase::extractTree(p4d)
  phy <- as(p4, "phylo")
  if(tree.ladderize){
    phy <- ladderize(phy)
  }
  new.order <- phy$edge[, 2][!phy$edge[, 2] %in% phy$edge[, 1]]
  tips <- phy$tip.label[new.order]
  n.tips <- length(tips)
  X <- tdata(p4d, type = "tip")
  X <- X[new.order, trait]
  X <- scale(X, center = center, scale = scale)
  X <- as.data.frame(X)
  colnames(X) <- trait
  n.traits <- ncol(X)
  if(is.numeric(trait)){
    trait <- names(tdata(p4d))[trait]
  }
  ####
  tree.type <- match.arg(tree.type, c("phylogram", "cladogram", "fan"))
  plot.type <- match.arg(plot.type, c("barplot", "dotplot", "gridplot"))
  
  if(!is.null(error.bar.inf)){
    error.bar.inf <- as.matrix(error.bar.inf)
    error.bar.inf[is.na(error.bar.inf)] <- 0
    error.bar.inf <- matchTipsAndTraits(error.bar.inf, p4d.tips = tips, p4d.traits = trait)
  }
  if(!is.null(error.bar.sup)){
    error.bar.sup <- as.matrix(error.bar.sup)
    error.bar.sup[is.na(error.bar.sup)] <- 0
    error.bar.sup <- matchTipsAndTraits(error.bar.sup, p4d.tips = tips, p4d.traits = trait)
  }
  
  if(!is.null(error.bar.inf)){
    arrow.inf <- X
    arrow.inf[X>0] <- X[X>0] - error.bar.inf[X>0]
    arrow.inf[X<0] <- X[X<0] + error.bar.inf[X<0]
  }
  if(!is.null(error.bar.sup)){
    arrow.sup <- X
    arrow.sup[X>0] <- X[X>0] + error.bar.sup[X>0]
    arrow.sup[X<0] <- X[X<0] - error.bar.sup[X<0]
  }
  
  if(is.null(data.xlim)){
    if(center & scale){
      data.xlim <- matrix(rep(NA, n.traits * 2), nrow = 2,
                          dimnames = list(c("xlim.min", "xlim.max"), trait))
      if(!is.null(error.bar.inf) & !is.null(error.bar.sup)){
        data.xlim[1, ] <- floor(min(arrow.inf, arrow.sup, na.rm = TRUE))
        data.xlim[2, ] <- ceiling(max(arrow.inf, arrow.sup, na.rm = TRUE))
      } else if(!is.null(error.bar.inf)){
        data.xlim[1, ] <- floor(min(arrow.inf, na.rm = TRUE))
        data.xlim[2, ] <- ceiling(max(arrow.inf, na.rm = TRUE))
      } else if(!is.null(error.bar.sup)){
        data.xlim[1, ] <- floor(min(arrow.sup, na.rm = TRUE))
        data.xlim[2, ] <- ceiling(max(arrow.sup, na.rm = TRUE))
      } else {
        data.xlim[1, ] <- floor(min(X, na.rm = TRUE))
        data.xlim[2, ] <- ceiling(max(X, na.rm = TRUE))
      }
    } else {
      data.xlim <- matrix(NA, nrow = 2, ncol = n.traits,
                          dimnames = list(c("xlim.min", "xlim.max"), trait))
      data.xlim[1, ] <- apply(X, 2, min, na.rm = TRUE)
      data.xlim[1, apply(X, 2, min, na.rm = TRUE) * apply(X, 2, max, na.rm = TRUE) > 0 & apply(X, 2, min, na.rm = TRUE) > 0] <- 0
      data.xlim[2, ] <- apply(X, 2, max, na.rm = TRUE)
      data.xlim[2, apply(X, 2, min, na.rm = TRUE) * apply(X, 2, max, na.rm = TRUE) > 0 & apply(X, 2, max, na.rm = TRUE) < 0] <- 0
      if(!is.null(error.bar.inf) & !is.null(error.bar.sup)){
        data.xlim[1, ] <- apply(cbind(apply(arrow.inf, 2, min, na.rm = TRUE), apply(arrow.sup, 2, min)), 1, min, na.rm = TRUE)
        data.xlim[2, ] <- apply(cbind(apply(arrow.inf, 2, max, na.rm = TRUE), apply(arrow.sup, 2, max)), 1, max, na.rm = TRUE)
      } else {
        if(!is.null(error.bar.inf)){
          data.xlim[1, ] <- apply(cbind(apply(arrow.inf, 2, min, na.rm = TRUE), data.xlim[1, ]), 1, min, na.rm = TRUE)
          data.xlim[2, ] <- apply(cbind(apply(arrow.inf, 2, max, na.rm = TRUE), data.xlim[2, ]), 1, max, na.rm = TRUE)
        }
        if(!is.null(error.bar.sup)){
          data.xlim[1, ] <- apply(cbind(apply(arrow.sup, 2, min, na.rm = TRUE), data.xlim[1, ]), 1, min, na.rm = TRUE)
          data.xlim[2, ] <- apply(cbind(apply(arrow.sup, 2, max, na.rm = TRUE), data.xlim[2, ]), 1, max, na.rm = TRUE)
        }
      }
    }
  } else if(is.vector(data.xlim) & length(data.xlim) == 2){
    data.xlim <- matrix(rep(data.xlim, n.traits), nrow = 2,
                        dimnames = list(c("xlim.min", "xlim.max"), trait))
  } else if(is.matrix(data.xlim)){
    if(isTRUE(all.equal(dim(data.xlim), c(2, n.traits)))){
      rownames(data.xlim) <- c("xlim.min", "xlim.max")
      colnames(data.xlim) <- trait
    } else{
      stop("Invalid 'data.xlim' argument: wrong matrix dimensions")
    }
  } else {
    stop("Invalid 'data.xlim' argument")
  }
  
  
  ylim <- c(1, n.tips)
  
  if(plot.type == "barplot"){
    bar.col <- .orderGrArg(bar.col, n.tips = n.tips, n.traits = n.traits,
                           new.order = new.order, tips = tips, default = "grey35")
  }
  
  if(plot.type == "dotplot"){
    
    dot.col <- .orderGrArg(dot.col, n.tips = n.tips, n.traits = n.traits,
                           new.order = new.order, tips = tips, default = 1)
    
    dot.pch <- .orderGrArg(dot.pch, n.tips = n.tips, n.traits = n.traits,
                           new.order = new.order, tips = tips, default = 1)
    
    dot.cex <- .orderGrArg(dot.cex, n.tips = n.tips, n.traits = n.traits,
                           new.order = new.order, tips = tips, default = 1)
    
  }
  
  if(!is.null(error.bar.inf) | !is.null(error.bar.sup)){
    error.bar.col <- .orderGrArg(error.bar.col, n.tips = n.tips, n.traits = n.traits,
                                 new.order = new.order, tips = tips, default = 1)
  }
  
  if(is.null(tip.labels)){
    tip.labels <- tips
  } else {
    tip.labels <- .orderGrArg(tip.labels, n.tips = n.tips, n.traits = n.traits,
                              new.order = new.order, tips = tips, default = "")
  }
  tip.col <- .orderGrArg(tip.col, n.tips = n.tips, n.traits = n.traits,
                         new.order = new.order, tips = tips, default = 1)
  tip.cex <- .orderGrArg(tip.cex, n.tips = n.tips, n.traits = n.traits,
                         new.order = new.order, tips = tips, default = 1)
  tip.font <- .orderGrArg(tip.font, n.tips = n.tips, n.traits = n.traits,
                          new.order = new.order, tips = tips, default = 3)
  
  
  
  if(is.null(trait.labels)){
    trait.labels <- trait
  }
  trait.labels <- rep(trait.labels, length.out = n.traits)
  trait.col <- rep(trait.col, length.out = n.traits)
  trait.cex <- rep(trait.cex, length.out = n.traits)
  trait.font <- rep(trait.font, length.out = n.traits)
  trait.bg.col <- rep(trait.bg.col, length.out = n.traits)
  
  if(is.null(tree.xlim)){
    tree.xlim <- .plotPhyloDisabled(phy, type = tree.type,
                                    show.tip.label = FALSE,
                                    x.lim = NULL, y.lim = NULL,
                                    no.margin = FALSE, direction = "rightwards",
                                    plot = FALSE, ...)$x.lim
  }
  
  ###
  par.mar0 <- par("mar")
  par.lend0 <- par("lend")
  par.xpd0 <- par("xpd")
  
  if(tree.type == "phylogram" | tree.type == "cladogram"){
    
    if(plot.type %in% c("barplot", "dotplot")){
      lay <- .layouterize(n.traits = n.traits, show.tip = show.tip)
      lay.w <- .layouterizeRatio(tree.ratio = tree.ratio, n.traits = n.traits, show.tip = show.tip)
    }
    if(plot.type == "gridplot"){
      lay <- .layouterize(n.traits = 1, show.tip = show.tip)
      lay.w <- .layouterizeRatio(tree.ratio = tree.ratio, n.traits = 1, show.tip = show.tip)
    }
    
    layout(lay, widths = lay.w)
    par(xpd = FALSE, mar=c(5, 1, 4, 0), lend = 1)
    fig.traits <- vector("list", n.traits)
    names(fig.traits) <- trait
    
    if(plot.type %in% c("barplot", "dotplot")){
      for(i in 1:n.traits){
        plot.new()
        plot.window(xlim = data.xlim[, i], ylim = ylim)
        fig.traits[[i]] <- par("fig")
        rect(par("usr")[1], par("usr")[3]-(3 * par("cxy")[2]), par("usr")[2], par("usr")[4],
             col = trait.bg.col[i], border = NA, xpd = TRUE)
        if(show.box){
          box()
        }
        if(grid.vertical){        
          grid(NULL, NA, col = grid.col, lty = grid.lty)
          abline(v = 0, lty = "solid", col = grid.col)
        } else {
          abline(v = 0, lty = "solid", col = grid.col)
        }
        if(grid.horizontal){
          abline(h= 1:n.tips, col = grid.col, lty = grid.lty)
        }
        if(plot.type == "barplot"){
          segments(x0 = 0, x1 = X[, i], y0 = 1:n.tips, lwd = bar.lwd, col = bar.col[, i])
        }
        if(plot.type == "dotplot"){
          points(x = X[, i], y = 1:n.tips, col = dot.col[, i], pch= dot.pch[, i], cex = dot.cex[, i])
        }
        options(warn = -1)
        if(!is.null(error.bar.inf)){
          arrows(x0 = X[, i], x1 = arrow.inf[, i], y0 = 1:n.tips,
                 lwd = 1, col = error.bar.col[, i], angle = 90, length = 0.04)
        }
        if(!is.null(error.bar.sup)){
          arrows(x0 = X[, i], x1 = arrow.sup[, i], y0 = 1:n.tips,
                 lwd = 1, col = error.bar.col[, i], angle = 90, length = 0.04)
        }
        options(warn = 1)
        if(show.data.axis){
          axis(1)
        }
        if(show.trait){
          mtext(trait.labels[i], side = 1, line = 3, las = par("las"),
                col = trait.col[i], cex = trait.cex[i],
                font = trait.font[i])
        }
      }
    }
    if(plot.type == "gridplot"){
      plot.new()
      rect(par("usr")[1], par("usr")[3] - (3 * par("cxy")[2]), par("usr")[2], par("usr")[4],
           col = trait.bg.col[1], border = NA, xpd = TRUE)
      data.xlim[1, ] <- 0
      data.xlim[2, ] <- n.traits
      plot.window(xlim = data.xlim[, 1], ylim = ylim)
      fig.traits[[1]] <- par("fig")
      image(x = 0:n.traits, y = 1:n.tips, z = t(X),
            col = cell.col, add = TRUE,
            xlab = "", ylab = "", yaxs = FALSE, xaxs = FALSE)
      
      if(show.box){
        box()
      }
      
      if(grid.horizontal){        
        abline(h = seq(1.5, n.tips - 0.5), col = grid.col, lty = grid.lty)
      }
      if(grid.vertical){
        abline(v = seq(1, n.traits - 1), col = grid.col, lty = grid.lty)
      }
      
      if(show.trait){
        mtext(trait.labels, at = seq(0.5, (n.traits - 0.5)),
              side = 1, line = 1, las = par("las"),
              col = trait.col, cex = trait.cex,
              font = trait.font)    
      }
    }
    
    if(show.tip){
      plot.new()
      tip.xlim <- c(-1, 1)
      if(tip.adj < 0.5) tip.xlim[1] <- -tip.adj / 0.5
      if(tip.adj > 0.5) tip.xlim[2] <- -2 * tip.adj + 2
      plot.window(xlim = tip.xlim, ylim = ylim)
      text(x = 0, y = 1:n.tips, labels = tip.labels,
           adj = tip.adj, col = tip.col, cex = tip.cex, font = tip.font)
      fig.tip <- par("fig")
    } else {
      fig.tip <- NULL
      tip.xlim <- NULL
    }
    
    plot.phylo(phy, type = tree.type, show.tip.label = FALSE,
               x.lim = tree.xlim, y.lim = NULL,
               no.margin = FALSE, direction = "rightwards", ...)
    fig.tree <- par("fig")
    
    if(plot.type == "gridplot" & show.color.scale){
      par(new = TRUE)
      plt.init <- par("plt")
      par(plt = c(par("plt")[1] + 0.05, par("plt")[2] - 0.2, 0.07, 0.1))
      plot.new()
      breaks <- seq(min(X), max(X), length.out = (length(cell.col) + 1))
      scale.xlim <- range(breaks)
      scale.ylim <- c(0, 1)
      plot.window(xlim = scale.xlim, ylim = scale.ylim)
      for(i in 1:length(cell.col)){
        polygon(c(breaks[i], breaks[i + 1], breaks[i + 1], breaks[i]), c(0, 0, 1, 1),
                col = cell.col[i], border = NA)
      }
      axis(1)
      par(plt = plt.init)
    }
    
    assign("last_barplotp4d", list(plot.type = plot.type,
                                   show.tip = show.tip,
                                   layout = lay,
                                   fig.tree = fig.tree,
                                   fig.traits = fig.traits,
                                   fig.tip = fig.tip,
                                   tree.xlim = tree.xlim,
                                   data.xlim = data.xlim,
                                   tip.xlim = tip.xlim,
                                   ylim = ylim, par.mar0 = par.mar0), 
           envir = ape::.PlotPhyloEnv)
    layout(1)
  }
  
  
  if(tree.type == "fan"){
    
    par(lend = 1)
    if(is.null(tree.ratio)){
      if(show.tip){
        tree.ratio <- 1 / (n.traits + 2)
      } else {
        tree.ratio <- 1 / (n.traits + 1)
      }
    }
    
    plot.phylo(phy, type = tree.type, show.tip.label = FALSE,
               x.lim = tree.xlim * (1/tree.ratio), y.lim = NULL,
               no.margin = TRUE, open.angle = tree.open.angle, rotate.tree = 0, ...)
    lp <- get("last_plot.phylo", envir = ape::.PlotPhyloEnv)
    
    length.phylo <- max(sqrt(lp$xx^2 + lp$yy^2))
    if(show.tip){
      length.gr0 <- (min(par("usr")[2] - par("usr")[1], par("usr")[4] - par("usr")[3]) / 2 - length.phylo) / (n.traits+1)
    } else {
      length.gr0 <- (min(par("usr")[2] - par("usr")[1], par("usr")[4] - par("usr")[3]) / 2 - length.phylo) / n.traits
    }
    length.intergr <- 0.2 * length.gr0
    length.gr <- length.gr0 - length.intergr
    
    theta <- atan2(lp$xx[1:n.tips], lp$yy[1:n.tips])[new.order]
    theta[theta > (pi/2)] <- -pi - (pi - theta[theta > (pi/2)])
    cos.t <- cos(pi/2 - theta)
    sin.t <- sin(pi/2 - theta)
    
    theta.real.open <- diff(c(min(theta), max(theta)))
    real.open <- theta.real.open * 180/pi
    
    if(tree.open.crown){
      theta.soft <- pi/2 - seq(-5, real.open+5, length.out=300) * pi/180
    } else {
      theta.soft <- pi/2 - seq(0, 360) * pi/180 + 10 * pi/180
    }
    cos.tsoft <- cos(pi/2 - theta.soft)
    sin.tsoft <- sin(pi/2 - theta.soft)
    
    for(i in 1:n.traits){
      # RINGS
      length.ring1 <- length.phylo + length.intergr*i + length.gr*(i-1) - 0.3*length.intergr
      length.ring2 <- length.phylo + length.intergr*i + length.gr*i + 0.3*length.intergr
      xx1 <- length.ring1 * cos.tsoft
      xx2 <- length.ring2 * cos.tsoft
      yy1 <- length.ring1 * sin.tsoft
      yy2 <- length.ring2 * sin.tsoft
      polygon(c(xx1, rev(xx2)), c(yy1, rev(yy2)), col = trait.bg.col[i], border = NA)
      
      # SCALING VALUES
      if(abs(sign(min(data.xlim[, i])) + sign(max(data.xlim[, i]))) == 2){
        scaling.factor <- length.gr / max(abs(min(data.xlim[, i])), abs(max(data.xlim[, i])))
      } else {
        scaling.factor <- length.gr / diff(c(min(data.xlim[, i]), max(data.xlim[, i])))
      }
      X.scale <- X[, i] * scaling.factor
      data.xlim.scale <- data.xlim[, i] * scaling.factor
      
      if(!is.null(error.bar.inf)){
        arrow.inf.scale <- arrow.inf[, i] * scaling.factor
      }
      if(!is.null(error.bar.sup)){
        arrow.sup.scale <- arrow.sup[, i] * scaling.factor
      }
      
      # FOR BARPLOT AND DOTPLOT
      # Baseline and Values
      if(plot.type == "barplot" | plot.type == "dotplot"){
        length.baseline <- (length.phylo + length.intergr*i + length.gr*(i-1) +
                              ifelse(min(data.xlim.scale) < 0, abs(min(data.xlim.scale)), 0))
        length.baseline <- rep(length.baseline, n.tips)
        length.values <- length.baseline + X.scale
        if(!is.null(error.bar.inf)){
          length.arrow.inf <- length.baseline + arrow.inf.scale
        }
        if(!is.null(error.bar.sup)){
          length.arrow.sup <- length.baseline + arrow.sup.scale
        }
        
        # Draw Baseline
        length.baseline <- rep(length.baseline, length.out = length(cos.tsoft))
        lines(length.baseline * cos.tsoft, length.baseline * sin.tsoft, lwd = 1)
        
        #Grid Species
        if(grid.horizontal){
          segments(x0 = length.ring1 * cos.t,
                   x1 = length.ring2 * cos.t,
                   y0 = length.ring1 * sin.t,
                   y1 = length.ring2 * sin.t,
                   col = grid.col, lty = grid.lty)
        }
        
        # Axis and circular grid
        if((show.data.axis | grid.vertical)){
          if(tree.open.crown){
            theta.ax <- theta.soft[1]
          } else {
            if(show.trait){
              theta.ax <-theta[1] + (360-real.open)*(1/3) * pi/180
            } else {
              theta.ax <-theta[1] + (360-real.open)*(1/2) * pi/180
            }
          }
          cos.tax <- cos(pi/2 - theta.ax)
          sin.tax <- sin(pi/2 - theta.ax)
          nint.ticks <- round((length.gr / min(par("usr")[2] - par("usr")[1], par("usr")[4] - par("usr")[3]))/3 *100) - 1
          if(min(data.xlim.scale) <= 0 & max(data.xlim.scale) >=0){
            ticks <- axisTicks(c(min(data.xlim.scale)/scaling.factor,
                                 max(data.xlim.scale)/scaling.factor), log = FALSE, nint = nint.ticks)
          } else {
            if(abs(min(data.xlim.scale)) > max(data.xlim.scale)){
              ticks <- axisTicks(c(0, min(data.xlim.scale)/scaling.factor), log = FALSE, nint = nint.ticks)
            } else {
              ticks <- axisTicks(c(0, max(data.xlim.scale)/scaling.factor), log = FALSE, nint = nint.ticks)
            }
          }
          ticks <- ifelse(ticks > max(data.xlim.scale)/scaling.factor, NA, ticks)
          
          length.ticks <- length.baseline[1] + ticks * scaling.factor
          
          # Circular Grid
          if(grid.vertical){
            for(j in 1:length(length.ticks)){
              lines(length.ticks[j] * cos.tsoft,
                    length.ticks[j] * sin.tsoft,
                    col = grid.col, lty = grid.lty)
            }
          }
          
          if(show.data.axis){
            # Frame for axis values
            segments(x0 = length.ring1 * cos.tax,
                     x1 = length.ring2 * cos.tax,
                     y0 = length.ring1 * sin.tax,
                     y1 = length.ring2 * sin.tax,
                     lwd = 20, col = trait.bg.col[i])
            
            # Axis Values
            text(x = length.ticks * cos.tax, 
                 y = length.ticks * sin.tax,
                 labels = ticks, cex = tip.cex)
          }
        }
      }
      
      
      
      # Draw Values
      if(plot.type == "barplot"){
        length.baseline <- rep(length.baseline, length.out = length(cos.t))
        segments(x0 = length.baseline * cos.t,
                 x1 = length.values * cos.t,
                 y0 = length.baseline * sin.t,
                 y1 = length.values * sin.t,
                 lwd = bar.lwd, col = bar.col[, i])
      }
      if(plot.type == "dotplot"){
        length.baseline <- rep(length.baseline, length.out = length(cos.t))
        points(x = length.values * cos.t,
               y = length.values * sin.t,
               col = dot.col[, i],
               pch= dot.pch[, i],
               cex = dot.cex[, i])
      }
      
      # GRIDPLOT
      if(plot.type == "gridplot"){
        nc <- length(cell.col)
        X.cut <- as.numeric(cut(as.matrix(X), nc))
        grid.colors <- cell.col[X.cut]
        grid.colors <- matrix(grid.colors, ncol = n.traits)
        
        theta.grid1 <- theta + pi/2 - (theta[1] + theta[2])/2
        theta.grid2 <- theta - pi/2 + (theta[1] + theta[2])/2
        
        for(k in 1:length(theta)){
          tile.x1 <- length.ring1 * cos(pi/2 - seq(theta.grid1[k], theta.grid2[k], length.out = 25))
          tile.x2 <- length.ring2 * cos(pi/2 - seq(theta.grid1[k], theta.grid2[k], length.out = 25))
          tile.y1 <- length.ring1 * sin(pi/2 - seq(theta.grid1[k], theta.grid2[k], length.out = 25))
          tile.y2 <- length.ring2 * sin(pi/2 - seq(theta.grid1[k], theta.grid2[k], length.out = 25))
          polygon(c(tile.x1, rev(tile.x2)), c(tile.y1, rev(tile.y2)), col = grid.colors[k, i] , border = NA)
        }
        if(grid.horizontal){
          segments(x0 = length.ring1 * cos(pi/2 - theta.grid1),
                   y0 = length.ring1 * sin(pi/2 - theta.grid1),
                   x1 = length.ring2 * cos(pi/2 - theta.grid1),
                   y1 = length.ring2 * sin(pi/2 - theta.grid1),
                   col = grid.col, lty = grid.lty)
          segments(x0 = length.ring1 * cos(pi/2 - theta.grid2[length(theta.grid2)]),
                   y0 = length.ring1 * sin(pi/2 - theta.grid2[length(theta.grid2)]),
                   x1 = length.ring2 * cos(pi/2 - theta.grid2[length(theta.grid2)]),
                   y1 = length.ring2 * sin(pi/2 - theta.grid2[length(theta.grid2)]),
                   col = grid.col, lty = grid.lty)
        }
        if(grid.vertical){
          if(i > 1){
            lines((length.ring1 - length.intergr/2 + 0.3*length.intergr) * cos.tsoft,
                  (length.ring1 - length.intergr/2 + 0.3*length.intergr) * sin.tsoft,
                  col = grid.col, lty = grid.lty)
          }
        }
      }
      
      
      # Draw Error Bars
      if(plot.type == "barplot" | plot.type == "dotplot"){
        options(warn = -1)
        if(!is.null(error.bar.inf)){
          arrows(x0 = length.values * cos.t,
                 x1 = length.arrow.inf * cos.t,
                 y0 = length.values * sin.t,
                 y1 = length.arrow.inf * sin.t,
                 lwd = 1, col = error.bar.col[, i],
                 angle = 90, length = 0.04)
        }
        
        if(!is.null(error.bar.sup)){
          arrows(x0 = length.values * cos.t,
                 x1 = length.arrow.sup * cos.t,
                 y0 = length.values * sin.t,
                 y1 = length.arrow.sup * sin.t,
                 lwd = 1, col = error.bar.col[, i],
                 angle = 90, length = 0.04)
        }
        options(warn = 1)
      }
      
      # Draw Box
      if(show.box){
        if(tree.open.crown){
          lines(c(xx1, rev(xx2)), c(yy1, rev(yy2)), col = 1)
        } else {
          lines(xx1, yy1, col = 1)
          lines(xx2, yy2, col = 1)          
        } 
      }
      
      # TRAITS LABELS
      if(show.trait){
        if(tree.open.crown){
          theta.trait <- theta.soft[length(theta.soft)]
        } else {
          if(show.data.axis & (plot.type == "barplot" | plot.type == "dotplot")){
            theta.trait <-theta[length(theta)] - (360-real.open)*(1/3) * pi/180
          } else {
            theta.trait <-theta[length(theta)] - (360-real.open)*(1/2) * pi/180
          }
        }
        cos.ttrait <- cos(pi/2 - theta.trait)
        sin.ttrait <- sin(pi/2 - theta.trait)
        # Frame for Traits Labels
        segments(x0 = length.ring1 * cos.ttrait,
                 x1 = length.ring2 * cos.ttrait,
                 y0 = length.ring1 * sin.ttrait,
                 y1 = length.ring2 * sin.ttrait,
                 lwd = 20, col = trait.bg.col[i])
        
        # Traits Labels
        text(x = (length.ring1 + length.ring2)/2 * cos.ttrait, 
             y = (length.ring1 + length.ring2)/2 * sin.ttrait,
             labels = trait.labels[i],
             col = trait.col[i], cex = trait.cex[i],
             font = trait.font[i],
             srt = ifelse(theta.trait > 0 | theta.trait < -pi, 
                          (pi/2 - theta.trait) * 180 /pi, (-pi/2 - theta.trait) * 180 /pi))
      }
      
    }
    
    # TIPS
    if(show.tip){
      length.tipsline <- (length.phylo+length.intergr*(n.traits+1)+length.gr*(n.traits))
      tip.xlim <- c(-1, 1)
      if(tip.adj < 0.5) tip.xlim[1] <- -tip.adj / 0.5
      if(tip.adj > 0.5) tip.xlim[2] <- -2 * tip.adj + 2
      
      for(i in 1:n.tips){
        text(x = length.tipsline * cos.t[i],
             y = length.tipsline * sin.t[i],
             labels = tip.labels[i],
             adj = ifelse(theta[i] > 0 | theta[i] < -pi, 0, 1), 
             col = tip.col[i], cex = tip.cex[i], font = tip.font[i],
             srt = ifelse(theta[i] > 0 | theta[i] < -pi, 
                          (pi/2 - theta[i]) * 180 /pi, (-pi/2 - theta[i]) * 180 /pi))
      }
    }
    
  }
  
  par(mar = par.mar0, xpd = par.xpd0, lend = par.lend0)
  invisible()
}



#' Focus on sub parts of a plot
#'
#' These functions can be used after \code{\link{barplot.phylo4d}}, \code{\link{dotplot.phylo4d}}
#' and \code{\link{gridplot.phylo4d}} when
#' \code{tree.type} is \code{"phylogram"} or \code{"cladogram"} to focus on
#' the different part of the plot and add graphical elements.
#' 
#' @param x the trait to focus on.
#' Can be a character string giving the name of the trait or an integer giving
#' the number of the trait in order of appearance in the plot.
#' 
#' @details #' Use \code{focusTree} to focus on the phylogenetic tree, \code{focusTraits}
#' to focus on a given trait and \code{focusTips} to focus on the tips labels.
#' Use \code{focusStop} to close the editing and restore graphical settings.
#' For each part of the plot, the coordinate system is restored, making edition easier.
#' For the phylogeny, post-editing functions of the package \code{ape} like \code{nodelabels} can be used.
#'
#' @examples
#' require(ape)
#' require(phylobase)
#' data(navic)
#' dat <- tdata(navic)
#' neidium.sp <- c("Neidium bisulcatum",
#'                 "Neidium affine",
#'                 "Neidium productum")
#' stauroneis.sp <- c("Stauroneis kriegeri",
#'                    "Stauroneis acuta",
#'                    "Stauroneis gracilior",
#'                    "Stauroneis phoenicenteron")
#' neidium.mean <- mean(dat[neidium.sp,])
#' stauroneis.mean <- mean(dat[stauroneis.sp, ])
#' 
#' dotplot(navic, center = FALSE, scale = FALSE, data.xlim= c(0, 6))
#' 
#' focusTree()
#' nodelabels(node=c(22, 32), pch = 20, cex = 3, col = c(2, 3))
#' 
#' focusTraits()
#' segments(x0 = neidium.mean, y0 = 14.5,
#'          x1 = neidium.mean, y1 = 17.5,
#'          col = 3, lty = "dashed", lwd = 2)
#' segments(x0 = stauroneis.mean, y0 = 2.5,
#'          x1 = stauroneis.mean, y1 = 7.5,
#'          col = 2, lty = "dashed", lwd = 2)
#' 
#' focusTips()
#' rect(xleft = 0, ybottom = 2.5,
#'      xright = 0.9, ytop = 7.5,
#'      col = "#FF000020", border = NA)
#' rect(xleft = 0, ybottom = 14.5,
#'      xright = 0.9, ytop = 17.5,
#'      col = "#00FF0020", border = NA)
#'
#' focusStop()
#'
#'@rdname focus
#'@export
focusTraits <- function(x){
  lp <- get("last_barplotp4d", envir = ape::.PlotPhyloEnv)
  if(lp$plot.type == "gridplot"){
    x <- 1
  }
  par(new = TRUE)
  plot.new()
  layout(lp$layout)
  par(mar = c(5, 1, 4, 0))
  fig <- unlist(lp$fig.traits[x])
  fig[fig < 0] <- 0
  par(fig = fig)
  plot.window(xlim = lp$data.xlim[, x], ylim = lp$ylim)
}

#'@rdname focus
#'@export
focusTree <- function(){
  lp <- get("last_barplotp4d", envir = ape::.PlotPhyloEnv)
  par(new = TRUE)
  plot.new()
  layout(lp$layout)
  par(mar = c(5, 1, 4, 0))
  fig <- unlist(lp$fig.tree)
  fig[fig < 0] <- 0
  par(fig = fig)
  plot.window(xlim = lp$tree.xlim, ylim = lp$ylim)
}

#'@rdname focus
#'@export
focusTips <- function(){
  lp <- get("last_barplotp4d", envir = ape::.PlotPhyloEnv)
  if(!lp$show.tip){
    stop("No tip labels on the figure")
  } else {
    par(new = TRUE)
    plot.new()
    layout(lp$layout)
    par(mar = c(5, 1, 4, 0))
    fig <- unlist(lp$fig.tip)
    fig[fig < 0] <- 0
    par(fig = fig)
    plot.window(xlim = lp$tip.xlim, ylim = lp$ylim)
  }
}

#'@rdname focus
#'@export
focusStop <- function(){
  lp <- get("last_barplotp4d", envir = ape::.PlotPhyloEnv)
  layout(1)
  par(mar = lp$par.mar0)
}




#' Barplot of Traits Values along a Phylogeny
#' 
#' @inheritParams multiplot.phylo4d
#' @param height a \code{phylo4d} object.
#' 
#' @examples
#' data(navic)
#' barplot(navic)
#' 
#' @method barplot phylo4d
#' @export
barplot.phylo4d <- function(height, trait = names(tdata(height)), center = TRUE, scale = TRUE, 
                            tree.ladderize = FALSE, tree.type = "phylogram", tree.ratio = NULL, tree.xlim = NULL,
                            tree.open.angle = 0, tree.open.crown = TRUE,
                            show.tip = TRUE, tip.labels = NULL, tip.col = "black", tip.cex = 1, tip.font = 3, tip.adj = 0,
                            data.xlim = NULL, bar.lwd = 10, bar.col = "grey35", show.data.axis = TRUE,
                            show.trait = TRUE, trait.labels = NULL, trait.col = "black", trait.cex = 1, trait.font = 1,
                            trait.bg.col = "grey90", error.bar.sup = NULL, error.bar.inf = NULL, error.bar.col = 1,
                            show.box = FALSE, grid.vertical = TRUE, grid.horizontal = FALSE, grid.col = "grey25",
                            grid.lty = "dashed", ...){
  
  multiplot.phylo4d(p4d = height, trait = trait, center = center, scale = scale, plot.type = "barplot",
                    tree.ladderize = tree.ladderize, tree.type = tree.type, tree.ratio = tree.ratio, tree.xlim = tree.xlim,
                    tree.open.angle = tree.open.angle, tree.open.crown = tree.open.crown,
                    show.tip = show.tip, tip.labels = tip.labels, tip.col = tip.col, tip.cex = tip.cex, tip.font = tip.font, tip.adj = tip.adj,
                    data.xlim = data.xlim, bar.lwd = bar.lwd, bar.col = bar.col, show.data.axis = show.data.axis,
                    show.trait = show.trait, trait.labels = trait.labels, trait.col = trait.col, trait.cex = trait.cex, trait.font = trait.font,
                    trait.bg.col = trait.bg.col, error.bar.sup = error.bar.sup, error.bar.inf = error.bar.inf, error.bar.col = error.bar.col,
                    show.box = show.box, grid.vertical = grid.vertical, grid.horizontal = grid.horizontal, grid.col = grid.col,
                    grid.lty = grid.lty, ...)
  
}



#' Dotplot of Traits Values along a Phylogeny
#' 
#' @inheritParams multiplot.phylo4d
#' @examples
#' data(navic)
#' dotplot(navic)
#' 
#' @export
dotplot.phylo4d <- function(p4d, trait = names(tdata(p4d)), center = TRUE, scale = TRUE,
                            tree.ladderize = FALSE, tree.type = "phylogram", tree.ratio = NULL, tree.xlim = NULL,
                            tree.open.angle = 0, tree.open.crown = TRUE,
                            show.tip = TRUE, tip.labels = NULL, tip.col = "black", tip.cex = 1, tip.font = 3, tip.adj = 0,
                            data.xlim = NULL, show.data.axis = TRUE,
                            dot.col = "black", dot.pch = 20, dot.cex = 2,
                            show.trait = TRUE, trait.labels = NULL, trait.col = "black", trait.cex = 1, trait.font = 1,
                            trait.bg.col = "grey90", error.bar.sup = NULL, error.bar.inf = NULL, error.bar.col = 1,
                            show.box = FALSE, grid.vertical = FALSE, grid.horizontal = TRUE, grid.col = "grey25",
                            grid.lty = "dashed", ...){
  
  multiplot.phylo4d(p4d, trait = trait, center = center, scale = scale, plot.type = "dotplot",
                    tree.ladderize = tree.ladderize, tree.type = tree.type, tree.ratio = tree.ratio, tree.xlim = tree.xlim,
                    tree.open.angle = tree.open.angle, tree.open.crown = tree.open.crown,
                    show.tip = show.tip, tip.labels = tip.labels, tip.col = tip.col, tip.cex = tip.cex, tip.font = tip.font, tip.adj = tip.adj,
                    data.xlim = data.xlim, show.data.axis = show.data.axis,
                    dot.col = dot.col, dot.pch = dot.pch, dot.cex = dot.cex,
                    show.trait = show.trait, trait.labels = trait.labels, trait.col = trait.col, trait.cex = trait.cex, trait.font = trait.font,
                    trait.bg.col = trait.bg.col, error.bar.sup = error.bar.sup, error.bar.inf = error.bar.inf, error.bar.col = error.bar.col,
                    show.box = show.box, grid.vertical = grid.vertical, grid.horizontal = grid.horizontal, grid.col = grid.col,
                    grid.lty = grid.lty, ...)
  
}

#' Dotplot
#' 
#' @param ... further arguments passed to or from other methods.
#' @export
dotplot <- function(...){
  UseMethod("dotplot")
}


#' Gridplot of Traits Values along a Phylogeny
#' @inheritParams multiplot.phylo4d
#' 
#' @examples
#' data(navic)
#' gridplot(navic)
#' 
#' # Multivariate data
#' require(phylobase)
#' tipData(navic) <- matrix(rnorm(170), nrow = 17)
#' gridplot(navic)
#' 
#' @export
gridplot.phylo4d <- function(p4d, trait = names(tdata(p4d)), center = TRUE, scale = TRUE,
                             tree.ladderize = FALSE, tree.type = "phylogram", tree.ratio = NULL, tree.xlim = NULL,
                             tree.open.angle = 0, tree.open.crown = TRUE,
                             show.tip = TRUE, tip.labels = NULL, tip.col = "black", tip.cex = 1, tip.font = 3, tip.adj = 0,
                             cell.col = white2red(100), show.color.scale = TRUE,
                             show.trait = TRUE, trait.labels = NULL, trait.col = "black", trait.cex = 0.7, trait.font = 1,
                             trait.bg.col = "grey90",
                             show.box = FALSE, grid.vertical = FALSE, grid.horizontal = FALSE, grid.col = "grey25",
                             grid.lty = "dashed", ...){
  
  multiplot.phylo4d(p4d, trait = trait, center = center, scale = scale, plot.type = "gridplot",
                    tree.ladderize = tree.ladderize, tree.type = tree.type, tree.ratio = tree.ratio, tree.xlim = tree.xlim,
                    tree.open.angle = tree.open.angle, tree.open.crown = tree.open.crown,
                    show.tip = show.tip, tip.labels = tip.labels, tip.col = tip.col, tip.cex = tip.cex, tip.font = tip.font, tip.adj = tip.adj,
                    cell.col = cell.col, show.color.scale = show.color.scale,
                    show.trait = show.trait, trait.labels = trait.labels, trait.col = trait.col, trait.cex = trait.cex, trait.font = trait.font,
                    trait.bg.col = trait.bg.col,
                    show.box = show.box, grid.vertical = grid.vertical, grid.horizontal = grid.horizontal, grid.col = grid.col,
                    grid.lty = grid.lty, ...)
  
}

#' Gridplot
#' 
#' @param ... further arguments passed to or from other methods.
#' @export
gridplot <- function(...){
  UseMethod("gridplot")
}