###############################
# linear stochastic simulator #
###############################

# a simulator for the model
# N(t+1)=r(t)*N(t)
# where r(t) is given by an empirical distribution

# INPUT: 
# r the empirical distribution of r values (choosen with equal likelihood)
# reps the number of simulations to run
# T the length of each simulation
# N0 initial density
# plottype 0 means none, 1 means population trajectories on regular scale, 2 means log densities and histogram


linear=function(r=c(2,0.2),reps=1,T=50,N0=100,plottype=1){
	N=matrix(N0,T,reps)
	for(i in 1:(T-1)){
     	N[i+1,]=sample(r,reps,replace=TRUE)*N[i,]}
    if(plottype==1){par(cex.lab=1.5,cex.axis=1.5); matplot(N,type="l",pch=21,bg="red",lwd=3,xlab="time",ylab="N")}
    if(plottype==2){par(cex.lab=1.5,cex.axis=1.5,mfrow=c(1,2),mar=c(4,4,1,1));
    	matplot(log(N),type="l",xlab="time",ylab="log(density)")
    	abline(a=log(N0),b=mean(log(r)),lwd=4,lty=2)
    	hist(log(N[T,]),40,col="red",yaxt="n",ylab="",xlab="log(density)",main="")
    }
}

########################
## Bethedging Function #
########################

# this function determines the germination fraction g that maximizes
# the expected log fitness log(gY+(1-g)s) where
# Y is a discrete random variable taking on the values Y=c(....)
# with probabilities p=c(....) and s is probability of a seed suviving

# INPUT
# Y list of yield values
# p list of probabilities of yield values
# s seed survival
# plotit - TRUE creates a plot, FALSE supresses the plot

# OUTPUT
# say out=bethedging(....) then
# out$maximum is the g value that maximizes the log-fitness
# out$objective is the maximal log-fitness

hedge=function(Y=c(100,0),p=c(0.5,0.5),s=0.9,plotit=TRUE){
  f=function(g)sum(p*log(g*Y+(1-g)*s));
  temp=optimize(f,interval=c(0,1),maximum=TRUE);
  if(plotit){
    g=seq(0,0.999999,length=100)
    out=sapply(g,f)
    par(cex.lab=1.25,cex.axis=1.25,mfrow=c(1,1),mar=c(4.5,4.5,1,1))
    plot(g,out,type="l",lwd=3,col="red",xlab="g",ylab="log(fitness)")
    abline(h=0,lty=2)
    abline(v=temp$maximum,lty=2,col="blue")
    abline(h=temp$objective,lty=2,col="blue")} 
  return(temp);}
  
  
###################################
## Stochastic Beverton-Holt model #
###################################

# a simulator for the model
# N(t+1)=r(t)*N(t)/(1+a*N[t])
# where r(t) is given by an empirical distribution
# and a is the strength of density dependence

# INPUT: 
# r the empirical distribution of r values (choosen with equal likelihood)
# a strength of intraspecific competition
# reps the number of simulations to run
# T the length of each simulation
# N0 initial density
# plottype 0 means none
#          1 means population trajectories on regular scale
#          2 means log densities and histogram of all trajectories at time T
#          3 means densties and historgram of first trajectrom from T/2 to T



nonlinear=function(r=c(1.2,0.95),reps=1,T=200,N0=1,plottype=1){
  N=matrix(N0,T,reps)
  for(i in 1:(T-1)){
    N[i+1,]=sample(r,reps,replace=TRUE)*N[i,]/(1+0.001*N[i,])}
  if(plottype==1){par(cex.lab=1.5,cex.axis=1.5); matplot(N,type="l",pch=21,bg="red",lwd=3,xlab="time",ylab="N")}
  if(plottype==2){par(cex.lab=1.5,cex.axis=1.5,mfrow=c(1,2),mar=c(4,4,1,1));
                  matplot(log(N),type="l",xlab="time",ylab="log(density)")
                  hist(log(N[T,]),40,col="red",yaxt="n",ylab="",xlab="log(density)",main="")
  }
  if(plottype==3){par(cex.lab=1.5,cex.axis=1.5,mfrow=c(1,2),mar=c(4,4,1,1));
                  matplot(N[,1],type="l",xlab="time",ylab="density")
                  hist(N[(T/2):T,1],40,xlim=c(0,max(N[,1])),col="red",yaxt="n",ylab="",xlab="density",main="")
  }
}