#' FSO
#'
#' usage
#' FSO(tree,y,ssdr.args=list(select.test=NULL,f=NULL),
#' odr.args=list(test.type=NULL,nsim=NULL,clus=NULL),BMrates=TRUE,RRrates=TRUE,
#' foldername,iter)
#' 
#' description 
#' tree: a polytomic phylogenetic tree. Needs not to be ultrametric.
#' y: either a single vector variable or a multivariate dataset. In any case, 'y' must be a named vector.
#' ssdr.args: a list of arguments to be passed to 'search.shiftDR', including 'list(select.test=NULL,f=NULL)'. If unspecified, default values will be used.
#' odr.args: a list of arguments to be passed to 'overfitDR', including 'list(test.type=NULL,nsim=NULL,clus=NULL)'. If unspecified, the function runs under 'test.type' = "double", with 'nsim' = 1 and 'clus' = 0.5.
#' BMrates: if 'TRUE' the function performs the analyses by computing phenotypic evolutionary rates by using package geomorph's 'compare.evol.rates'.
#' RRrates: if 'TRUE' the function performs the analyses by computing phenotypic evolutionary rates by using 'RRphylo'.
#' foldername: the path of the folder where figures and results are to be stored.
#' iter: the number of iterations.
#' 
#' 
#' details 
#' The function stores a "FSO.rda" file into the 'foldername'. The file contains:
#' tree.list: the list of phylogenetic trees dichotomized by 'fix.poly'.
#' SSdr.list: and/or 'SSdrRR.list' the list of 'search.shiftDR' outputs as derived by setting 'BMrates'=TRUE and 'RRrates'=TRUE, respectively.
#' Odr.list: and/or 'OdrRR.list' the list of 'overfitDR' outputs as derived by setting 'BMrates'=TRUE and 'RRrates'=TRUE, respectively.
#' RR.list: the list of 'RRphylo' outputs.
#' The function also returns a summary output for clade-associated diversification testing, under both 'BMrates'=TRUE and 'RRrates'=TRUE. Refer to 'search.shiftDR' output for further details about the tables.


FSO<-function(tree,y,ssdr.args=list(select.test=NULL,f=NULL),odr.args=list(test.type=NULL,nsim=NULL,clus=NULL),
              BMrates=TRUE,RRrates=TRUE,foldername,iter){
  require(geiger)
  require(phytools)
  require(ape)
  require(RColorBrewer)
  require(RRphylo)
  require(scales)
  
  ssdr.args$select.test->seltest
  ssdr.args$f->f
  if(is.null(seltest)) seltest<-"both"
  
  odr.args$clus->clus
  odr.args$nsim->ns
  odr.args$test.type->testy
  if(is.null(clus)) clus<-0.5
  if(is.null(testy)) testy<-"double"
  if(is.null(ns)) ns<-1
  
  tree.list<-SSdr.list<-Odr.list<-SSdrRR.list<-OdrRR.list<-RR.list<-list()
  for(i in 1:iter){
    print(i)
    
    fix.poly(tree,"resolve")->tree.f->tree.list[[i]]
    if (is.null(nrow(y))) yf <- treedata(tree, y, sort = TRUE)[[2]][, 1] else yf <- treedata(tree, y, sort = TRUE)[[2]]
    
    if(BMrates){
      search.shiftDR(tree=tree.f,y=yf,select.test = seltest,f=f,foldername=tempdir())->SSdrf->SSdr.list[[i]]
      overfitDR(ssDR=SSdrf,y=yf,test.type = testy,nsim=ns,clus=clus,foldername=tempdir())->Odr.list[[i]]
    }
    
    if(RRrates){
      RRphylo(tree=tree.f,y=yf)->RRf->RR.list[[i]]
      search.shiftDR(RR=RRf,y=yf,select.test = seltest,f=f,foldername=tempdir())->SSdrRRf->SSdrRR.list[[i]]
      overfitDR(ssDR=SSdrRRf,y=yf,RR=RRf,test.type = testy,nsim=ns,clus=clus,foldername=tempdir())->OdrRR.list[[i]]
    }
  }
  
  if(BMrates&RRrates) save(tree.list=tree.list,SSdr.list=SSdr.list,Odr.list=Odr.list,RR.list=RR.list,SSdrRR.list=SSdrRR.list,OdrRR.list=OdrRR.list,file=paste(foldername,"FSO results.rda",sep="/"))
  if(!RRrates) save(tree.list=tree.list,SSdr.list=SSdr.list,Odr.list=Odr.list,file=paste(foldername,"FSO results.rda",sep="/"))
  if(!BMrates) save(tree.list=tree.list,RR.list=RR.list,SSdrRR.list=SSdrRR.list,OdrRR.list=OdrRR.list,file=paste(foldername,"FSO results.rda",sep="/"))
  
  if(BMrates){
    if(any(!sapply(lapply(Odr.list,"[[",1),is.null))){    
      Odr.list->odr.list
      plot.node<-shift.list<-list()
      for(i in 1:length(which(!sapply(lapply(odr.list,"[[",1),is.null)))){
        which(!sapply(lapply(odr.list,"[[",1),is.null))[i]->h
        if(nrow(lapply(odr.list,"[[",1)[[h]])>1){
          sapply(lapply(odr.list,"[[",1)[[h]]$node,function(w) getMRCA(tree,tips(tree.list[[h]],as.numeric(as.character(w)))))->plot.node[[i]]
          data.frame(iter=h,nsim=lapply(odr.list,"[[",1)[[h]][,1],original=plot.node[[i]],lapply(odr.list,"[[",1)[[h]][,-1])->shift.list[[i]]
        }else{
          getMRCA(tree,tips(tree.list[[h]],as.numeric(as.character(lapply(odr.list,"[[",1)[[h]]$node))))->plot.node[[i]]
          data.frame(iter=h,nsim=lapply(odr.list,"[[",1)[[h]][,1],original=plot.node[[i]],lapply(odr.list,"[[",1)[[h]][,-1])->shift.list[[i]]
        }
      }
      unlist(plot.node)->plot.node
      do.call(rbind,lapply(shift.list,function(x) 
        if(length(grep("ssdr",colnames(x)))>0) x[,-grep("ssdr",colnames(x))] else x))->shift.df
      
      tiff(paste(getwd(),"treeFSO BMrates.tiff",sep="/"),
           width=3524,
           height=4286,
           pointsize=8,
           res=600,
           compression="lzw")
      
      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+"),pch=16,
                  col=colo[c(1,2,4,3)],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+"),pch=16,
             col=colo[c(1,2,4,3)],bty="n",pt.cex=1.7,x.intersp=0.7,cex=.8,inset=c(0,ii),xpd=TRUE)
      
      
      
      unique(plot.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
      }
      
      percentages<-list()
      for(s in 1:length(pn.list)){
        pn.list[[s]]->pn
        shift.df[which(shift.df[,3]%in%pn),]->shift.pn
        
        apply(shift.pn,1,function(w){
          if(length(which(w[6]>0.975&w[8]>w[9]))>0) return(colo[3])
          if(length(which(w[6]>0.975&w[8]<w[9]))>0) return(colo[4])
          if(length(which(w[6]<0.025&w[8]>w[9]))>0) return(colo[2])
          if(length(which(w[6]<0.025&w[8]<w[9]))>0) return(colo[1])
        })->colo.rect
        
        sort(unique(pn), decreasing = TRUE)->pn
        perc<-array()
        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[k], 0.1), border=NA, lty=NULL, lwd=par("lwd"), xpd=FALSE)
          
          if (k==1) text(xdx,coord.tip[1]+(coord.tip[2]-coord.tip[1])/2,adj=0,font=2,label=pn[k],xpd=TRUE)
          
          length(which(sapply(as.numeric(shift.df$original), 
                              function(x) length(which(tips(tree,x)%in%tips(tree,as.numeric(pn[k]))))/length(tips(tree,as.numeric(pn[k]))))==1))/(ns*iter)*100->perc[k]
          
        }
        names(perc)<-pn
        perc->percentages[[s]]
      }
      
      paste(round(do.call(c,percentages),3),"%",sep="")->node.perc
      names(node.perc)<-names(do.call(c,percentages))
      dev.off()
      
      shift.df[,-4]->shift.df
      colnames(shift.df)[3]<-"node"
      list(sum.tab=shift.df,node.percent=node.perc)->sum.BMrates.res
    }else NULL->sum.BMrates.res
  }else NULL->sum.BMrates.res
  
  if(RRrates){
    if(any(!sapply(lapply(OdrRR.list,"[[",1),is.null))){ 
      OdrRR.list->odr.list
      plot.node<-shift.list<-list()
      for(i in 1:length(which(!sapply(lapply(odr.list,"[[",1),is.null)))){
        which(!sapply(lapply(odr.list,"[[",1),is.null))[i]->h
        if(nrow(lapply(odr.list,"[[",1)[[h]])>1){
          sapply(lapply(odr.list,"[[",1)[[h]]$node,function(w) getMRCA(tree,tips(tree.list[[h]],as.numeric(as.character(w)))))->plot.node[[i]]
          data.frame(iter=h,nsim=lapply(odr.list,"[[",1)[[h]][,1],original=plot.node[[i]],lapply(odr.list,"[[",1)[[h]][,-1])->shift.list[[i]]
        }else{
          getMRCA(tree,tips(tree.list[[h]],as.numeric(as.character(lapply(odr.list,"[[",1)[[h]]$node))))->plot.node[[i]]
          data.frame(iter=h,nsim=lapply(odr.list,"[[",1)[[h]][,1],original=plot.node[[i]],lapply(odr.list,"[[",1)[[h]][,-1])->shift.list[[i]]
        }
      }
      unlist(plot.node)->plot.node
      do.call(rbind,lapply(shift.list,function(x) 
        if(length(grep("ssdr",colnames(x)))>0) x[,-grep("ssdr",colnames(x))] else x))->shift.dfRR
      
      tiff(paste(foldername,"treeFSO RRrates.tiff",sep="/"),
           width=3524,
           height=4286,
           pointsize=8,
           res=600,
           compression="lzw")
      
      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 RR rate-","real shift in DR- and RR rate+",
                                     "real shift in DR+ and RR rate-","real shift in DR+ and RR rate+"),pch=16,
                  col=colo[c(1,2,4,3)],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 RR rate-","real shift in DR- and RR rate+",
                                "real shift in DR+ and RR rate-","real shift in DR+ and RR rate+"),pch=16,
             col=colo[c(1,2,4,3)],bty="n",pt.cex=1.7,x.intersp=0.7,cex=.8,inset=c(0,ii),xpd=TRUE)
      
      
      
      unique(plot.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
      }
      
      percentagesRR<-list()
      for(s in 1:length(pn.list)){
        pn.list[[s]]->pn
        shift.dfRR[which(shift.dfRR[,3]%in%pn),]->shift.pn
        
        apply(shift.pn,1,function(w){
          if(length(which(w[6]>0.975&w[8]>0.975))>0) return(colo[3])
          if(length(which(w[6]>0.975&w[8]<0.025))>0) return(colo[4])
          if(length(which(w[6]<0.025&w[8]>0.975))>0) return(colo[2])
          if(length(which(w[6]<0.025&w[8]<0.025))>0) return(colo[1])
        })->colo.rect
        
        sort(unique(pn), decreasing = TRUE)->pn
        percRR<-array()
        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[k], 0.1), border=NA, lty=NULL, lwd=par("lwd"), xpd=FALSE)
          
          if (k==1) text(xdx,coord.tip[1]+(coord.tip[2]-coord.tip[1])/2,adj=0,font=2,label=pn[k],xpd=TRUE)
          
          length(which(sapply(as.numeric(shift.dfRR$original), 
                              function(x) length(which(tips(tree,x)%in%tips(tree,as.numeric(pn[k]))))/length(tips(tree,as.numeric(pn[k]))))==1))/(ns*iter)*100->percRR[k]
          
        }
        names(percRR)<-pn
        percRR->percentagesRR[[s]]
      }
      
      paste(round(do.call(c,percentagesRR),3),"%",sep="")->node.percRR
      names(node.percRR)<-names(do.call(c,percentagesRR))
      dev.off()
      
      shift.dfRR[,-4]->shift.dfRR
      colnames(shift.dfRR)[3]<-"node"
      list(sum.tab=shift.dfRR,node.percent=node.percRR)-> sum.RRrates.res
    }else NULL->sum.RRrates.res
  }else NULL->sum.RRrates.res
  
  
  if(!BMrates){ 
    if(any(!sapply(lapply(OdrRR.list,"[[",2),is.null))){
      do.call(rbind,lapply(1:length(OdrRR.list),function(x){
        data.frame(iter=x,nsim=rownames(OdrRR.list[[x]][[2]]),OdrRR.list[[x]][[2]])
      }))->sum.trend
    }else sum.trend<-NULL
  }else{
    if(any(!sapply(lapply(Odr.list,"[[",2),is.null))){
      do.call(rbind,lapply(1:length(Odr.list),function(x){
        data.frame(iter=x,nsim=rownames(Odr.list[[x]][[2]]),Odr.list[[x]][[2]])
      }))->sum.trend 
    }else sum.trend<-NULL
  }
  
  return(list(shift.BMrates=sum.BMrates.res,shift.RRrates=sum.RRrates.res,trend=sum.trend))
}
