#################Packages and custom functions

library(vegan)
library(letsR)
library(phangorn)
library(paleotree)
library(iNEXT)
library(geoscale)
library(divvy)
library(chronosphere)

Alroy_Forbes<-function(occs)
{
	output<-matrix(nrow=ncol(occs),ncol=ncol(occs))
	rownames(output)<-colnames(occs)
	colnames(output)<-colnames(occs)
	for(i in 1:nrow(output))
	{
		for (j in 1:ncol(output))
		{
			species.A<-rownames(occs)[which(occs[,rownames(output)[i]]!=0)]
			species.B<-rownames(occs)[which(occs[,colnames(output)[j]]!=0)]

			n<-length(unique(c(species.A,species.B)))

			a<-length(intersect(species.A,species.B))
			b<-length(which(species.A%in%species.B==F))
			c<-length(which(species.B%in%species.A==F))
			
			output[i,j]<-1-((a*(n+sqrt(n)))/((a*(n+sqrt(n)))+(3/2*b*c)))
			
		}
	}
	output<-as.dist(output)
	return(output)
}


define.bioregions<-function(tax.dist,geo.dist,method="average",plot=TRUE)
{
	tax.tree<-hclust(tax.dist,method=method)
	geo.tree<-hclust(geo.dist,method=method)

	if(plot==T)
	{
		par(mfrow=c(1,2))
		plot(tax.tree)
		plot(geo.tree)
	}


	geo.phylo<-as.phylo(geo.tree)
	tax.phylo<-as.phylo(tax.tree)
	geo.phylo$edge.length<-geo.phylo$edge.length*2

	tax.subtrees<-subtrees(tax.phylo)
	tax.clusters<-vector(length=length(tax.subtrees),mode="list")
	for(i in 1:length(tax.clusters))
	{
		tax.clusters[[i]]<-sort(tax.subtrees[[i]]$tip.label)
	}





	node.heights<-sort(dateNodes(geo.phylo),decreasing=F)
	test.node.heights<-node.heights[which(names(node.heights)%in%as.character(1:length(geo.phylo$tip.label))==F)]

	for(i in 1:geo.phylo$Nnode)
	{

		if(i==1)
		{
			level<-1		
			height.ranges<-vector(length=1,mode="list")
			height.ranges[[1]]<-vector(length=2)
			height.ranges[[1]][1]<-0
	
			bioregions<-vector(length=1,mode="list")
			bioregions[[1]]<-1:length(geo.phylo$tip.label)
			names(bioregions[[1]])<-sort(geo.phylo$tip.label)

		}
		else
		{
			clust.node<-as.numeric(names(test.node.heights)[i-1])
			children<-sort(geo.phylo$tip.label[Descendants(geo.phylo, clust.node, type ="tips")[[1]]])
	
			for(j in 1:length(tax.clusters))
			{
				if(identical(tax.clusters[[j]],children)==T)
				{
					level<-level+1
					
					new.bioregions<-vector(length=level,mode="list")
					new.bioregions[1:length(bioregions)]<-bioregions
					new.bioregions[[level]]<-new.bioregions[[level-1]]				
	
					change<-children[which(new.bioregions[[level]][children]!=new.bioregions[[level]][children[1]])]
					
					for(k in 1:length(change))
					{
						new.bioregions[[level]][which(new.bioregions[[level]]>new.bioregions[[level]][change][k])]<-new.bioregions[[level]][which(new.bioregions[[level]]>new.bioregions[[level]][change][k])]-1
					}
					new.bioregions[[level]][change]<-new.bioregions[[level]][children[1]]				
					bioregions<-new.bioregions
	
					new.height.ranges<-vector(length=level,mode="list")
					new.height.ranges[1:length(height.ranges)]<-height.ranges
					new.height.ranges[[level]]<-vector(length=2)
					new.height.ranges[[level-1]][2]<-test.node.heights[i-1]
					new.height.ranges[[level]][1]<-test.node.heights[i-1]
					height.ranges<-new.height.ranges

		
				}
			}
		}
	}
	height.ranges[[length(height.ranges)]][2]<-max(test.node.heights)
		
	results<-list(bioregions,height.ranges)
	names(results)<-c("bioregions","height.ranges")	
	return(results)
}	



################	
################
################ Data import and prep
################
################



data<-read.csv(file.choose())#####Read in either supplementary Data 3, 4 or 5


collections<-unique(data$collection_no)#######identify unique collections
ages<-read.csv(file.choose(),row.names=1)#######Read in supplementary data 7


#######Age range of collections
collection.ages<-matrix(nrow=length(collections),ncol=2)


for(i in 1:length(collections))
{
	rows<-data[which(data$collection_no==collections[i]),]
	collection.bins<-rows[1,c("early_interval","late_interval")]

	if(collection.bins[1,2]=="")
	{
		collection.ages[i,1]<-ages[collection.bins[1,1],1]
		collection.ages[i,2]<-ages[collection.bins[1,1],2]

	}
	else
	{
		collection.ages[i,1]<-ages[collection.bins[1,1],1]
		collection.ages[i,2]<-ages[collection.bins[1,2],2]
	}
}
rownames(collection.ages)<-collections

time.bins<-read.csv(file.choose(),row.names=1)######Read in supplementary Data 8
plot.time<-rowMeans(time.bins)




############
############
############ Identify bioregions in each time bin
############
############


bin.bioregs<-vector(length=nrow(time.bins),mode="list")
spatial.scale<-100

for(i in 1:nrow(time.bins))
{
	bin.collect<-collections[intersect(which(collection.ages[,1]>time.bins[i,2]),
		which(collection.ages[,2]<time.bins[i,1]))]
	bin.rows<-data[which(data$collection_no%in%bin.collect),]	

	bin.coords<-matrix(nrow=length(bin.collect),ncol=2)
	rownames(bin.coords)<-bin.collect


	if(is.null(nrow(bin.coords))==F && nrow(bin.coords)>2)
	{
		

		higher.tax<-bin.rows[,c("class","order","family",
			"genus","accepted_name")]

		higher.tax$order[which(higher.tax$order=="NO_ORDER_SPECIFIED")]<-""
		higher.tax$family[which(higher.tax$family=="NO_FAMILY_SPECIFIED")]<-""
		higher.tax$accepted_name[which(bin.rows$accepted_rank!="species")]<-""

		higher.tax<-higher.tax[order(higher.tax$accepted_name,higher.tax$genus,
			higher.tax$family,higher.tax$order,higher.tax$class,decreasing=T),]



		if(nrow(higher.tax)>1)
		{

			tax.incl<-higher.tax$accepted_name[
				which(higher.tax$accepted_name!="")]
			next.level<-higher.tax[which(higher.tax$accepted_name!=""),
				"genus"]
			current.tax<-higher.tax[which(higher.tax$accepted_name==""),]

			for(k in 4:1)
			{
				if(nrow(current.tax)!=0)
				{
					tax.to.test<-current.tax[which(current.tax[,k]!=""),k]
					tax.to.add<-which(tax.to.test
						%in%next.level==F)
					tax.incl<-c(tax.incl,tax.to.test[tax.to.add])
					if(k!=1)
					{
						next.level<-current.tax[tax.incl,(k-1)]
					}
					current.tax<-current.tax[which(current.tax[,k]==""),]
				}
			}
			collection.occs<-tax.incl
		}
		else(collection.occs<-bin.rows$accepted_name)


	}




	bin.spec<-unique(collection.occs)


	bin.pres.abs<-matrix(nrow=length(bin.spec),ncol=nrow(bin.coords),data=0)
	rownames(bin.pres.abs)<-bin.spec
	colnames(bin.pres.abs)<-bin.collect

	for(j in 1:length(bin.collect))
	{
		collect.row<-which(bin.rows$collection_no==bin.collect[j])
		bin.coords[j,1]<-bin.rows[collect.row[1],"paleolng"]
		bin.coords[j,2]<-bin.rows[collect.row[1],"paleolat"]

		collect.spec<-which(bin.spec%in%bin.rows[collect.row,"accepted_name"])
		bin.pres.abs[collect.spec,j]<-1
	}

	bin.pres.abs<-bin.pres.abs[,which(colSums(bin.pres.abs)>0)]
	if(is.null(ncol(bin.pres.abs))==F )
	{
		bin.coords<-bin.coords[colnames(bin.pres.abs),]
		bin.coords<-bin.coords[is.na(bin.coords[,1])==F,]
		if(is.null(nrow(bin.coords))==F && nrow(bin.coords)>2)
		{


			bin.pres.abs<-bin.pres.abs[,rownames(bin.coords)]

			geo.dist<-lets.distmat(bin.coords)
			tax.dist<-Alroy_Forbes(bin.pres.abs)

			bioregions<-define.bioregions(tax.dist,geo.dist)
			for(j in 1:length(bioregions$height.ranges))
			{
				if(bioregions$height.ranges[[j]][1]<spatial.scale && 
					bioregions$height.ranges[[j]][2]>spatial.scale)
				{
					bin.bioregs[[i]]<-bioregions$bioregions[[j]]
				}
			}
		}
	}	
}



#######################
####################### Find diversity of all synapsids, stem mammals and 
####################### crown mammals in each bioregion
#######################


bioreg.div<-vector(length=length(bin.bioregs),mode="list")
ther.div<-vector(length=length(bin.bioregs),mode="list")
mam.div<-vector(length=length(bin.bioregs),mode="list")


ther.count<-vector(length=length(bin.bioregs),mode="list")
mam.count<-vector(length=length(bin.bioregs),mode="list")

for(i in 1:length(bin.bioregs))
{
	if(length(bin.bioregs[[i]])>0)
	{

		no.bioregs<-max(bin.bioregs[[i]])
	
		for(j in 1:no.bioregs)
		{
			bio.collect<-as.numeric(names(bin.bioregs[[i]])[which(bin.bioregs[[i]]==j)])
	
			rows<-which(data$collection_no%in%bio.collect)
			collect.data<-data[rows,]

	
			higher.tax<-collect.data[,c("class","order","family",
			"genus","accepted_name")]

			higher.tax$order[which(higher.tax$order=="NO_ORDER_SPECIFIED")]<-""
			higher.tax$family[which(higher.tax$family=="NO_FAMILY_SPECIFIED")]<-""
			higher.tax$accepted_name[which(collect.data$accepted_rank!="species")]<-""
		
			higher.tax<-higher.tax[order(higher.tax$accepted_name,higher.tax$genus,
				higher.tax$family,higher.tax$order,higher.tax$class,decreasing=T),]



			if(nrow(higher.tax)>1)
			{

				tax.incl<-higher.tax$accepted_name[
					which(higher.tax$accepted_name!="")]
				next.level<-higher.tax[which(higher.tax$accepted_name!=""),
					"genus"]
				current.tax<-higher.tax[which(higher.tax$accepted_name==""),]
	
				for(k in 4:1)
				{
					if(nrow(current.tax)!=0)
					{
						tax.to.test<-current.tax[which(current.tax[,k]!=""),k]
						tax.to.add<-which(tax.to.test
							%in%next.level==F)
						tax.incl<-c(tax.incl,tax.to.test[tax.to.add])
						if(k!=1)
						{
							next.level<-current.tax[tax.incl,(k-1)]
						}
						current.tax<-current.tax[which(current.tax[,k]==""),]
					}
				}
				collection.occs<-tax.incl
			}
			else(collection.occs<-collect.data$accepted_name)

			occ.counts<-data.frame(table(collection.occs))
				
			groups<-vector(length=length(occ.counts))
			for(k in 1:nrow(occ.counts))
			{
				groups[k]<-collect.data[which(collect.data$accepted_name==
					occ.counts$collection.occs[k])[1],"class"]
			}



			if("Osteichthyes"%in%groups)
			{
				ther.occs<-occ.counts[which(groups=="Osteichthyes"),]
				ther.count[[i]]<-c(ther.count[[i]],nrow(ther.occs))
			}
			if("Mammalia"%in%groups)
			{
				mam.occs<-occ.counts[which(groups=="Mammalia"),]
				mam.count[[i]]<-c(mam.count[[i]],nrow(mam.occs))
			}

					
			if(length(collection.occs)>1 && length(unique(collection.occs))>1 )
			{
				bioreg.div[[i]]<-c(bioreg.div[[i]],
					estimateD(occ.counts$Freq,base="coverage"
					,level=0.9)[1,6])

				

				if("Osteichthyes"%in%groups)
				{

					if(ncol(ther.occs)>1 &&
					 	length(unique(ther.occs$collection.occs))>1)
					{


						ther.div[[i]]<-c(ther.div[[i]],
							estimateD(ther.occs$Freq,
							base="coverage",level=0.9)[1,6])


					}

				}
				if("Mammalia"%in%groups)
				{

					if(ncol(mam.occs)>1 &&
					 	length(unique(mam.occs$collection.occs))>1)
					{


						mam.div[[i]]<-c(mam.div[[i]],
							estimateD(mam.occs$Freq,
							base="coverage",level=0.9)[1,6])


					}

				}



			}

		}
		
	}
	print(i)
}	



div<-vector(length=length(bioreg.div))
for(i in 1:length(bioreg.div))
{
	if(length(bioreg.div[[i]])>0)
	{
		div[i]<-mean(bioreg.div[[i]],na.rm=T)
	}
}




mam.med.div<-vector(length=length(mam.div))
for(i in 1:length(mam.div))
{
	if(length(mam.div[[i]])>0)
	{
		mam.med.div[i]<-10^mean(log10(mam.div[[i]]),na.rm=T)
	}
}

mam.med.div[which(mam.med.div==0)]<-NA

ther.med.div<-vector(length=length(ther.div))
for(i in 1:length(ther.div))
{
	if(length(ther.div[[i]])>0)
	{
		ther.med.div[i]<-10^mean(log10(ther.div[[i]]),na.rm=T)
	}
}

ther.med.div[which(ther.med.div==0)]<-NA





############
############ Relative richness
############


bioreg.prop<-vector(length=length(bin.bioregs),mode="list")


for(i in 1:length(bin.bioregs))
{
	if(length(bin.bioregs[[i]])>0)
	{

		no.bioregs<-max(bin.bioregs[[i]])
	
		for(j in 1:no.bioregs)
		{
			bio.collect<-as.numeric(names(bin.bioregs[[i]])[which(bin.bioregs[[i]]==j)])
	
			rows<-which(data$collection_no%in%bio.collect)
			collect.data<-data[rows,]

	
			higher.tax<-collect.data[,c("class","order","family",
			"genus","accepted_name")]

			higher.tax$order[which(higher.tax$order=="NO_ORDER_SPECIFIED")]<-""
			higher.tax$family[which(higher.tax$family=="NO_FAMILY_SPECIFIED")]<-""
			higher.tax$accepted_name[which(collect.data$accepted_rank!="species")]<-""
		
			higher.tax<-higher.tax[order(higher.tax$accepted_name,higher.tax$genus,
				higher.tax$family,higher.tax$order,higher.tax$class,decreasing=T),]



			if(nrow(higher.tax)>1)
			{

				tax.incl<-higher.tax$accepted_name[
					which(higher.tax$accepted_name!="")]
				next.level<-higher.tax[which(higher.tax$accepted_name!=""),
					"genus"]
				current.tax<-higher.tax[which(higher.tax$accepted_name==""),]
	
				for(k in 4:1)
				{
					if(nrow(current.tax)!=0)
					{
						tax.to.test<-current.tax[which(current.tax[,k]!=""),k]
						tax.to.add<-which(tax.to.test
							%in%next.level==F)
						tax.incl<-c(tax.incl,tax.to.test[tax.to.add])
						if(k!=1)
						{
							next.level<-current.tax[tax.incl,(k-1)]
						}
						current.tax<-current.tax[which(current.tax[,k]==""),]
					}
				}
				collection.occs<-tax.incl
			}
			else(collection.occs<-collect.data$accepted_name)

			occ.counts<-data.frame(table(collection.occs))
				
			groups<-vector(length=length(occ.counts))
			for(k in 1:nrow(occ.counts))
			{
				groups[k]<-collect.data[which(collect.data$accepted_name==
					occ.counts$collection.occs[k])[1],"class"]
			}


			if("Mammalia" %in% groups && "Osteichthyes" %in% groups)
			{
				prop<-length(which(groups=="Mammalia"))/length(groups)
				bioreg.prop[[i]]<-c(bioreg.prop[[i]],prop)
			}
		}
	}
}


prop.med<-lapply(bioreg.prop,mean,na.rm=T)
	



########
######## Geographic range
########


syn.mst<-vector(length=nrow(time.bins))
mam.mst<-vector(length=nrow(time.bins))
ther.mst<-vector(length=nrow(time.bins))


for(i in 1:nrow(time.bins))
{
	bin.collect<-collections[intersect(which(collection.ages[,1]>time.bins[i,2]),
		which(collection.ages[,2]<time.bins[i,1]))]
	bin.rows<-data[which(data$collection_no%in%bin.collect),]	


	bin.coords<-matrix(nrow=length(bin.collect),ncol=2)
	rownames(bin.coords)<-bin.collect

	for(j in 1:length(bin.collect))
	{
		collect.row<-which(bin.rows$collection_no==bin.collect[j])
		bin.coords[j,1]<-bin.rows[collect.row[1],"paleolng"]
		bin.coords[j,2]<-bin.rows[collect.row[1],"paleolat"]
	}

	bin.coords<-bin.coords[which(is.na(bin.coords[,1])==F),]

	if(is.null(nrow(bin.coords))==F)
	{
		geo.dist<-as.matrix(lets.distmat(bin.coords))

		tree<-mst(geo.dist)
		lengths<-tree
		for(j in 1:nrow(tree))
		{
			for(k in j:ncol(tree))
			{
				if(k!=j)
				{
					if(tree[j,k]==1)
					{
						lengths[j,k]<-lengths[k,j]<-geo.dist[j,k]
					}
				}
			}
		}
		syn.mst[i]<-sum(as.dist(lengths))

		groups<-bin.rows$class

		if("Mammalia"%in%groups)
		{
			mammal.rows<-bin.rows[which(groups=="Mammalia"),]
			mammal.collect<-as.character(unique(
				mammal.rows$collection_no))
			mammal.collect<-mammal.collect[which(
				mammal.collect%in%rownames(geo.dist))]

			mammal.dist<-geo.dist[mammal.collect,mammal.collect]
			if(is.null(nrow(mammal.dist))==F)
			{


				tree<-mst(mammal.dist)
				lengths<-tree
				for(j in 1:nrow(tree))
				{
					for(k in j:ncol(tree))
					{
						if(k!=j)
						{
							if(tree[j,k]==1)
							{
								lengths[j,k]<-lengths[k,j]<-
									mammal.dist[j,k]
							}
						}
					}
				}
				mam.mst[i]<-sum(as.dist(lengths))
			}
		}
		if("Osteichthyes"%in%groups)
		{
			ther.rows<-bin.rows[which(groups=="Osteichthyes"),]
			ther.collect<-as.character(unique(
				ther.rows$collection_no))
			ther.collect<-ther.collect[which(
				ther.collect%in%rownames(geo.dist))]

			ther.dist<-geo.dist[ther.collect,ther.collect]

			if(is.null(nrow(ther.dist))==F)
			{

				tree<-mst(ther.dist)
				lengths<-tree
				for(j in 1:nrow(tree))
				{
					for(k in j:ncol(tree))
					{
						if(k!=j)
						{
							if(tree[j,k]==1)
							{
								lengths[j,k]<-lengths[k,j]<-
									ther.dist[j,k]
							}
						}
					}
				}
				ther.mst[i]<-sum(as.dist(lengths))
			}
		}
	}	
}

syn.mst[syn.mst==0]<-NA
mam.mst[mam.mst==0]<-NA
ther.mst[ther.mst==0]<-NA