#rooted=FALSE if the root (external environment) has not been added already (in which case it is done internally)
#r00,qq and bb are the three parameters of the eco-FOB model. They can be a single value or a vector of one value per species
#tstep is the step size, currently only used in the g-FOB part of the simulations
#if max_time=NULL, it will run each simulation until the last species goes extinct. If a value max_time>0 is given, extinction times of some species will not be known
#nsims is the number of replicate simulations wanted
#output is a character vector with some combination of:
###"eco_G_meanvar" (run both models, store only mean & variances of richness through time
####"eco_richness_time" and/or "G_richness_time" (store species richness at each time step in each simulation)
####,"eco_species_extimes" and/or "G_species_extimes" (store extinction times of all species in each simulation)

simulate_ecoFOB<-function(web,rooted=FALSE,r00=0.2,qq=0.8,bb=1,tstep=0.01,max_time=NULL,nsims=10,output=c("eco_G_meanvar"),mean_model="conditional",G_sampling=FALSE,print_output=FALSE,print_freq=10){

wsum_calc<-function(x)sum(x*extant)
calc_richness<-function(x)sum(!is.na(x))
calc_meanRT<-function(x)ifelse(sum(!is.na(x))>0,mean(x,na.rm=T),NA)
find_extinction<-function(x)which(x==0)[1]

if(all(output%in%c("eco_G_meanvar","eco_richness_time","G_richness_time","eco_species_extimes","G_species_extimes"))==F)stop("improperly specified output")

sources<-c(1,1+which(rowSums(web)==0))

if(print_output)print("Sources identified as:")
if(print_output)print(as.numeric(sources))
	n<-dim(web)[1]+1
	if(rooted==F){
		web<-cbind(rep(0,n-1),web);web<-rbind(rep(0,n),web)
		web[rowSums(web)==0,1]<-1
			}
	n<-dim(web)[1]
	p<-web/rowSums(web)
	wsum<-apply(p,1,sum)
if(print_output)print("Sum of diet proportions (should all =1)")
if(print_output)print(wsum)	
	
if(print_output)print("Beginning eco-FOB simulation")

rlist<-list();maxtimes<-numeric()
xtimes<-matrix(NA,n,nsims,dimnames=list(paste("s",1:n,sep="")))

for(k in 1:nsims){
if(print_output)if(k%%print_freq==0)print(k)

	if(length(r00)==1){
		r0<-c(NA,rep(r00,n-1))
		q<-c(NA,rep(qq,n-1))
		}else{
			r0<-c(NA,r00);q<-c(NA,qq)
			}
	b<-c(NA,rep(bb,n-1))
			
if(print_output)if(k==1){print("r0:");print(r0)}
if(print_output)if(k==1){print("q:");print(q)}
if(print_output)if(k==1){print("b:");print(b)}
	
r<-r02<-r0
extant<-rep(1,n)
	
ttemp<-tstep
mx<-matrix(r,n,1,dimnames=list(paste("s",1:n,sep=""),ttemp))
while(sum(extant)>1){
	#determine next species to go extinct given current rates (As far as I can tell, this is equivalent to iterating time steps)
	dtimes<-max(ttemp)+tstep+rexp(n,r)#gets a warning, but its ok
	ttemp<-seq(max(ttemp)+tstep,min(dtimes,na.rm=T),by=tstep)
	mx<-cbind(mx,matrix(r,n,length(ttemp),dimnames=list(NULL,ttemp)))
	extant[which.min(dtimes)]<-0;r0[which.min(dtimes)]<-NA
	xtimes[which.min(dtimes),k]<-dtimes[which.min(dtimes)]
	wsum<-apply(p,1,wsum_calc)
	r<-r0+(q-r0)*(1-wsum)^b
}
if(!is.null(max_time))mx<-mx[,which(as.numeric(colnames(mx))<=max_time)]
rlist[[k]]<-mx
maxtimes[k]<-dim(mx)[2];names(maxtimes)[k]<-colnames(mx)[dim(mx)[2]]
}

if(print_output)print("eco-FOB simulations completed (ignore warnings)")

if(print_output)print("computing richness thru time")
maxt<-as.numeric(names(maxtimes))
tt<-seq(tstep,max(maxt),by=tstep)
richness<-list()
for(k in 1:nsims)richness[[k]]<-apply(rlist[[k]],2,calc_richness)
richmx<-matrix(0,nsims,max(maxtimes))
for(i in 1:nsims)richmx[i,1:maxtimes[i]]<-richness[[i]]
richness<-richmx

if("eco_G_meanvar"%in%output|"G_richness_time"%in%output|"G_species_extimes"%in%output){

meanrichness<-apply(richness,2,mean,na.rm=T)
varrichness<-apply(richness,2,var,na.rm=T)
	
if(print_output)print("building species' rate matrix")
spratemx<-list()
spratemx[[1]]<-matrix(NA,nsims,max(maxtimes))
for(i in 2:n){
	spratemx[[i]]<-matrix(NA,nsims,max(maxtimes))
	for(k in 1:nsims){
		spratemx[[i]][k,1:maxtimes[k]]<-rlist[[k]][i,]
		}
	}	

if(mean_model=="unconditional"){
for(i in 2:n)spratemx[[i]][is.na(spratemx[[i]])]<-ifelse(i%in%sources,r02[i],q[i])
	}

if(print_output)print("calculating species' rates thru time")
	
meanRT<-matrix(0,n,max(maxtimes),dimnames=list(paste("s",1:n,sep=""),tt))
for(i in 2:n){
	meanRT[i,]<-apply(spratemx[[i]],2,calc_meanRT)
	meanRT[i,is.na(meanRT[i,])]<-ifelse(i%in%sources,r02[i],q[i])
	}
	
richnessG<-matrix(0,nsims,length(tt),dimnames=list(paste("r",1:nsims,sep=""),tt))
dimnames(richness)<-dimnames(richnessG)
xtimesG<-matrix(NA,n,nsims,dimnames=list(paste("s",1:n,sep=""),paste("r",1:nsims,sep="")))

if(print_output)print("Beginning g-FOB simulation")

for(k in 1:nsims){
if(print_output)if(k%%print_freq==0)print(k)
	
extant<-rep(1,n)
for(i in 1:length(tt)){
if(sum(extant)>1){
	r<-c(0,meanRT[,i][-1])*extant
	extinctions<-rbinom(n,1,r*tstep)
	if(sum(extinctions)>0)xtimesG[which(extinctions==1),k]<-i
	extant<-extant-extinctions
	richnessG[k,i]<-sum(extant)-1
	}
}
}
if(print_output)print("g-FOB simulations completed (ignore warnings)")
	
meanrichnessG<-apply(richnessG,2,mean,na.rm=T)
varrichnessG<-apply(richnessG,2,var,na.rm=T)

}

out<-list()
if("eco_G_meanvar"%in%output)out[[which(output=="eco_G_meanvar")]]<-data.frame(tt,meanrichness,meanrichnessG,varrichness,varrichnessG)
if("eco_richness_time"%in%output)out[[which(output=="eco_richness_time")]]<-richness
if("G_richness_time"%in%output)out[[which(output=="G_richness_time")]]<-richnessG
if("eco_species_extimes"%in%output)out[[which(output=="eco_species_extimes")]]<-xtimes
if("G_species_extimes"%in%output)out[[which(output=="G_species_extimes")]]<-xtimesG*tstep

if(length(out)==1)out<-out[[1]]

out
}

web<-read.csv("chesapeake.csv",row.names=1)
web[web>0]<-1

out1<-simulate_ecoFOB(web,nsims=5,output=c("eco_G_meanvar"))

out2<-simulate_ecoFOB(web,nsims=5,output=c("eco_richness_time","G_richness_time"),print_output=TRUE)

out3<-simulate_ecoFOB(web,nsims=5,output=c("eco_species_extimes","G_species_extimes"),print_output=TRUE)

out4<-simulate_ecoFOB(web,nsims=5,output=c("eco_G_meanvar","eco_richness_time","G_richness_time","eco_species_extimes","G_species_extimes"),print_output=TRUE)