#may.R   
#R worksheet to accompany Robert M. May, 1974 ,  Biological populations with non-overlapping generations:  
#  stable points, stable cycles and chaos , Science, 186: 645-647.

#This worksheet has three parts.  
#  In the first two parts  we look at the model itself, and especially at the long term behaviour 
#  of the trajectories and how this may depend on initial conditions, growth parameters, etc. 
#  In section 1 we consider the dynamics of a single species; 
#  in section 2 we consider two competing species.  (Not done yet)
#  The third part introduces the idea of a bifurcation diagram, which is a record of the repeated 
#  "stable points" that occur for different growth parameter values.

#Section 1.  Stable cycles and chaos.

#May examines the behavior of difference equation analogues of the logistic equation.  
#  He looks at two equations:
#       1.  N(t+1) = N(t) * exp[r * (1 - N(t)/K]
#       2.  N(t+1) = N(t) * [ 1 + r * (1 - N(t)/K)]  
#Both are non-linear equations, and both have a stable equilibrium at N = K, when 2 > r > 0. 
# The second, although more clearly "logistic-based" has the unfortunate tendancy to permit overshoot 
# in which the population can dip negative under extreme oscillation, even with reasonably small stepsize. 
# The reader is invited to experiment with this model as well, 
#  but in this worksheet all computations in this model have been "commented out".

#We can demonstrate this stability by calculating the trajectory of populations from an initial size of 10 individuals,
# over 30 time steps.  We will try several values of r from 0.5 to 2.0.  
#We will use equation 1 as May did in his paper.  
#To avoid having extraneous "junk" getting printed, we suppress all printing and plotting in the main loop 
#(under control of r), give the plots names, and plot them afterwards.  
# Two plots are created: plot_list gives N as a function of t; pairs gives N(t+1) against N(t). 
# Clustering of points in the latter is an indication of stable cyclic behaviour.

N<-vector(length=31)
N[1]<- 10 ;  K<-1000 ; # initial condition and `carrying capacity`
 title1<-"N(t+1) vs N(t) when r=" 
 title2<-"N(t) vs t when r= " # title blanks for plots
 sample_r<-c (0.4, 0.7, 1.0, 1.3, 1.6, 1.9);
 P_frames<- NULL; R_frames<- NULL;
 par(mfrow=c(2,3))
 
 for (r in sample_r)
 {	
	
    for (t in 1:30)
		{
		#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] );  # logistic model
		N[t+1]<- N[t] * exp(r * (1 - N[t] / K))   #ricker model
		}
  plot(N[1:(length(N)-1)],N[2:length(N)], col="green", main=paste(title1, r),
			xlab="N[t]", ylab="N[t+1]")
  }

  for (r in sample_r)
 {	
	
    for (t in 1:30)
		{
		#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] );  # logistic model
		N[t+1]<- N[t] * exp(r * (1 - N[t] / K))   #ricker model
		}
  plot(1:length(N),N, col="blue",type="l",main=paste(title1, r),
			xlab="t", ylab="N[t]")
  }
#You should now have six graphs of the trajectory of the population reaching a stable equilibrium 
#   of K, at growth rates r of 0.4, 0.7, 1.0, 1.3, 1.6, and 1.9.  
#   Let us now try values of r close to 2.0 and see if the behavior indeed changes 
#   where May says it does.  

N[1]<- 10 ;  K<-1000 ; # initial condition and `carrying capacity`
 title1<-"N(t+1) vs N(t) when r=" 
 title2<-"N(t) vs t when r= " # title blanks for plots
 sample_r<-c ( 1.99, 1.999, 2.001, 2.01);
 P_frames<- NULL; R_frames<- NULL;
 par(mfrow=c(2,2))
 
 for (r in sample_r)
 {
	
    for (t in 1:2000)
		{
		#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] );  # logistic model
		N[t+1]<- N[t] * exp(r * (1 - N[t] / K))   #ricker model
		}
  plot(N[1:(length(N)-1)],N[2:length(N)], col="green", main=paste(title1, r),
			xlab="N[t]", ylab="N[t+1]")
  }

  for (r in sample_r)
 {
	
    for (t in 1:2000)
		{
		#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] );  # logistic model
		N[t+1]<- N[t] * exp(r * (1 - N[t] / K))   #ricker model
		}
  plot(1:length(N),N, col="blue",type="l",main=paste(title2, r),
			xlab="t", ylab="N[t]",xlim=c(1900,2000),ylim=c(0,1200))
  }

#We have plotted from time step 1900 to time step 2000 on the graphs to see whether 
# the populations have reached an equilibrium or whether they have stabilized.  
# The model with r=1.999 appears to be a damped oscillation which is still damping, 
# whereas the model with r=2.001 appears to have entered into a stable oscillation.
#Let us now try examining the regions where May says there is a 2 point cycle (2.526 > r > 2.0),
# a 4 point cycle (2.656 > r > 2.526), an 8 point cycle (2.685 > r > 2.656), 
# a cycle of 16 or 32 points (2.692>r>2.685).  
# We will also plot graphs of N(t+1) vs N(t) to see the underlying structure in the data:


N[1]<- 10 ;  K<-1000 ; # initial condition and `carrying capacity`
 title1<-"N(t+1) vs N(t) when r=" 
 title2<-"N(t) vs t when r= " # title blanks for plots
 sample_r<-c ( 2.4,2.6,2.67,2.69);
 P_frames<- NULL; R_frames<- NULL;
 par(mfrow=c(2,2))
 
 for (r in sample_r)
 {
	
    for (t in 1:60)
		{
		#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] );  # logistic model
		N[t+1]<- N[t] * exp(r * (1 - N[t] / K))   #ricker model
		}
  plot(N[1:(length(N)-1)],N[2:length(N)], col="green", main=paste(title1, r),
			xlab="N[t]", ylab="N[t+1]")
  }

  for (r in sample_r)
 {	
	
    for (t in 1:60)
		{
		#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] );  # logistic model
		N[t+1]<- N[t] * exp(r * (1 - N[t] / K))   #ricker model
		}
  plot(1:length(N),N, col="blue",type="l",main=paste(title2, r),
			xlab="t", ylab="N[t]",xlim=c(1,60),ylim=c(0,2000))
  }
  
  
#Now let us examine May's assertion that in the region of chaotic behavior, 
#  the initial conditions determine the pattern of behavior.  
#  Let us try four simulations with r = 3.3, each with a different starting population.
#  The starting conditions used in the second and fourth graphs correspond to May's conditions 
#  of No/K = 0.02 (No = 20) and No/K = 1.5 (No = 1500); see Fig 1(d) and Fig 1(e).  
#  We run the time up to 50 instead of 30 (why stop when you're having fun?).

r<- 3.3; K<- 1000; sample_inits<- c(10, 20, 500, 1500); # parameters and initial conditions
title1<-"N(t+1) vs N(t) when N[0]=" 
 title2<-"N(t) vs t when N[0]= " # title blanks for plots
 for (j in sample_inits)
 {	plotnum<-plotnum+1
	N[1]<-j
    for (t in 1:50)
		{
		#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] );  # logistic model
		N[t+1]<- N[t] * exp(r * (1 - N[t] / K))   #ricker model
		}
  plot(1:length(N),N, col="blue",type="l",main=paste(title1, j),
			xlab="t", ylab="N[t]",xlim=c(1,60),ylim=c(0,3000))
  }
for (j in sample_inits)
 {
	N[1]<-j
    for (t in 1:50)
		{
		#N[t+1]<- N[t] + ( r / k ) * N[t] * ( K - N[t] );  # logistic model
		N[t+1]<- N[t] * exp(r * (1 - N[t] / K))   #ricker model
		}
  plot(N[1:(length(N)-1)],N[2:length(N)], col="green", main=paste(title2, r),
			xlab="N[t]", ylab="N[t+1]")
  }

# Section 3. Bifurcation analysis.

#We recall that May examines the behavior of difference equation analogues of the 
#  logistic equation.  He looks at two equations:
#       1.  N(t+1) = N(t) * exp[r * (1 - N(t)/K]
#       2.  N(t+1) = N(t) * [ 1 + r * (1 - N(t)/K)]
#Both are non-linear equations, and both have a stable equilibrium at N = K, 
#  when 2 > r > 0.  In the long term the population N(t) exhibits cyclic behaviour, 
#  which becomes more pronounced  as the value of the growth parameter r increases. 
#  The extremes of the oscillations, known as stable points, recur periodically. 
#  Here we examine the distribution of stable points as a function of the growth rate 
#  parameter r.  For the purpose of this exercise we will assume the system becomes 
#  stable after 20 time steps, so the last 10 time steps will be used to estimate the 
#  locations of the 'stable points'.  This kind of plot is called a bifurcation plot 
#  
 number_of_steps<- 30;  record_after_step<- 21;
 rvals<-seq(1.0,3.0,0.02)
 par(mfrow=c(1,1))
 for (r in rvals)
	{
     for (t in 1:number_of_steps)
		{
     	#N[t+1]<- N[t] + ( r / K ) * N[t] * ( K - N[t] );  #  discrete logistic model
     	N[t+1]<- N[t] * exp(r * (1-N[t] / K));  # May's eq 1    
#Save the "stable points"  in a data list called plot_list.  Only use data from time steps  21 - 30 in each run.
		}
		 if(r==1.0){ plot(rep(r,10),N[record_after_step:number_of_steps],
						xlab="r",ylab="N",xlim=c(1,3),ylim=c(0,3000),
						main="Bifurcation Map of N vs r")}
		 if(r>1.0){ points(rep(r,10),N[record_after_step:number_of_steps])}
	}