######### a call to fi.sel will
######### generate a cobweb plot of evolutionary dynamics
######### for frequency independent selection at a biallelic locus
######### for n steps starting from x[0] = x0
######### By R. Gomulkiewicz (26 Sep 2013)

fi.sel<- function(wAA = 1, wAa = 1, waa = 0.9, p0 = runif(1,0,1), t = 100,...) 
{
p.prime <- function(wAA, wAa, waa, p){
	q <- 1-p
	w.A <- p*wAA + q*wAa
	w.bar <- p*p*wAA + 2*p*q*wAa + q*q*waa
	p*w.A/w.bar
	}

p<-seq(from=0,to=1,length.out=100)
plot(p, p.prime(wAA, wAa, waa, p),type='l',xlab=expression(p[t]),ylab=expression(p[t+1]), ylim=c(0,1))
abline(a=0,b=1, col = "gray")


start = p0
end = p.prime(wAA, wAa, waa, start)

points(start,start, pch = 19)  #show the initial point
lines(x=c(start,start),y=c(start,end) )  #first draw line from (x0, x0) to (x0, f(x0)

vert=FALSE  #2nd line will be horizonatal from (x0,f(x0) to (f(x0),f(x0))
for(i in 1:(2*t))
{
if(vert) #draw vertical line from (x,x) to (x,f(x))
{
lines(x=c(start,start),y=c(start,end) )
vert=FALSE
}
else  #draw horizontal line from (x,f(x)) to (f(x),f(x))
{
lines(x=c(start,end),y=c(end,end) )
vert=TRUE
start= p.prime(wAA, wAa, waa, start)
end <- p.prime(wAA, wAa, waa, start)
}
}
}

########################################
#### Spread of a recessive advantageous allele
fi.sel(wAA =1.4, wAa = 1, waa = 1, p0 = 0.05)

###############################
#### Overdominant fitnesses
### starting below equilibrium
quartz()
fi.sel(wAA =0.7, wAa = 1, waa = 0.8, p0 = 0.1)
### starting above equilibrium
quartz()
fi.sel(wAA =0.7, wAa = 1, waa = 0.8, p0 = 0.9)

#######################
### Underdominance
### starting below equilibrium
quartz()
fi.sel(wAA =1.2, wAa = 1, waa = 1.3, p0 = 0.59)

### starting above equilibrium
quartz()
fi.sel(wAA =1.2, wAa = 1, waa = 1.3, p0 = 0.61)