#name: eta_over_time

#title: Phylogenetic \eta over time

#description: The peak height, \eta, of the spectral density profile (Lewitus & Morlon, Syst Bio, 2016; Lewitus & Rolland, Virus Evolution, 2019) at subtrees extending from each internal node of a given tree.

#usage: eta_over_time(tr,ntips,threshold,depth,plot.tree=F,log.plot=F,return.tree=F)

#arguments:
#tr=an object of type 'phylo'
#ntips=the minimum number of tips in a subtree to calculate \eta for
#threshold=the top fraction (0.1=10%) of \eta values for subtrees that will be circled in the tree (if plot.tree=T) and in the plot of \eta~nodes (if log.plot=T)
#depth=the depth of the tree to start extracting subtrees from (1=root;2=tips)
#plot.tree=if TRUE, tree will be plotted
#log.plot=if TRUE, log \eta for subtrees extracted for each node will be plotted
#return.tree=if TRUE, extracted subtrees will be returned

#value:
#1. \eta values by subtree
#2. subtrees (if return.tree=T) 

#example:
#tr=rtree(100)
#res<-(tr,5,0.1,1.1,plot.tree=T,log.plot=T,return.tree=F)

library(RPANDA)
library(igraph)

eta_over_time<-function(tr,ntips,threshold,depth,plot.tree=F,log.plot=F,return.tree=F){
			#MGL functions
			integr<-function(x,f){
				#get lengths of var and integrand
       			n=length(x)

       			#trapezoidal integration
       			int=0.5*sum((x[2:n]-x[1:(n-1)])*
       							(f[2:n]+f[1:(n-1)]))

		      return(int)				
				}
       			
			gete<-function(phy){
				eigen(
				graph.laplacian(
				graph.adjacency(
				data.matrix(dist.nodes(phy)),
				weighted=T),
				normalized=F),
				only.values=T)$values	}	
			getStats<-function(e){
				density(e)->d
				d$y/integr(d$x,d$y)->dsc
				max(dsc)->height
				return(height)	}
				
			
			#check for polytomies, negative branch-lengths
			if(min(tr$edge.length)<=0){
				stop("phylogeny must have resolved polytomies")
			}
				else{tr<-tr}
			
			#extract subtrees at each node
			q1<-round(length(tr$tip.label)*depth)+1
			q2<-2*length(tr$tip.label)-2
			subtr<-lapply(c(q1:q2),function(a){
				extract.clade(tr,a)})
			
			#estimate subtree ages	
			ages<-sapply(1:length(subtr),function(m){
				max(node.age(subtr[[m]])$ages)	})
				
			#MGL stats	
			es<-lapply(subtr,gete)
			stats<-sapply(es,getStats)
			
			#create table
			tab<-cbind(stats,sapply(subtr,Ntip),c(c(q1:q2)-1),ages)
			colnames(tab)<-c('height','tips','node','age')
			tab<-as.data.frame(tab)
			
			#plot tree
			if(plot.tree==T){
				tabt<-subset(tab,tab$tips>ntips)
					tabh<-tabt[order(tabt$height,decreasing=T),]
					sel=tabh[1:round(threshold*dim(tabh)[1]),c(1,3)]
						nodes<-sel[,2];heights<-sel[,1]
					col<-colorRampPalette(c('hotpink','chartreuse3','cornflowerblue'))
				par(mfrow=c(2,1),mar=c(2,4,1,1))
				plot(ladderize(tr),show.tip.label=F,edge.width=0.3)
					nodelabels('o',node=nodes,frame='none',
						col=col(length(nodes)))
				
				#plot \eta~nodes
					if(log.plot==T){
						plot(log(tabt$height)~tabt$node,type='l',axes=F)
						points(nodes,log(heights),
						col=col(length(nodes)))}
					else{plot(tabt$height~tabt$node,type='l',axes=F)
						points(nodes,heights,col=col(length(nodes)))}
						axis(1,at=round(seq(min(tabt$node),
						max(tabt$node),
						diff(range(tabt$node)/20))))	
						axis(2,las=2)
							}
			
			#return subtrees
			if(return.tree==T){				
			return(list(tab,subtr))}
			else{return(tab)}
			}