function evolve4()

evolve4= function( n.pop1=1000, n.pop2=1000, n.gen=60, n.samp=100, f.init=0.05, s.adv=0.2, dom = FALSE, histogram=FALSE, plot.mean=TRUE, color=FALSE, title = "") 

{

#Based on code by Lex Flagel in his blog A Hopeful Monster

#

#      -allows for a linear time-varying population

#      -uses simplified forms of the Hardy-Weinberg expressions for calculating the change of allele frequency  (see Falconer and Mackay 4th ed., p.29, equations 2.7 and 2.8)

#      -case of either no dominance or complete dominance

#      -assumes: s.adv is small, or allele frequency is small.  Usually I use 0 < s.adv < 0.2, as in F&M page 29.

#  Usage: 

#     Recreate Fig 2.3 (upper panel) in F&M

#        evolve4(f.init=0.95, s.adv= -0.2, title="No dominance, neg. advantage", histogram=TRUE, color=TRUE) 

#       evolve4(f.init=0.05, s.adv=  0.2, title="No dominance, pos. advantage", histogram=TRUE, color=FALSE)  

#       evolve4(f.init=0.005, s.adv=  0.2, n.gen=150, title="No dominance, pos. advantage", histogram=TRUE, color=TRUE) 

#       Drift in 5 generations#     evolve4(n.pop1=100, n.pop2=100, n.gen=5, n.samp=1000, f.init=0.5, s.adv=0.0)

#       evolve4(f.init=0.5, s.adv=  0.1, n.pop1=200, n.pop2=200, n.gen=200, title="No dominance, pos. advantage", histogram=TRUE, color=TRUE, dom=FALSE) 

#      evolve4(f.init=0.5, s.adv=  -0.1, n.pop1=1000, n.pop2=1000, n.gen=200, title="No dominance, neg. advantage", histogram=TRUE, color=TRUE, dom=FALSE) 

#      evolve4(f.init=0.1, s.adv=  0.005, n.pop1=5000, n.pop2=5000, n.gen=200, title="", histogram=FALSE, color=TRUE, dom=FALSE) 

# n.pop1 = starting population size

# n.pop2 = ending population size

# n.gen = number of generations to run

# n.samp = number of individual lines (monte carlo realizations) to calculate

# s.adv = reproductive advantage (=0.0 for just random drift)

# f.init = initial allele frequency in population

#

    k = matrix(-999, nr=n.samp, nc=n.gen-1) #matrix to store results

#  graphics.off()

#  setwd("/Users/paulreasenberg/Desktop/evolution/plots")

#

for (i in 1:n.samp)

 {

#  Falcon & Mackay parameters

      q = f.init  # allele frequency

      p = 1-q

      s = s.adv

#

   for (j in 1:(n.gen-1) )   # n.gen generations require n.gen - 1 samples

{     

      dq = 0.5 * s * q * p                  #  No dominance

      if (dom) {dq = s * q * q * p}      #  Complete dominance

      prob= min((q + dq), 1.0)    #truncate prob from above

      prob= max(prob, 0.0)         #truncate prob from below

      sz = as.integer( 2* (n.pop1 + ( (n.pop2-n.pop1) *j / (n.gen-1) ) ) )  

      q = rbinom(1, size=sz, prob) / sz

      p = 1 - q

#      

      k[i, j] = q  #update matrix

}

}

#

k = cbind(rep(f.init, n.samp), k) #bind first gen. first column of matrix

#

#

#######

##       that’s it for sampling, code below is plotting

#

# Plot all n.samp random walks

dev.new(width=5,height=4)

#jpeg(“randomwalks.jpg”)

plot(1,1, xlim=c(1,n.gen), ylim=c(0,1), ty='n', xlab='generation', ylab='allele freq')

abline(h=0.5) #

mtext(line=0, cex=0.7, paste('pop=', n.pop1, n.pop2, 'gen=', n.gen, 'samp=', n.samp, 'f.init=', f.init, 'adv=', s.adv, 'dom=', dom ))

mtext(line=1, cex=1.0, title)

#

for (i in 1:n.samp){

ic=1

if (color) {ic = round(502 * i/n.samp, digits=0)}

  lines(1:n.gen, k[i, ], lwd=.3, col=ic)  #draw n.samp random walks as lines in color col=ic

  points(n.gen+0.1, k[i, n.gen-1], pch='-'    )  #add red rug on right

}

#

if (histogram) {

#  Histogram of final allele frequencies

dev.new(width=5,height=4)

hist (k[, n.gen], xlim=c(0,1),  main='', xlab="allele freq") 

mtext(line=0, cex=0.7, paste('Gen', n.gen, 'mean=', round(mean(k[, n.gen]),3), 'sigma=', round(sd(k[, n.gen]),3) ) )

mtext(line=1, cex=1, "Final allele frequency distribution" )

}

if (plot.mean)  {

# Plot the mean and std. dev. as function of generation

dev.new(width=5,height=4)

q.mean =  rep(-999, n.gen)

q.std =   rep(-999, n.gen)

q.sqmean = rep(-999, n.gen)

for (j in 1:(n.gen))

   {    # accumulate mean and sd of freq at each generation

          q.mean[j] = mean(k[,j])

          q.sqmean[j] = mean( (k[,j])*(k[,j]) )

          q.std[j] = sd(k[,j])

   }

plot(1,1, xlim=c(1,n.gen-1), ylim=c(0, 1 ), ty='n'  , xlab='Generation', ylab='Means of q and q^2')

lines (1:n.gen,  q.mean,  lwd=2)

lines (1:n.gen,  q.mean-q.std, lwd=0.5, lty=3)

lines (1:n.gen,  q.mean+q.std, lwd=0.5, lty=3)

lines (1:n.gen,  q.sqmean, lwd=2)

}

} ###### END 

###############################

#  Turn off all graphics windows

#  graphics.off()

#