

#######Disparity analysis
trees<-read.nexus(file.choose())#####Read in posterior tree file
sample.trees<-trees[sample(5000:20000,100,replace=F)]####sample 100 random trees

matrix<-ReadMorphNexus(file.choose()) ###Read in matrix
data<-read.csv(file.choose(),row.name=1)#### Read in data file
time<-358:248 #### create vector of time slices


#####Empty matrices to store results. each column is a time slice, each row is an evolutioanry history

carn.grad.disp<-matrix(ncol=length(time), nrow=10000)#####Carnivore disparity
herb.grad.disp<-matrix(ncol=length(time), nrow=10000)#####herbivore disparity
amn.grad.disp<-matrix(ncol=length(time), nrow=10000)######amniote disparity






for(w in 1:100)
{



	#####################################################################
	##Rescale Tree (Posterior tree file is not scaled to absolute time)##
	#####################################################################



	tree<-ladderize(sample.trees[[w]])

	eocae<-199.3####age of eocaecilia
	capt<-289###age of captorhinus magnus
	true.diff<-capt-eocae
	obs.diff<-dateNodes(tree)[which(tree$tip.label=="Captorhinus_magnus")]
	scale<-true.diff/obs.diff

	tree$edge.length<-tree$edge.length*scale
	tree$root.time<-max(dateNodes(tree))+199.3

	###################################################
	##distance matrix, including estimating ancestors##
	###################################################

	if(0%in%tree$edge.length)
	{
		zero.tree<-tree
		zero.tree$edge.length<-zero.tree$edge.length+0.0000001
		ancs<-AncStateEstMatrix(matrix,zero.tree)
	}
		

	else(ancs<-AncStateEstMatrix(matrix,tree))

	
	dist<-MorphDistMatrix(ancs,Distance="MORD")
	dist<-dist$"DistanceMatrix"

	
	tree$node.label<-540:1077#####label nodes to maintain identity of nodes when assigning diet. Change to 514:1025 when using Ford tree


	for(z in 1:100)
	{
		#################################
		##ancestral herb reconstruction##
		#################################
		herb<-data[1:(nrow(data)-4),"Herb"]
		names(herb)<-rownames(data)[1:(nrow(data)-4)]
		herb[which(herb==0)]<-sample(c(1,2),length(which(herb==0)),replace=T)#####randomly assign diet to those for which it is uncertain
	
		sim<-make.simmap(tree,herb,nsim=1)
		sim.herb<-as.numeric(summary(sim)[[4]])




		#####################################################
		##drop outgroups and taxa with no characters scored##
		#####################################################

		seym.amn.node<-getMRCA(tree,c("Leptoropha","Pyozia"))
		root.age<-dateNodes(tree)[seym.amn.node]
		test.tree<-drop.tip(tree,tree$tip.label[which(apply(matrix[[2]][[3]], 1, function(x) all(is.na(x))))])
		test.tree$root.time<-root.age
		geoscalePhylo(test.tree,direction="upwards", boxes="User")

		

		test.herb<-c(herb[test.tree$tip.label],sim.herb[test.tree$node.label-539])####change to 513 when using Ford matrix
		names(test.herb)[(length(test.tree$tip.label)+1):length(test.herb)]<-test.tree$node.label


		###################
		##distance matrix##
		###################

		test.dist<-dist[c(test.tree$tip.label,test.tree$node.label),c(test.tree$tip.label,test.tree$node.label)]


		##################
		##Add tip ranges##
		##################

		tip.ages<-data[test.tree$tip.label,2:5]
		tip.dates<-dateNodes(test.tree)[1:length(test.tree$tip.label)]
		disp.tree<-test.tree
	
		for(i in 1:length(test.tree$tip.label))
		{
			x<-100000
			y<-0
			while(x>tip.dates[i] && y<100)
			{
				x<-runif(1,max=tip.ages[i,3],min=tip.ages[i,4])
				y<-y+1
			}
			if(y!=100)
			{
				range<-tip.dates[i]-x

				edge<-which(disp.tree$edge[,2]==i)
				disp.tree$edge.length[edge]<-disp.tree$edge.length[edge]+range
			}
		}



		####################
		##Find edge ranges##
		####################
	
		dates<-dateNodes(disp.tree)
		edge.ages<-disp.tree$edge
		
		for(i in 1:nrow(edge.ages))
		{
			edge.ages[i,1]<-dates[edge.ages[i,1]]
			edge.ages[i,2]<-dates[edge.ages[i,2]]
		}
		

		##################
		##Find disparity##
		##################
	
		for(i in 1:length(time))
		{
			bin.edge<-intersect(which(edge.ages[,2] <= time[i]),which(edge.ages[,1] >= time[i]))
			if(length(bin.edge)>1)
			{

				#####identify taxa in bin	
				bin.tax<-vector(length=length(bin.edge))
				for(j in 1:length(bin.tax))
				{
					branch.age<-edge.ages[bin.edge[j],]
					prop.along<-(branch.age[1]-time[i])/(branch.age[1]-branch.age[2])
					probs<-c((1-prop.along),prop.along)
					bin.tax[j]<-disp.tree$edge[bin.edge[j],sample(1:2,1,prob=probs)]
				}	
				bin.dist<-dist[bin.tax,bin.tax] ####distances in bin
		
				bin.herb<-which(test.herb[bin.tax]==2)####which are herbivores
				bin.carn<-which(test.herb[bin.tax]==1)####which are carnivores

							

				herb.dist<-bin.dist[bin.herb,bin.herb]
				carn.dist<-bin.dist[bin.carn,bin.carn]
			
				amn.grad.disp[(((w-1)*10)+z),i]<-median(as.dist(bin.dist),na.rm=T)

				carn.grad.disp[(((w-1)*10)+z),i]<-median(as.dist(carn.dist),na.rm=T)
				herb.grad.disp[(((w-1)*10)+z),i]<-median(as.dist(herb.dist),na.rm=T)


			}
		}

	}

}









##########rates######
mcc.tree<-read.beast(file.choose())#####read in MCC tree


####extract rates
rates<-as.numeric(mcc.tree@data$effectivebrlenIgrBrlens_mean)
rates<-rates[order(as.numeric(mcc.tree@data$node))]
rates<-rates[which(rates!=1)]
rates/tree$edge.length




###Empty matrices to store results
amn.rates<-matrix(nrow=100,ncol=length(time))

carn.rates<-matrix(nrow=100,ncol=length(time))
herb.rates<-matrix(nrow=100,ncol=length(time))





for(w in 1:100)
{


	tree$node.label<-540:1077 #####label nodes to maintain identity of nodes when assigning diet. Change to 514:1025 when using Ford tree


	
	#################################
	##ancestral herb reconstruction##
	#################################
	herb<-data[1:(nrow(data)-4),"Herb"]
	names(herb)<-rownames(data)[1:(nrow(data)-4)]
	herb[which(herb==0)]<-sample(c(1,2),length(which(herb==0)),replace=T)
	
	sim<-make.simmap(tree,herb,nsim=1)
	sim.herb<-as.numeric(summary(sim)[[4]])




	####################################################
	##drop outgroups and taa with no characters scored##
	####################################################


	seym.amn.node<-getMRCA(tree,c("Leptoropha","Pyozia"))
	root.age<-dateNodes(tree)[seym.amn.node]
	test.tree<-drop.tip(tree,c("TEMNOSPONDYLI","BAPHETIDAE","ANTHRACOSAURIA"))
	test.tree$root.time<-root.age
	geoscalePhylo(test.tree,direction="upwards", boxes="User")
	

	test.herb<-c(herb[test.tree$tip.label],sim.herb[test.tree$node.label-539])####change to 513 when using Ford matrix  
	names(test.herb)[(length(test.tree$tip.label)+1):length(test.herb)]<-test.tree$node.label

	####################
	##Find edge ranges##
	####################
	
	dates<-dateNodes(test.tree)
	edge.ages<-test.tree$edge
	
	for(i in 1:nrow(edge.ages))
	{
		edge.ages[i,1]<-dates[edge.ages[i,1]]
		edge.ages[i,2]<-dates[edge.ages[i,2]]
	}


	for(i in 1:length(time))
	{
		bin.edge<-intersect(which(edge.ages[,2] <= time[i]),which(edge.ages[,1] >= time[i]))
		if(length(bin.edge)>1)
		{


			bin.tax<-test.tree$edge[bin.edge,2]


			bin.herb<-which(test.herb[bin.tax]==2)
			bin.carn<-which(test.herb[bin.tax]==1)
							
			carn.rates[w,i]<-median(rates[bin.carn])

			herb.rates[w,i]<-median(rates[bin.herb])

			amn.rates[w,i]<-median(rates[bin.herb])
		}

	}
	print(w)
}






#########Saturation test#########


morph.tree<-read.nexus(file.choose())####read in tree where brnch lengths = morphological changes
morph.tree<-drop.tip(morph.tree,"TEMNOSPONDYLI")####drop outgroup




morph.dist<-dist[[2]][morph.tree$tip.label,morph.tree$tip.label]#####morphological distanc
pat.dist<-as.matrix(distTips(morph.tree,method = "patristic"))####patristic distance


affinities<-as.vector(data[morph.tree$tip.label,"Clade"])#####affinities
names(affinities)<-morph.tree$tip.label


herb<-data[1:(nrow(data)-4),"Herb"]#####herbivorous taxa
names(herb)<-rownames(data)[1:(nrow(data)-4)]
test.herb<-herb[morph.tree$tip.label]


######pairwise comparisons
pairs<-matrix(nrow=0,ncol=8)
colnames(pairs)<-c("Taxon 1","Taxon 2","Herb 1","Herb 2","Clade 1","Clade 2","Morph dist","Patrist dist")

for(i in 1:nrow(morph.dist))
{
	for(j in i:ncol(morph.dist))
	{
		tax.1<-rownames(morph.dist)[i]
		tax.2<-colnames(morph.dist)[j]
		herb.1<-test.herb[tax.1]
		herb.2<-test.herb[tax.2]
		clade.1<-affinities[tax.1]
		clade.2<-affinities[tax.2]
		pair.morph<-morph.dist[tax.1,tax.2]
		pair.pat<-pat.dist[tax.1,tax.2]
		
		pairs<-rbind(pairs,c(tax.1,tax.2,herb.1,herb.2,clade.1,clade.2,pair.morph,pair.pat))
	}
}





