#Project sizes of population growing at finite rate R = 1 +r
#By R. Gomulkiewicz

initial.size <- 10 #starting populatio size

#r <- 0.01 #intrinsic rate of increase

#create a list to hold all the population sizes at all the steps
steps <- 200
n <- matrix(0, steps+1, 2)

#store the initial population size
n[1,] <- initial.size

#this is a function that computes the population size at the next time step
next.size.function <- function(n, r) {n*(1+r)}

#this is a "loop" that projects the population sizes at all steps
#use r = .01
for(i in 1:steps){ n[i+1,1] <- next.size.function(n[i,1],.01) }

#do the loop for the second value of r=-.01
for(i in 1:steps){ n[i+1,2] <- next.size.function(n[i,2],-.01) }


#plot the results for BOTH trajectories
matplot(0:steps,n,type="l",xlab="steps",ylab="population size")

#########################
#plot the relative frequency of "clone 1" versus "clone 2"
freq <- n[,1]/(n[,1]+n[,2])

#evolutionary dynamics of clone 1
plot(0:steps,freq)

#evolutionary dynamics of clone 2
plot(0:steps,1-freq)

##################################
#project sizes of a population with randomness
#this is a function that computes the population size at the next time step 
#assuming each individual leaves a binomially-distributed number of descendants,
# with mean R = 1+r and maximum number of offspring = "max.off"
max.off <-2


next.size.rand <- function(n, max.off,r) {sum(rbinom(n,max.off,(1+r)/max.off))}

reps = 20
nstoch <- matrix(0,steps+1,reps)
nstoch[1,] <- initial.size

for(j in 1:reps){
for(i in 1:steps){ nstoch[i+1,j] <- next.size.rand(nstoch[i,j],max.off,.01) }
}

matplot(0:steps,nstoch,type="l",xlab="steps",ylab="population size")