#############################################################################
# Adaptive dynamics with ecology simulation based on Boots, M. and Y. Haraguchi. 1999. 
# The evolution of costly resistance in host-parasite systems. Am. Nat. 153: 359-370
# Coded by D. George and C. Webb May 28, 2008, updated May 13, 2009
#############################################################################
rm(list=ls()) 	# clear memory

### Parameters #########################################
q= 0.003		# intraspecific competition coefficient
g= 1			# background plus disease mortality
a= 2			# shape parameter for cost function
b= 0.012		# shape parameter for cost function
m= 0.3			# mutation probability
N= 10000		# Number of events that will occur

### Set Initial Values ################
S <- numeric(N)	# create a vector of zeros of length N
W <- S
X <- S
Y <- S
X[1] <- 7		# Initial susceptible population size
Y[1] <- 2		# Initial infected population size
W[1] <- 0		# Initial wait time
loopTime= 0

B <- rep(0.4, X[1])  # Initial susceptibility values  

trait <- array(0, 1000) 
trait[1:X[1]]= B
traitCumulative <- matrix(0, 1000, 1000)

for (i in 2:N){
    # Birth and Death rates for X and Y
    Xborn <- a*sum(B) + b*X[i-1] - q*(X[i-1]+Y[i-1])*X[i-1]
    Xdie <- sum(B)*Y[i-1]
    Yborn <- sum(B)*Y[i-1]
    Ydie <- g*Y[i-1]

    sumrates = Xborn + Xdie + Yborn + Ydie

    if (sumrates <= 0) break
    
    # Determine time to next event from exponential distribution (Sojurn time)
    #S[i] = rexp(1, rate= 1/sumrates) ???
    S[i]= rexp(1, rate=sumrates)
    W[i]= W[i-1] + S[i]   			# Update wait time

    #Decide what event occurred and update populations
    u <- runif(1, min=0, max=1)

    c = length(B)        			# Determine size of B
    
    if (u < (Xborn)/sumrates){ 		# Susceptible born
       X[i]= X[i-1] + 1
       Y[i]= Y[i-1]

       # Decide who gives birth
       s1 = runif(1, min=0, max=1)   # Gives random number to determine event
       selection = 0
       cumtotal = 0
       
		while ( s1 > cumtotal/Xborn){
           selection = selection + 1
           cumtotal = cumtotal + (a*B[selection]+b) - q*(X[i-1]+Y[i-1])
       } # End while (s1 > cumtotal/Xborn)

       # Mutation process
       
		s2 = runif(1, min=0, max=1)
		
		if (s2 < m) {B[c+1] = B[selection] + rnorm(1, mean=0, sd=0.01)}  # Mutation
		else {B[c+1]= B[selection]}		                                 # No mutation
       
       # End susceptible is born
       } 
       
		else{ if (u < (Xborn + Xdie)/sumrates){  # Susceptible dies
          X[i]= X[i-1] - 1
          Y[i]= Y[i-1]
      
          # Decide who dies
          s1 = runif(1, min=0, max=1)			# Gives random number to determine event
          selection = 0
          cumtotal = 0
          while (s1 > cumtotal/Xdie){
                selection = selection + 1
                cumtotal = cumtotal + (B[selection]*Y[i-1])
            }	# End while (s 1> cumtotal/Xdie)
            
            B[selection] = -200
        
		} 	# End susceptible dies

		else{ if (u < (Xborn + Xdie + Yborn)/sumrates){    # Infected born
          X[i]= X[i-1]
          Y[i]= Y[i-1] + 1
    
		} 	# end infected born
    
		else{ if (u >= (Xborn + Xdie + Yborn)/sumrates){   # Infected dies
          X[i]= X[i-1]
          Y[i]= Y[i-1] - 1
		} # End infected dies

    } # end third else if  "infected dies"
    } # end second else if "infected born"
    } # end first else if "susceptible dies"
    
            
   B = B[B > -1] # Shortens B array to the same size as the population (ie remove dead individuals)
   #trait[i,]= rep(0, 1000)
   #trait[i,1:X[i]]= B     # Record trait values
   
   if (length(B) == 0) break         
   traitCumulative[i,]= rep(0,1000) # Line may not be necessary as already zeros, reset to be sure
   traitCumulative[i, 1:X[i]]= B        
   loopTime = loopTime + 1
            
   if (X[i] <= 0) break # end if susceptibles go extinct
   if (Y[i] < 0)  break # end if infecteds go negative (but can be zero)
                                         
} # end for-loop                   


### Plotting ################################
# Clip vectors to appropriate sizes...
clip = length(W[W > 0]) + 1
W= W[1:clip]
X= X[1:clip]
Y= Y[1:clip]

par(mfrow=c(2,1))
plot(c(W,W), c(X,Y), xlab="Time", main="Ecological Dynamics", ylab="Population Size", type="n")
lines(W, X, type="l", col="navy", lwd=1.5)            # Susceptibles
lines(W, Y, type="l", col="dark red", lwd=1.5)        # Infecteds
legend("topleft", c("Susceptibles","Infecteds"), cex=0.7, lwd= rep(1.5,2), col=c("navy","dark red"), bg="gray88")

cut = length(W)
numPlots = seq(from=1, to=cut, by=10)
traitCumulative = traitCumulative[1:cut,]
plot(traitCumulative[,cut], W, type= "n", xlab="Susceptibility, B", ylab="Time", main="Evolutionary Dynamics", col=1, pch=19, xlim= c(0.05,1))
for(p in 1:length(numPlots)){
  points(traitCumulative[,numPlots[p]], W, cex= 1.25, pch= 19, col= rainbow(length(numPlots)))
}