#' Testing the effect of sampling on clade-associated and character-associated diversification patterns
#'
#' usage
#' overfitDR(ssDR,y,RR=NULL,test.type=c("phantom","anagenesis","double"),
#' ph=.2,nsim=10,clus=0.5,foldername)
#' 
#' description 
#' The function overfitDR tests the robustness of search.shiftDR results to the effect of incomplete sampling and/or anagenesis.
#' ssDR: the output of 'search.shiftDR' to be tested.
#' y: either a single vector variable or a multivariate dataset as supplied to 'search.shiftDR'. In any case, 'y' must be a named vector.
#' RR: the output of 'RRphylo' to be supplied if 'search.shiftDR' is performed by using 'RRphylo' rates as phenotypic evolutionary rates.
#' test.type: whether to test the effect of incomplete-sampling by adding 'phantom' tips (test.type = "phantom"), the effect of anagenesis by removing species (test.type = "anagenesis"), or both (test.type = "both"). See details for further explanations.
#' ph: the proportion of species to be added or removed or both relative to the tree size. It is 20% of the tree size by default.
#' nsim: the number of simulations to be performed. It is set at 10 by default.
#' clus: the proportion of clusters to be used in parallel computing. To run the single-threaded version of 'overfitDR' set 'clus' = 0.
#' foldername: the path of the folder where plots are to be stored.

#' details 
#' The function returns a list object including:
#' $shift results of clade-associated diversification by accounting for the effect of sampling. For each simluation, results for shifting clades are returned. If the 'ssDR' includes significant results ($single.clades object), the dataframe also includes a column indicating the proportion of species within the 'ssDR' clade which are shared with the clades located by 'overfitDR' at each simulation.
#' $trend results of character-associated diversification by accounting for the effect of sampling. For each simluation, the slope and p-value (p.real) of DR versus phenotype regression, and the p-value derived through randomization (p.random) are registered.


overfitDR<-function(ssDR,
                    y,
                    RR=NULL,
                    test.type=c("phantom","anagenesis","double"),
                    ph=.2,
                    nsim=10,
                    clus=0.5,
                    foldername){
  
  require(phytools)
  require(geiger)
  require(picante)
  require(Rphylopars)
  require(RRphylo)
  require(doParallel)
  require(parallel)
  require(scales)
  require(graphics)
  require(RColorBrewer)
  require(plotrix)
  
  ssDR$tree->tree
  if(ph==0) stop("setting ph=0 is the same as performing no test at all, revert to search.shiftDR or use higher ph value")
  if(length(test.type)>1) stop("Please indicate the argument test.type. It must be either one of 'phantom', 'anagenesis', or 'double'")
  if(!is.binary(tree)) stop("Tree includes polytomies, it is impossible to test for overfit. Consider resolving by function fix.poly")
  if(is.null(foldername)) stop()
  
  if (is.null(nrow(y))) y <- treedata(tree, y, sort = TRUE)[[2]][, 1] else y <- treedata(tree, y, sort = TRUE)[[2]]
  if(!is.null(ssDR$shift.res)&!is.null(ssDR$trend.res)) "both"->st.arg else{
    if(!is.null(ssDR$shift.res)) "shift"->st.arg else "trend"->st.arg
  }
  
  if(!is.null(ssDR$shift.res)){
    if(!is.null(ssDR$shift.res$focal.nodes)) {
      as.numeric(rownames(ssDR$shift.res$focal.nodes))->node
      10^nchar(strsplit(as.character(ssDR$shift.res$focal.nodes[1,2]), "\\.")[[1]][2])->nrep
      f=NULL
    }else{
      node<-NULL
      10^nchar(strsplit(format(ssDR$shift.res$all.clades[1,2],nsmall=3), "\\.")[[1]][2])->nrep
      min(sapply(as.numeric(rownames(ssDR$shift.res$all.clades)),function(x) length(tips(tree,x))))->f
    }
  }
  
  switch(test.type, phantom = {
    res<-list()
    if(round((detectCores() * clus), 0)==0) cl<-makeCluster(1) else cl <- makeCluster(round((detectCores() * clus), 0))
    registerDoParallel(cl)
    res<-foreach(z = 1:nsim, .packages = c("ape","geiger","phytools","doParallel","Rphylopars","RRphylo"),
                 .export="search.shiftDR") %dopar% {
                   
      tree$tip.label[sample(1:Ntip(tree),Ntip(tree)*ph)]->samtips
      tree->tree.ph

      for(i in 1:length(samtips)){
        tree.ph$edge[,2]->eds
        eds[which(eds<=Ntip(tree.ph))]<-tree.ph$tip.label[eds[which(eds<=Ntip(tree.ph))]]
        tree.ph$edge.length[which(eds==samtips[i])]/2->pos
        bind.tip(tree.ph,paste("phantom",i,sep=""),edge.length = pos,
                 position=pos,where=which(tree.ph$tip.label==samtips[i]))->tree.ph
        
      }
      
      fastBM(tree.ph)->xx
      if(min(tree.ph$edge.length)==0) tree.ph$edge.length+1e-03->tree.ph$edge.length    
      
      if(length(y)>Ntip(tree)){
        y[match(names(xx),rownames(y)),]->yy
        names(xx)[which(is.na(match(names(xx),rownames(y))))]->rownames(yy)[which(is.na(match(names(xx),rownames(y))))]
        data.frame(species=rownames(yy),yy,xx)->yy
        phylopars(yy,tree.ph)->ress
        ress$anc_recon[1:Ntip(tree.ph),1:ncol(y)]->y.ph
      }else{
        y[match(names(xx),names(y))]->yy
        names(xx)[which(is.na(match(names(xx),names(y))))]->names(yy)[which(is.na(match(names(xx),names(y))))]
        data.frame(species=names(yy),yy,xx)->yy
        phylopars(yy,tree.ph)->ress
        ress$anc_recon[1:Ntip(tree.ph),1]->y.ph
      }
      if(st.arg!="trend"){
        if(is.null(f)) f.ph<-NULL else{
          if(f<Ntip(tree.ph)/10) f.ph<-f else f.ph<-NULL
        }
        
        if(is.null(node)) node.ph<-NULL  else 
          sapply(as.numeric(node),function(x) getMRCA(tree.ph,tips(tree,x)[which(tips(tree,x)%in%tree.ph$tip.label)]))->node.ph
      }else{
        f.ph<-NULL
        node.ph<-NULL
      }
      
      if(!is.null(RR)){
        RRphylo(tree.ph,y.ph,clus=0)->RR.ph
        search.shiftDR(RR=RR.ph,y=y.ph,select.test = st.arg,node=node.ph,f=f.ph,DR.log=ssDR$DR.log,foldername = tempdir())->ssdr
      }else search.shiftDR(tree=tree.ph,y=y.ph,select.test = st.arg,node=node.ph,f=f.ph,DR.log=ssDR$DR.log,foldername = tempdir())->ssdr

      list(tree.ph,ssdr)
    }
    stopCluster(cl)
    
    lapply(res,"[[",1)->tree.list
    lapply(res,"[[",2)->res.par
    
  }, anagenesis = {
    
    makeL(tree)->L
    apply(L[,1:Nnode(tree)],2,function(x) length(which(x!=0)))->R
    names(which(R==3))->RRs
    as.numeric(RRs)->RRs
    
    tree$edge[tree$edge[,2]<(Ntip(tree)+1),1]->papa
    tree$edge[tree$edge[,2]%in%papa,1]->papa1
    papa1[which(papa1%in%papa)]->papa2
    papa2[-which(papa2%in%RRs)]->papa3
    
    c(papa3,RRs)->papaRR
    
    node.cond<-array()
    tip.cond<-array()
    for (i in 1:length(papaRR)){
      if (tree$edge.length[tree$edge[,1]== papaRR[i]][which(tree$edge[tree$edge[,1]== papaRR[i],2]<(Ntip(tree)+1))]<
          tree$edge.length[tree$edge[,1]== papaRR[i]][which(tree$edge[tree$edge[,1]== papaRR[i],2]>(Ntip(tree)+1))]){ 
        papaRR[i]->node.cond[i]
        tree$edge[tree$edge[,1]== papaRR[i],2][which(tree$edge[tree$edge[,1]== papaRR[i],2]<(Ntip(tree)+1))]->tip.cond[i]
      }
    }
    
    
    if(length(which(is.na(node.cond)))>0) node.cond[which(!is.na(node.cond))]->node.cond
    if(length(which(is.na(tip.cond)))>0) tip.cond[which(!is.na(tip.cond))]->tip.cond
    names(tip.cond)<-tree$tip.label[c(tip.cond)]

    res<-list()
    if(round((detectCores() * clus), 0)==0) cl<-makeCluster(1) else cl <- makeCluster(round((detectCores() * clus), 0))
    registerDoParallel(cl)
    res<-foreach(z = 1:nsim, .packages = c("ape","geiger","phytools","doParallel","RRphylo"),
                 .export="search.shiftDR") %dopar% {
                   
      sample(tip.cond,length(tip.cond)*ph)->tip.cond1
      drop.tip(tree,tip.cond1)->tree.an
      if(length(y)>Ntip(tree)) y[which(rownames(y)%in%tree.an$tip.label),]->y.an else y[which(names(y)%in%tree.an$tip.label)]->y.an

      if(st.arg!="trend"){
        if(is.null(f)) f.an<-NULL else{
          if(f<Ntip(tree.an)/10) f.an<-f else f.an<-NULL
        }
        
        if(is.null(node)) node.an<-NULL  else 
          sapply(as.numeric(node),function(x) getMRCA(tree.an,tips(tree,x)[which(tips(tree,x)%in%tree.an$tip.label)]))->node.an
      }else{
        f.an<-NULL
        node.an<-NULL
      } 
      
      if(!is.null(RR)){
        RRphylo(tree.an,y.an,clus=0)->RR.an
        search.shiftDR(RR=RR.an,y=y.an,select.test = st.arg,node=node.an,f=f.an,DR.log=ssDR$DR.log,foldername = tempdir())->ssdr
      }else search.shiftDR(tree=tree.an,y=y.an,select.test = st.arg,node=node.an,f=f.an,DR.log=ssDR$DR.log,foldername = tempdir())->ssdr
      
      list(tree.an,ssdr)
    }
    stopCluster(cl)
    
    lapply(res,"[[",1)->tree.list
    lapply(res,"[[",2)->res.par
    
  }, double = {
    
    res<-list()
    if(round((detectCores() * clus), 0)==0) cl<-makeCluster(1) else cl <- makeCluster(round((detectCores() * clus), 0))
    registerDoParallel(cl)
    res<-foreach(z = 1:nsim, .packages = c("ape","geiger","phytools","doParallel","Rphylopars","RRphylo"),
                 .export="search.shiftDR") %dopar% {
      
                   tree->tree.ph
      tree$tip.label[sample(1:Ntip(tree),Ntip(tree)*ph)]->samtips
      
      for(i in 1:length(samtips)){
        tree.ph$edge[,2]->eds
        eds[which(eds<=Ntip(tree.ph))]<-tree.ph$tip.label[eds[which(eds<=Ntip(tree.ph))]]
        tree.ph$edge.length[which(eds==samtips[i])]/2->pos
        bind.tip(tree.ph,paste("phantom",i,sep=""),edge.length = pos,
                 position=pos,where=which(tree.ph$tip.label==samtips[i]))->tree.ph
        
      }
      
      fastBM(tree.ph)->xx
      if(min(tree.ph$edge.length)==0) tree.ph$edge.length+1e-03->tree.ph$edge.length    
      
      if(length(y)>Ntip(tree)){
        y[match(names(xx),rownames(y)),]->yy
        names(xx)[which(is.na(match(names(xx),rownames(y))))]->rownames(yy)[which(is.na(match(names(xx),rownames(y))))]
        data.frame(species=rownames(yy),yy,xx)->yy
        phylopars(yy,tree.ph)->ress
        ress$anc_recon[1:Ntip(tree.ph),1:ncol(y)]->y.ph
      }else{
        y[match(names(xx),names(y))]->yy
        names(xx)[which(is.na(match(names(xx),names(y))))]->names(yy)[which(is.na(match(names(xx),names(y))))]
        data.frame(species=names(yy),yy,xx)->yy
        phylopars(yy,tree.ph)->ress
        ress$anc_recon[1:Ntip(tree.ph),1]->y.ph
      }
      
      
      makeL(tree.ph)->L
      apply(L[,1:Nnode(tree.ph)],2,function(x) length(which(x!=0)))->R
      names(which(R==3))->RRs
      as.numeric(RRs)->RRs
      
      tree.ph$edge[tree.ph$edge[,2]<(Ntip(tree.ph)+1),1]->papa
      tree.ph$edge[tree.ph$edge[,2]%in%papa,1]->papa1
      papa1[which(papa1%in%papa)]->papa2
      papa2[-which(papa2%in%RRs)]->papa3
      
      c(papa3,RRs)->papaRR
      
      node.cond<-array()
      tip.cond<-array()
      for (i in 1:length(papaRR)){
        if (tree.ph$edge.length[tree.ph$edge[,1]== papaRR[i]][which(tree.ph$edge[tree.ph$edge[,1]== papaRR[i],2]<(Ntip(tree.ph)+1))]<tree.ph$edge.length[tree.ph$edge[,1]== papaRR[i]][which(tree.ph$edge[tree.ph$edge[,1]== papaRR[i],2]>(Ntip(tree.ph)+1))]){ 
          papaRR[i]->node.cond[i]
          tree.ph$edge[tree.ph$edge[,1]== papaRR[i],2][which(tree.ph$edge[tree.ph$edge[,1]== papaRR[i],2]<(Ntip(tree.ph)+1))]->tip.cond[i]
        }
      }
      
      
      if(length(which(is.na(node.cond)))>0) node.cond[which(node.cond!="NA")]->node.cond
      if(length(which(is.na(tip.cond)))>0) tip.cond[which(tip.cond!="NA")]->tip.cond
      names(tip.cond)<-tree.ph$tip.label[c(tip.cond)]
      #if(length(grep("phantom",names(tip.cond)))>0) tip.cond[-grep("phantom",names(tip.cond))]->tip.cond
      
      
      sample(tip.cond,length(tip.cond)*ph)->tip.cond
      drop.tip(tree.ph,tip.cond)->tree.an
      if(length(y)>Ntip(tree)) y.ph[which(rownames(y.ph)%in%tree.an$tip.label),]->y.an else y.ph[which(names(y.ph)%in%tree.an$tip.label)]->y.an
      
      if(st.arg!="trend"){
        if(is.null(f)) f.an<-NULL else{
          if(f<Ntip(tree.an)/10) f.an<-f else f.an<-NULL
        }
        
        if(is.null(node)) node.an<-NULL  else 
          sapply(as.numeric(node),function(x) getMRCA(tree.an,tips(tree,x)[which(tips(tree,x)%in%tree.an$tip.label)]))->node.an
      }else{
        f.an<-NULL
        node.an<-NULL
      }
      
      if(!is.null(RR)){
        RRphylo(tree.an,y.an,clus=0)->RR.an
        search.shiftDR(RR=RR.an,y=y.an,select.test = st.arg,node=node.an,f=f.an,DR.log = ssDR$DR.log,foldername = tempdir())->ssdr
      }else search.shiftDR(tree=tree.an,y=y.an,select.test = st.arg,node=node.an,f=f.an,DR.log = ssDR$DR.log,foldername = tempdir())->ssdr
      list(tree.an,ssdr)
    }
    stopCluster(cl)
    
    lapply(res,"[[",1)->tree.list
    lapply(res,"[[",2)->res.par
    
  })
  
  if(st.arg=="both"|st.arg=="shift"){
    if(!is.null(node)){
      do.call(rbind,lapply(1:length(res.par), function(x) data.frame(nsim=x,real.node=rownames(ssDR$shift.res$focal.nodes),res.par[[x]]$shift.res$focal.nodes, stringsAsFactors = FALSE)))->shift.df  
      rownames(shift.df)<-NULL
      if(length(unique(shift.df$real.node))>1) lapply(split(shift.df, shift.df$real.node), function(x) x[,-2])->shift else shift.df->shift
      shift.df[which(shift.df$p.DR<=0.025&shift.df$p.sigma<=0.05&shift.df$sigma.ratio>3|shift.df$p.DR>=0.975&shift.df$p.sigma<=0.05&shift.df$sigma.ratio>3),]->shift.df
      
      #### pdf ####
      pdf(file=paste(foldername, "Uncertainty on clade associated diversification plot.pdf",sep="/"))
      if(Ntip(tree)>100) plot(tree,edge.col="gray50", show.tip.label = FALSE) else plot(tree,edge.col="gray50", cex=.8)
      
      brewer.pal(8,"Paired")[c(2,4,6,8)]->colo
      
      plot.info<-get("last_plot.phylo",envir =ape::.PlotPhyloEnv)
      data.frame(plot.info$edge, xx=plot.info$xx[match(plot.info$edge[,2],seq(1:(Ntip(tree)+Nnode(tree))))],
                 yy=plot.info$yy[match(plot.info$edge[,2],seq(1:(Ntip(tree)+Nnode(tree))))])->plot.df
      
      tree$tip.label[plot.df[,2][which(plot.df[,2]<=Ntip(tree))]]->plot.df[,2][which(plot.df[,2]<=Ntip(tree))]
      
      leg<-legend("topleft",legend=c("real shift in DR- and BM rate-","real shift in DR- and BM rate+",
                                     "real shift in DR+ and BM rate-","real shift in DR+ and BM rate+",
                                     "no match with real shift"),pch=16,
                  col=c(colo[c(1,2,4,3)], "gray30"),bty="n",pt.cex=1.7,x.intersp=0.7,cex=.8,plot=FALSE)
      -0.8*leg$rect$h/(max(plot.info$yy)-min(plot.info$yy))->ii
      legend("topleft",legend=c("real shift in DR- and BM rate-","real shift in DR- and BM rate+",
                                "real shift in DR+ and BM rate-","real shift in DR+ and BM rate+",
                                "no match with real shift"),pch=16,
             col=c(colo[c(1,2,4,3)], "gray30"),bty="n",pt.cex=1.7,x.intersp=0.7,cex=.8,inset=c(0,ii),xpd=TRUE)
      
      if(nrow(shift.df)!=0){
        unique(shift.df$real.node)->plot.node
        
        for(s in 1:length(plot.node)){
          plot.node[s]->pn
          
          ssDR$shift.res$focal.nodes[match(pn,rownames(ssDR$shift.res$focal.nodes)),]->ssdr.tips
          shift.df[which(shift.df[,2]==pn),]->shift.pn
          
          
          if(length(which(ssdr.tips$p.DR>0.975&ssdr.tips$sigma.node>ssdr.tips$sigma.back))>0){
            if(length(which(shift.pn$p.DR>0.975&shift.pn$sigma.node>shift.pn$sigma.back))>0){
              shift.pn[which(shift.pn$p.DR>0.975&shift.pn$sigma.node>shift.pn$sigma.back),]->shift.pn
              colo.rect<-colo[3]
            }else colo.rect<-"gray30"
          } else if(length(which(ssdr.tips$p.DR>0.975&ssdr.tips$sigma.node<ssdr.tips$sigma.back))>0){
            if(length(which(shift.pn$p.DR>0.975&shift.pn$sigma.node<shift.pn$sigma.back))>0){
              shift.pn[which(shift.pn$p.DR>0.975&shift.pn$sigma.node<shift.pn$sigma.back),]->shift.pn
              colo.rect<-colo[4]
            }else colo.rect<-"gray30"
          } else if(length(which(ssdr.tips$p.DR<0.025&ssdr.tips$sigma.node>ssdr.tips$sigma.back))>0){
            if(length(which(shift.pn$p.DR<0.025&shift.pn$sigma.node>shift.pn$sigma.back))>0){
              shift.pn[which(shift.pn$p.DR<0.025&shift.pn$sigma.node>shift.pn$sigma.back),]->shift.pn
              colo.rect<-colo[2]
            }else colo.rect<-"gray30"
          } else if(length(which(ssdr.tips$p.DR<0.025&ssdr.tips$sigma.node<ssdr.tips$sigma.back))>0){
            if(length(which(shift.pn$p.DR<0.025&shift.pn$sigma.node<shift.pn$sigma.back))>0){
              shift.pn[which(shift.pn$p.DR<0.025&shift.pn$sigma.node<shift.pn$sigma.back),]->shift.pn
              colo.rect<-colo[1]
            }else colo.rect<-"gray30"
          } else colo.rect<-NULL
          
          if(nrow(shift.pn)!=0&!is.null(colo.rect)) { 
            
            tips(tree, as.numeric(pn))[c(1, length(tips(tree, as.numeric(pn))))]->plot.tip
            plot.df[which(plot.df[,2]%in%plot.tip),4]->coord.tip
            
            sort(coord.tip)->coord.tip
            
            plot.df[which(plot.df[,2]==pn),3]->xsx
            plot.df[which(plot.df[,2]%in%names(which.max(diag(vcv(tree)[tips(tree, as.numeric(pn)),tips(tree, as.numeric(pn))])))),3]->xdx
            
            rect(xsx-0.01,coord.tip[1]-0.01,xdx,coord.tip[2]+0.01,
                 col=alpha(colo.rect, 0.1), border=NA, lty=NULL, lwd=par("lwd"), xpd=FALSE)
            
            nrow(shift.pn)/nsim*100->labe
            text(xsx, coord.tip[2],labels = paste(labe, "%", sep = ""), adj=c(0,1))
            
          }
        }
      }
      dev.off()
      #### End of pdf ####
    }else{
      lapply(res.par, function(x) x$shift.res$single.clades)->shift
      rownames(ssDR$shift.res$single.clades)->ssdr.node
      if(!is.null(ssdr.node)) as.numeric(ssdr.node)->ssdr.node
      
      for(i in 1:length(shift)){
        if(!is.null(shift[[i]])) {
          lapply(as.numeric(rownames(shift[[i]])), function(x) tips(tree.list[[i]],x))->shtips
          sapply(shtips, function(x)  getMRCA(tree,x[which(x%in%tree$tip.label)]))->real.node
          data.frame(nsim=i,node=rownames(shift[[i]]),real.node,shift[[i]])->shift[[i]]
        }

      }
      
      do.call(rbind,shift)->shift.df
      rownames(shift.df)<-NULL
      
      if(!is.null(shift.df)&!is.null(ssdr.node)){
        sapply(shift.df[,3],function(x){
          sapply(ssdr.node,function(w){
            if(length(which(tips(tree,x)%in%tips(tree,as.numeric(w))))>0) length(tips(tree,x))/length(tips(tree,as.numeric(w))) else 0
          })
        })->prop.node
        
        if(!is.null(ncol(prop.node))) {
          t(prop.node)->prop.node
          ncol(prop.node)->dprop
        }else dprop<-1
        data.frame(shift.df,prop.node)->shift.df
        colnames(shift.df)[(ncol(shift.df)-dprop+1):ncol(shift.df)]<-paste("ssdr",ssdr.node,sep="")
      }
      #### pdf ####
      pdf(file=paste(foldername, "Uncertainty on clade associated diversification plot.pdf",sep="/"))
      if(Ntip(tree)>100) plot(tree,edge.col="gray50", show.tip.label = FALSE) else plot(tree,edge.col="gray50", cex=.8)
      brewer.pal(8,"Paired")[c(2,4,6,8)]->colo
      
      plot.info<-get("last_plot.phylo",envir =ape::.PlotPhyloEnv)
      data.frame(plot.info$edge, xx=plot.info$xx[match(plot.info$edge[,2],seq(1:(Ntip(tree)+Nnode(tree))))],
                 yy=plot.info$yy[match(plot.info$edge[,2],seq(1:(Ntip(tree)+Nnode(tree))))])->plot.df
      
      tree$tip.label[plot.df[,2][which(plot.df[,2]<=Ntip(tree))]]->plot.df[,2][which(plot.df[,2]<=Ntip(tree))]
      
      leg<-legend("topleft",legend=c("real shift in DR- and BM rate-","real shift in DR- and BM rate+",
                                     "real shift in DR+ and BM rate-","real shift in DR+ and BM rate+",
                                     "no match with real shift"),pch=16,
                  col=c(colo[c(1,2,4,3)], "gray30"),bty="n",pt.cex=1.7,x.intersp=0.7,cex=.8,plot=FALSE)
      -0.8*leg$rect$h/(max(plot.info$yy)-min(plot.info$yy))->ii
      legend("topleft",legend=c("real shift in DR- and BM rate-","real shift in DR- and BM rate+",
                                "real shift in DR+ and BM rate-","real shift in DR+ and BM rate+",
                                "no match with real shift"),pch=16,
             col=c(colo[c(1,2,4,3)], "gray30"),bty="n",pt.cex=1.7,x.intersp=0.7,cex=.8,inset=c(0,ii),xpd=TRUE)
      
      if(!is.null(shift.df)){  
        unique(shift.df$real.node)->plot.node
        pn.list<-list()
        i=1
        while(length(plot.node)>0){
          plot.node[i]->man.temp
          c(pn.list,list(plot.node[which(plot.node%in%c(getMommy(tree,man.temp), 
                                                        getDescendants(tree,man.temp), man.temp))]))->pn.list
          plot.node[-which(plot.node%in%c(getMommy(tree,man.temp), 
                                          getDescendants(tree,man.temp), man.temp))]->plot.node
        }
        
        for(s in 1:length(pn.list)){
          pn.list[[s]]->pn
          if(!is.null(ssdr.node)){
            shift.df[which(shift.df[,3]==min(pn)),][(ncol(shift.df)-dprop+1):ncol(shift.df)]->leng.tips
            if(ncol(shift.df)>(ncol(shift.df)-dprop+1)) apply(shift.df[which(shift.df[,3]==min(pn)),(ncol(shift.df)-dprop+1):ncol(shift.df)],2,unique)->leng.tips else unique(leng.tips)->leng.tips
            #if(!is.null(ssdr.node)) sapply(ssdr.node, function(x) length(which(tips(tree,min((pn)))%in%tips(tree,x))))->leng.tips else leng.tips<-0
          }else leng.tips<-0
          
          if(any(leng.tips>0)){
            ssDR$shift.res$single.clades[which(rownames(ssDR$shift.res$single.clades)==gsub("ssdr","",names(leng.tips))),]->ssdr.tips
            shift.df[which(shift.df[,3]%in%pn),]->shift.pn
            
            if(is.null(RR)){
              if(length(which(ssdr.tips$p.DR>0.975&ssdr.tips$sigma.node>ssdr.tips$sigma.back))>0){
                if(length(which(shift.pn$p.DR>0.975&shift.pn$sigma.node>shift.pn$sigma.back))>0){
                  shift.pn[which(shift.pn$p.DR>0.975&shift.pn$sigma.node>shift.pn$sigma.back),3]->pn
                  colo.rect<-colo[3] 
                }else colo.rect<-"gray30"
              }
              
              if(length(which(ssdr.tips$p.DR>0.975&ssdr.tips$sigma.node<ssdr.tips$sigma.back))>0){
                if(length(which(shift.pn$p.DR>0.975&shift.pn$sigma.node<shift.pn$sigma.back))>0){
                  shift.pn[which(shift.pn$p.DR>0.975&shift.pn$sigma.node<shift.pn$sigma.back),3]->pn
                  colo.rect<-colo[4] 
                }else colo.rect<-"gray30"
              }
              
              if(length(which(ssdr.tips$p.DR<0.025&ssdr.tips$sigma.node>ssdr.tips$sigma.back))>0){
                if(length(which(shift.pn$p.DR<0.025&shift.pn$sigma.node>shift.pn$sigma.back))>0){
                  shift.pn[which(shift.pn$p.DR<0.025&shift.pn$sigma.node>shift.pn$sigma.back),3]->pn
                  colo.rect<-colo[2] 
                }else colo.rect<-"gray30"
              }
              
              if(length(which(ssdr.tips$p.DR<0.025&ssdr.tips$sigma.node<ssdr.tips$sigma.back))>0){
                if(length(which(shift.pn$p.DR<0.025&shift.pn$sigma.node<shift.pn$sigma.back))>0){
                  shift.pn[which(shift.pn$p.DR<0.025&shift.pn$sigma.node<shift.pn$sigma.back),3]->pn
                  colo.rect<-colo[1] 
                }else colo.rect<-"gray30"
              }
              
            }else{
              
              if(length(which(ssdr.tips$p.DR>0.975&ssdr.tips$p.phen>0.975))>0){
                if(length(which(shift.pn$p.DR>0.975&shift.pn$p.phen>0.975))>0){
                  shift.pn[which(shift.pn$p.DR>0.975&shift.pn$p.phen>0.975),3]->pn
                  colo.rect<-colo[3] 
                }else colo.rect<-"gray30"
              }
              
              if(length(which(ssdr.tips$p.DR>0.975&ssdr.tips$p.phen<0.025))>0){
                if(length(which(shift.pn$p.DR>0.975&shift.pn$p.phen<0.025))>0){
                  shift.pn[which(shift.pn$p.DR>0.975&shift.pn$p.phen<0.025),3]->pn
                  colo.rect<-colo[4] 
                }else colo.rect<-"gray30"
              }
              
              if(length(which(ssdr.tips$p.DR<0.025&ssdr.tips$p.phen>0.975))>0){
                if(length(which(shift.pn$p.DR<0.025&shift.pn$p.phen>0.975))>0){
                  shift.pn[which(shift.pn$p.DR<0.025&shift.pn$p.phen>0.975),3]->pn
                  colo.rect<-colo[2] 
                }else colo.rect<-"gray30"
              }
              
              if(length(which(ssdr.tips$p.DR<0.025&ssdr.tips$p.phen<0.025))>0){
                if(length(which(shift.pn$p.DR<0.025&shift.pn$p.phen<0.025))>0){
                  shift.pn[which(shift.pn$p.DR<0.025&shift.pn$p.phen<0.025),3]->pn
                  colo.rect<-colo[1] 
                }else colo.rect<-"gray30"
              }
            }
            
            
            
          }else colo.rect<-"gray30"
          
          sort(unique(pn), decreasing = TRUE)->pn
          
          for(k in 1:length(pn)){
            
            tips(tree, as.numeric(pn[k]))[c(1, length(tips(tree, as.numeric(pn[k]))))]->plot.tip
            plot.df[which(plot.df[,2]%in%plot.tip),4]->coord.tip
            
            sort(coord.tip)->coord.tip
            
            plot.df[which(plot.df[,2]==pn[k]),3]->xsx
            plot.df[which(plot.df[,2]%in%names(which.max(diag(vcv(tree)[tips(tree, as.numeric(pn[k])),tips(tree, as.numeric(pn[k]))])))),3]->xdx
            
            rect(xsx-0.01,coord.tip[1]-0.01,xdx,coord.tip[2]+0.01,
                 col=alpha(colo.rect, 0.1), border=NA, lty=NULL, lwd=par("lwd"), xpd=FALSE)
            
            length(which(sapply(as.numeric(shift.df$real.node), 
                                function(x) length(which(tips(tree,x)%in%tips(tree,as.numeric(pn[k]))))/length(tips(tree,as.numeric(pn[k]))))==1))/nsim*100->labe
            
            if(k%%2==0) text(xsx, coord.tip[1],labels = paste(labe, "%", sep = ""), adj=c(0,0)) else text(xsx, coord.tip[2],labels = paste(labe, "%", sep = ""), adj=c(0,1))
            
            
          } 
        }
        
        
        shift.df[,-2]->shift.df
        colnames(shift.df)[2]<-"node"
      } 
      dev.off()
      
      shift.df->shift
      #### End of pdf ####
      
    #   if(!is.null(prop.node)&length(prop.node)!=0){
    #     do.call(cbind,prop.node)->prop.node
    #     colnames(prop.node)<-which(!sapply(shift,is.null))
    #     rownames(prop.node)<-rownames(ssDR$shift.res$single.clades)
    #     list(shift=shift.df,node.proportion=prop.node)->shift
    #   }else shift.df->shift
    }
    
  } else shift<-NULL
  
  
  if(st.arg=="both"|st.arg=="trend"){
    # do.call(rbind,lapply(res.par, function(x) x$trend.res))->trend
    # rownames(trend)<-seq(1,nsim)
    
    evol.distinct(tree)->ES
    DR<-1/ES[,2]
    names(DR)<-ES[,1]
    if(isTRUE(ssDR$DR.log)) log(DR)->DR
    
    
    # RRphylo:::range01(y)->yR
    # coef(summary(lm(yR~DR)))[2]->corr
    # predict(lm(yR~DR))[which.min(DR)]->miny
    # range(trend[,1])->ran
    # data.frame(min=(miny-ran[1]*min(DR))+ran[1]*DR,max=(miny-ran[2]*min(DR))+ran[2]*DR,x=DR)->CI
    # CI[order(CI$x),]->CI
    
    pdf(file=paste(foldername, "Uncertainty on character associated diversification plot.pdf",sep="/"))
    if(length(y)>Ntip(tree)){
      if(ncol(y)<7) par(mfrow=c(ceiling(ncol(y)/2),2),mar = c(3.5, 3.5, 1.5, 1)) else par(mfrow=c(3,2),mar = c(3.5, 3.5, 1.5, 1))
      if(!is.null(colnames(y))) ynams<-colnames(y) else ynams<-paste("y",seq(1:ncol(y)),sep="")
    }
    
    trend<-list()
    for(w in 1:ncol(as.matrix(y))){
      if(length(y)>Ntip(tree)) {
        do.call(rbind,lapply(res.par, function(x) x$trend.res[w,]))->trendy 
        colnames(trendy)<-colnames(lapply(res.par, function(x) x$trend.res)[[1]])
      }else{ 
        do.call(rbind,lapply(res.par, function(x) x$trend.res))->trendy
        #colnames(trendy)<-colnames(lapply(res.par, function(x) x$trend.res))
      }
      rownames(trendy)<-seq(1,nsim)
      
      trendy->trend[[w]]
      
      RRphylo:::range01(as.matrix(y)[,w])->yR
      DR[match(names(yR),names(DR))]->DR
      coef(summary(lm(yR~DR)))[2]->corr
      predict(lm(yR~DR))[which.min(DR)]->miny
      range(trendy[,1])->ran
      data.frame(min=(miny-ran[1]*min(DR))+ran[1]*DR,max=(miny-ran[2]*min(DR))+ran[2]*DR,x=DR)->CI
      CI[order(CI$x),]->CI
      
      if(isTRUE(ssDR$DR.log)) xxlab="log real DR" else xxlab="real DR"
      plot(DR,yR,mgp = c(2, 0.5, 0),xlab="",ylab="rescaled real phenotype",cex.lab=1.2)
      if(length(y)>Ntip(tree)) title(main=ynams[w],xlab=xxlab,mgp=c(1.5,0.5,0),cex.lab=1.2) else title(xlab=xxlab,mgp=c(1.5,0.5,0),cex.lab=1.2)
      polygon(c(CI[,3],rev(CI[,3])),c(CI[,1],rev(CI[,2])),border=NA,col=rgb(0.5,0.5,0.5,0.4))
      
      if(length(which(trendy[,3]<=0.025|trendy[,3]>=0.975))>0){
        range(trendy[which(trendy[,3]<=0.025|trendy[,3]>=0.975),1])->ran.sign
        if(length(unique(ran.sign))==1) 
          if(length(y)>Ntip(tree)) sort(c(unique(ran.sign),ssDR$trend.res[w,1]))->ran.sign else sort(c(unique(ran.sign),ssDR$trend.res[1]))->ran.sign
        data.frame(min=(miny-ran.sign[1]*min(DR))+ran.sign[1]*DR,max=(miny-ran.sign[2]*min(DR))+ran.sign[2]*DR,x=DR)->CI.sign
        CI.sign[order(CI.sign$x),]->CI.sign
        polygon(c(CI.sign[,3],rev(CI.sign[,3])),c(CI.sign[,1],rev(CI.sign[,2])),border=NA,col=alpha("aquamarine",0.7))
      }
      ablineclip(lm(yR~DR),x1=min(DR),x2=max(DR),col="aquamarine3",lwd=3,lend=0)
      
    }
    dev.off()
    
    
    # pdf(file=paste(foldername, "Uncertainty on character associated diversification plot.pdf",sep="/"))
    # if(isTRUE(ssDR$DR.log)) xxlab="log real DR" else xxlab="real DR"
    # plot(DR,yR,mgp = c(2, 0.5, 0),xlab=xxlab,ylab="rescaled real phenotype",cex.lab=1.2)
    # polygon(c(CI[,3],rev(CI[,3])),c(CI[,1],rev(CI[,2])),border=NA,col=rgb(0.5,0.5,0.5,0.4))
    # if(length(which(trend[,3]<=0.025|trend[,3]>=0.975))>0){
    #   range(trend[which(trend[,3]<=0.025|trend[,3]>=0.975),1])->ran.sign
    #   if(length(unique(ran.sign))==1) sort(c(unique(ran.sign),ssDR$trend.res[1,1]))->ran.sign
    #   data.frame(min=(miny-ran.sign[1]*min(DR))+ran.sign[1]*DR,max=(miny-ran.sign[2]*min(DR))+ran.sign[2]*DR,x=DR)->CI.sign
    #   CI.sign[order(CI.sign$x),]->CI.sign
    #   polygon(c(CI.sign[,3],rev(CI.sign[,3])),c(CI.sign[,1],rev(CI.sign[,2])),border=NA,col=alpha("aquamarine",0.7))
    # }
    # ablineclip(lm(yR~DR),x1=min(DR),x2=max(DR),col="aquamarine3",lwd=3,lend=0)
    # dev.off()
    
    if(length(y)>Ntip(tree)){
      names(trend)<-ynams
      #do.call(rbind,trend)->trend
    }else do.call(rbind,trend)->trend
    
  } else trend<-NULL
  
  res.tot<-list(shift=shift,trend=trend)
  
  
  return(res.tot)
  
}
