#' Make a phylomorphospace plot

#Made by W.Andrew Barr, modified by V Fischer in 2018-2019 to add density and some plotting options

#' @param tree object of class \code{phylo}
#' @param tipinfo a dataframe containing two continuous variables (typically PC1 and PC2 from a PCA) as well as a grouping variable and a variable of labels for each point on the plot
#' @param xvar the unquoted name of the column containing the x variable in the tipinfo dataframe. Defaults to \code{PC1}
#' @param yvar the unquoted name of the column containing the y variable in the tipinfo dataframe. Defaults to \code{PC2}
#' @param factorvar the unquoted name of the column containing the grouping variable in the tipinfo dataframe. Defaults to \code{group}
#' @param labelvar the unquoted name of the column containing the labeling variable in the tipinfo dataframe. Defaults to \code{taxon}
#' @param title the text string to use as the title of the phylomorphospace plot
#' @param xlab the text label for the x axis
#' @param ylab the text label for the y axis
#' @param repel boolean indicating whether or not to use the \code{ggrepel} package to prevent overplotting of labels
#' @param edge.width size parameter for edges drawn on phylomorphospace. Default is 1
#' @param fontface font face for labels. Default is "italic".
#' @param tree.alpha alpha parameter for the main tree. default is 0.7
#' @return the ggplot object representing the phylomorphospace


ggphylomorpho <- function(tree,

                          tipinfo,

                          xvar=PC1,

                          yvar=PC2,
                          
                          sizedata=1,

                          factorvar=group,

                          labelvar=taxon,

                          title="Phylomorphospace",

                          xlab="Axis1",

                          ylab="Axis2",
                          
                          names=TRUE, #whether to show taxon names or not

                          repel=TRUE, #spreads taxon names

                          edge.width=1,

                          fontface="italic",

                          tree.alpha = 0.7,
                          
                          h=1,
                          
                          lim=10)

  {


  require(ggplot2)
  require(phytools)
  require(ggrepel)


  ## create matrix for use in phytools::fastAnc()

  mat <- cbind(eval(substitute(xvar), tipinfo),eval(substitute(yvar), tipinfo))
  rownames(mat) <- eval(substitute(labelvar), tipinfo)
  stopifnot(length(setdiff(tree$tip.label, rownames(mat))) == 0)

  xAnc <- fastAnc(tree, mat[,1])
  yAnc <- fastAnc(tree, mat[,2])

  all_node_coords <-
    data.frame(
      #put PC values for all nodes and tips in a dataframe
      #tips go first in order of tip labels, then numerical order for nodes
      x=c(mat[tree$tip.label,1], xAnc),
      y=c(mat[tree$tip.label,2], yAnc),
      nodeid=1:(tree$Nnode + length(tree$tip.label))
      )

  #get edge list from tree object
  edges <- data.frame(tree$edge)
  names(edges) <- c("node1", "node2")
  #translate tip/node numbers into PC coordinates in all_node_coords dataframe
  edgecoords <- merge(
    merge(edges, all_node_coords, by.x="node1", by.y="nodeid"),
    all_node_coords, by.x="node2", by.y="nodeid")

  pointsForPlot <-
    data.frame(x=eval(substitute(xvar), tipinfo),
               y=eval(substitute(yvar), tipinfo),
               color=eval(substitute(factorvar), tipinfo),
               label=eval(substitute(labelvar), tipinfo))

  theplot <-
    ggplot() +
    stat_density_2d(data = pointsForPlot,aes(x =x, y =y,alpha=..level..),h=h,geom="polygon")+
    scale_alpha(guide="none")+
    geom_segment(data=edgecoords,aes(x=x.x,xend=x.y, y=y.x, yend=y.y), size=edge.width, alpha=tree.alpha) +
    geom_point(data=pointsForPlot, aes(x=x, y=y,color=color,size=log(sizedata))) +
    expand_limits(x=c(-lim,lim),y=c(-lim,lim))+
    coord_fixed(ratio=1)
if (names){
  if(repel){
    theplot <- theplot + geom_text_repel(data=pointsForPlot, aes(x=x, y=y, label=label), segment.alpha=0.5, fontface=fontface)
  } else{
    theplot <- theplot + geom_text(data=pointsForPlot, aes(x=x, y=y, label=label), fontface=fontface)
  }
 }
  return(theplot)
}