##############################################
##############################################
########
######## Arguments:
######## tree - object of class phylo; the tree over which evolution is 
########        simulated 
########
######## max.teeth - class numeric; the maximum number of teeth allowed
########
######## param.H.1 - class numeric; instantaneous transition rate of level 1 of 
########           the developmental controls: toothlessness or not
########
######## param.H.2 - class numeric; instantaneous transition rate of level 2 of 
########           the developmental controls: presence or abence of teeth in
########           each bone (Maxilla or premaxilla)
########
######## param.H.3 - class numeric; instantaneous transition rate of level 3 of 
########           the developmental controls: presence or abence of teeth in
########           each dental region (incisor, caniniform, premolar or 
########           molar)
########
######## param.T - class numeric; instantaneous transition rate of level 4 of 
########           the developmental controls: addition or subtraction of 
########           teeth one by one
########
######## param.R - class numeric; sigma2 (i.e. rate) parameter of the brownian
########           motion models used to simulate changes in relative size of 
########           the regions (incisor, caniniform, premolar or molar)
########
##############################################
##############################################
 






tooth.sim<-function(tree,max.teeth=40,param.H.1=0.002,param.H.2=0.002,param.H.3=0.002,param.T=0.5,param.R=0.04)
{



	###################################
	###################################
	#######                      ######
	####### q matrix for level 4 ######
	#######                      ######
	###################################
	###################################




	q.4<-matrix(nrow=max.teeth,ncol=max.teeth,data=0)
	
	for(i in 1:ncol(q.4))
	{
		if(i!=1)
		{
			q.4[i-1,i]<-param.T
		}
		if(i!=ncol(q.4))
		{
			q.4[i+1,i]<-param.T
		}
		q.4[i,i]<- -sum(q.4[,i])
	}
	rownames(q.4)<-colnames(q.4)<-1:max.teeth


	###########################################
	###########################################
	#######                              ######
	####### simulate tooth row evolution ######
	#######                              ######
	###########################################
	###########################################




	tooth.row<-sim.Mk(tree, q.4, anc=NULL, nsim=1)
	tooth.row<-as.numeric(as.vector(tooth.row))
	names(tooth.row)<-tree$tip.label

	##############################################
	##############################################
	#######                                 ######
	####### Total tooth development control ######
	#######                                 ######
	##############################################
	##############################################

	if(is.null(param.H.1)==F)
	{
		q.1<-matrix(nrow=2,ncol=2,data=c(-param.H.1,param.H.1,param.H.1,-param.H.1))
		rownames(q.1)<-colnames(q.1)<-1:2

		tooth.1<-sim.Mk(tree, q.1, anc=1, nsim=1)
	}
	else
	{
		tooth.1<-rep("1",length(tree$tip.label))
		names(tooth.1)<-tree$tip.label
	}

	#####################################################
	#####################################################
	#######                                        ######
	####### tooth development control in each bone ######
	#######                                        ######
	#####################################################
	#####################################################


	if(is.null(param.H.2)==F)
	{
		q.2<-matrix(nrow=2,ncol=2,data=c(-param.H.2,param.H.2,param.H.2,-param.H.2))
		rownames(q.2)<-colnames(q.2)<-1:2
		premax.2<-sim.Mk(tree, q.2, anc="1", nsim=1)
		max.2<-sim.Mk(tree, q.2, anc="1", nsim=1)
	}
	else
	{
		premax.2<-rep("1",length(tree$tip.label))
		max.2<-rep("1",length(tree$tip.label))
		names(premax.2)<-names(max.2)<-tree$tip.label
	}
	
	#############################################################
	#############################################################
	#######                                                ######
	####### tooth development control in maxillary regions ######
	#######                                                ######
	#############################################################
	#############################################################



	if(is.null(param.H.3)==F)
	{
		q.3<-matrix(nrow=2,ncol=2,data=c(-param.H.3,param.H.3,param.H.3,-param.H.3))
		rownames(q.3)<-colnames(q.3)<-1:2
		precan.3<-sim.Mk(tree, q.3, anc="1", nsim=1)
		can.3<-sim.Mk(tree, q.3, anc="1", nsim=1)
		postcan.3<-sim.Mk(tree, q.3, anc="1", nsim=1)
	}
	else
	{
		precan.3<-rep("1",length(tree$tip.label))
		can.3<-rep("1",length(tree$tip.label))
		postcan.3<-rep("1",length(tree$tip.label))
		names(precan.3)<-names(can.3)<-names(postcan.3)<-tree$tip.label
	}



	####################################
	####################################
	#######                       ######
	####### Relative region sizes ######
	#######                       ######
	####################################
	####################################


	if(is.null(param.R)==F)
	{
		region.bar<-matrix(nrow=length(tree$tip.label),ncol=3)
		region.bar[,1]<-fastBM(tree,a=0.25,sig2=param.R,bounds=c(0,1))
		region.bar[,2]<-fastBM(tree,a=0.5,sig2=param.R,bounds=c(0,1))
		region.bar[,3]<-fastBM(tree,a=0.75,sig2=param.R,bounds=c(0,1))
		
		region.props<-matrix(nrow=length(tree$tip.label),ncol=4)


		for(i in 1:nrow(region.props))
		{
			bar<-c(0,sort(region.bar[i,]),1)
			
			for(j in 2:length(bar))
			{
				region.props[i,j-1]<-bar[j]-bar[j-1]
			}
		}
		region.no<-round(tooth.row*region.props)
		colnames(region.no)<-c("Incisor","Caniniform","Premolar","Molar")
		rownames(region.no)<-tree$tip.label
	}
	else(region.no<-tooth.row)

	###################################
	###################################
	#######                      ######
	####### Suppress development ######
	#######                      ######
	###################################
	###################################

	if(is.matrix(region.no))
	{
		for(j in 1:nrow(region.no))
		{
			if(tooth.1[j]=="2")
			{
				region.no[j,]<-0
			}
			else 
			{
				if(premax.2[j]=="2")
				{
					region.no[j,1]<-0
				}
				if(max.2[j]=="2")
				{
					region.no[j,2:4]<-0
				}
				else
				{
					if(precan.3[j] =="2")
					{
						region.no[j,2]<-0
					}
					if(can.3[j] =="2")		
					{
						region.no[j,3]<-0
					}
					if(postcan.3[j] =="2")		
					{
						region.no[j,4]<-0
					}
				}
			}
		}
	}
	return(region.no)
}