###################################################
### chunk number 1: CandG
###################################################
setwd("~/Documents/Projects/Book/draft/Chap04")
rm(list=ls())
badfonts = FALSE
.detach.stuff = function() {
 s1 = grep("package|Autoloads",search())
 nattach = length(search())
 xattach = search()[c(-1,-s1)]
 for (i in xattach)
   eval(substitute(detach(i),list(i=i)))
}
.detach.stuff()
## ps.options(family="NimbusSan")
ps.options(colormodel="cmyk")
options(width=75, SweaveHooks=list(fig=function() par(mar=c(5,4,.5,1.5))) ) 

## palette(gray((0:8)/8))
library(lattice)
trellis.par.set(canonical.theme(color=FALSE))
myQ <- function(w=4,h=4, mr = c(5,4,.5,1.5)) {
  quartz(width=w,height=h)
  par(mar=mr)
}
CandG <- c(32,16,10,9,8,7,8,6,4,5,7,4,2,1,1,1,3,8,12)
NoS <- 1:19
par(mar=c(3,3,1,1), mgp=c(1.5,.5,0))
bp <- barplot(CandG,ylim=c(0,35), 
        ylab="No. of Species", xlab="Number of Sites Occupied")
axis(1, at=bp[seq(1,19, by=2)], labels=seq(1,19, by=2) )



###################################################
### chunk number 2: 
###################################################
L1 <- 2; L2 <- .4
A <- matrix(c(1,    0, 
              L1-1, L2), nrow=2, byrow=TRUE)


###################################################
### chunk number 3: 
###################################################
eigen(A)


###################################################
### chunk number 4: sourceP
###################################################
L1s <- seq(1, 3, by=0.01)
p1 <- sapply(L1s, function(l1) {A[2,1] <- l1-1;
                                eigen(A)$vectors[1,1]/sum(eigen(A)$vectors[,1])
                                })
plot(L1s, p1, type='l', ylab="Source Population", xlab=expression(lambda[1]) ) 


###################################################
### chunk number 5: 
###################################################
levins <- function(t, y, parms) {
       p <- y[1]
       with(as.list(parms), {
         ## Here is the equation
       dp <- ci * p * (1-p) - e * p
       return(list(dp))
       }) }


###################################################
### chunk number 6: 
###################################################
library(deSolve)
prms <- c(ci=0.15, e=0.05)
Initial.p <- 0.01
out.L <- data.frame(ode(y=Initial.p, times=1:100, func=levins, parms=prms ))


###################################################
### chunk number 7: levins
###################################################
plot(out.L[,2] ~ out.L[,1],   type='l', ylim=c(0,1), ylab="p", xlab="time")  


###################################################
### chunk number 8: 
###################################################
gotelli <- function( t, y, parms) {
       p <- y[1]
       with(as.list(parms), {
         ## Here is the equation
       dp <- ce *  (1-p) - e * p
       return(list(dp))
       })
       }



###################################################
### chunk number 9: 
###################################################
hanski <- function( t, y, parms) {
       p <- y[1]
       with(as.list(parms), {
         ## Here is the equation
       dp <- ci * p *  (1-p) - e * p *(1-p)
       return(list(dp))
       })
       }


###################################################
### chunk number 10: 
###################################################
prms <- c(ci<- 0.15, ce <- 0.15, e=0.05)
out.IMH <- data.frame(ode(y=Initial.p, times=1:100, func=gotelli, parms=prms ))
out.IMH[["pH"]]<- ode(y=Initial.p, times=1:100, func=hanski, parms=prms )[,2]


###################################################
### chunk number 11: IMH
###################################################
matplot(out.IMH[,1], out.IMH[,2:3], type='l', col=1, ylab="p", xlab="time")  
legend("topleft", c("Hanski", "Propagule Rain"), lty=2:1, bty='n')


###################################################
### chunk number 12: 
###################################################
dpdtCS <- expression((ci-e)*p*(1-p))
ci <- 0.15; e <- 0.05
p <- seq(0,1, length=50)


###################################################
### chunk number 13: hanski
###################################################
plot(p, eval(dpdtCS), type='l', ylab="dp/dt")


###################################################
### chunk number 14: 
###################################################
lande <- function(t, y, parms) {
       p <- y[1]
       with(as.list(parms), {
         ## Here is the equation
       dp <- ci * p * (1-D-p) - e * p
       return(list(dp))
       }) }


###################################################
### chunk number 15: 
###################################################
library(deSolve)
prmsD <- c(ci=0.15, e=0.05, D=0)
Ds <- c(0,.2,.5)
Initial.p <- 0.01
t <- 1:200


###################################################
### chunk number 16: 
###################################################
ps <- sapply(Ds, function(d) {
  prmsD["D"] <- d
  ode(y=Initial.p, times=t, func=lande, parms=prmsD )[,2]
})


###################################################
### chunk number 17: lande
###################################################
matplot(t, ps, type='l', ylab="p", xlab="time")
text(c(200,200,200), ps[200,], paste("D = ", Ds, sep=""),
     adj=c(1,0))


###################################################
### chunk number 18: 
###################################################
C1 <- ode(y=0.999, times=t, func=hanski, parms=c(ci=0.2, e=0.01))
C2 <- ode(y=0.999, times=t, func=hanski, parms=c(ci=0.2, e=0.25))
L2 <- ode(y=0.95, times=t, func=levins, parms=c(ci=0.2, e=0.25))


###################################################
### chunk number 19: lande2
###################################################
matplot(t, cbind(C1[,2],C2[,2],L2[,2]), type='l', ylab="p", xlab="Time", col=1)
legend('right', c("c > e", "c < e", "c < e (Levins)"), lty = 1:3, bty='n' )


###################################################
### chunk number 20: 
###################################################
args(MetaSim)


###################################################
### chunk number 21: CoreSatSim10
###################################################
out.CS.10 <- MetaSim(method="hanski", NSims=10)
matplot(out.CS.10$t, out.CS.10$Ns, type="l", xlab="Time", ylab="Occupancy", 
        sub=out.CS.10$method)



###################################################
### chunk number 22: 
###################################################
system.time(
out.CS.Lots <- MetaSim(method="hanski",NSims=50, Time=1000)
)


###################################################
### chunk number 23: CoreSatSimHist
###################################################
hist(out.CS.Lots$Ns[501,], breaks=10, main=NULL,
     xlab=expression("Occupancy ("*italic("p")*")"), ylab="Number of Species",sub=paste(out.CS.Lots$method, " Model", sep=""))



###################################################
### chunk number 24: 
###################################################
system.time(
out.L.Lots <- MetaSim(NSims=50, Time=500, method="levins")
)


###################################################
### chunk number 25: LevinsSimHist
###################################################
hist(out.L.Lots$Ns[501,], breaks=10,
     xlab=expression("Occupancy ("*italic("p")*")"), 
     ylab="Number of Species", 
     main=NULL, sub=paste(out.L.Lots$method, " Model", sep=""))