################################################### ### chunk number 1: ################################################### setwd("~/Documents/Projects/Book/draft/Chap07") 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)) #source("DMBpplane.R") require(primer); ################################################### ### chunk number 2: MayCriteria1 ################################################### S <- 40; C=.3 curve(C*x^ -2, .1, 1, ylab=expression("Number of Species ("*italic(S)*")"), xlab = expression("Interaction Strength ("*italic(I)*")") ) text(.6, 20, "Unstable Webs", srt=45) text(.2, 3.75, "Stable Webs" , srt=-45, cex=.8) ################################################### ### chunk number 3: ################################################### Aq = matrix(c( -1, -10, 0, 0, 0.1, 0, -10, 0, 0, 0.1, 0, -10, 0, 0, 0.1, 0), nrow=4, byrow=TRUE) ################################################### ### chunk number 4: ################################################### S <- nrow(Aq) ################################################### ### chunk number 5: ################################################### M <- Aq * runif(S^2) ################################################### ### chunk number 6: ################################################### eM <- eigen(M)[["values"]] ################################################### ### chunk number 7: ################################################### deM <- max( Re(eM) ) ################################################### ### chunk number 8: ################################################### intraNDD <- sqrt(sum(diag(M)^2)/S) ################################################### ### chunk number 9: ################################################### diag(M) <- 0 IS <- sqrt( sum(M^2)/(S*(S-1)) ) ################################################### ### chunk number 10: PL ################################################### pimmlawton <- function(mat, N=1, omni.i=NA, omni.j=NA, omega=NULL){ S <- nrow(mat) if(is.na(omni.i)) { out <- matrix(NA, nrow=N, ncol=4) colnames(out) <- c("DomEig", "Im", "IntraDD", "I") for(n in 1:N) out[n,] <- { M = runif(S^2) * mat ## Do eigenanalysis eigs <- eigen(M)[["values"]] mx <- which.max(Re(eigs)) deM = Re(eigs)[mx] deMi = Im(eigs)[mx] intraNDD <- sqrt(sum(diag(M)^2)/S) diag(M) <- 0 IS = sqrt( sum(M^2)/(S*(S-1)) ) c(deM, deMi, intraNDD, IS) } } else { out <- matrix(NA, nrow=N, ncol=5) colnames(out) <- c("DomEig", "Im", "IntraDD", "I", "I.omni") for(n in 1:N) out[n,] <- { M = runif(S^2) * mat ## Adjust for McCann type omnivory if(!is.null(omega)) {M[omni.i,omni.j] <- omega*M[omni.i+1,omni.j] M[omni.i+1,omni.j] <- (1-omega)*M[omni.i+1,omni.j] M[omni.j,omni.i] <- omega*M[omni.j,omni.i+1] M[omni.j,omni.i+1] <- (1-omega)*M[omni.j,omni.i+1]} ## Do eigenanalysis eigs <- eigen(M)[["values"]] mx <- which.max(Re(eigs)) deM = Re(eigs)[mx] deMi = Im(eigs)[mx] intraNDD <- sqrt(sum(diag(M)^2)/S) diag(M) <- 0 IS = sqrt( sum(M^2)/(S*(S-1)) ) omnivory <- sqrt(mean(c(M[omni.i,omni.j],M[omni.j,omni.i])^2)) c(deM, deMi,intraNDD, IS, omnivory) } } return(as.data.frame(out)) } ################################################### ### chunk number 11: ################################################### args(pimmlawton) ################################################### ### chunk number 12: ################################################### set.seed(1) pimmlawton(Aq) ################################################### ### chunk number 13: ################################################### out.A <- pimmlawton(Aq, N=2000) ################################################### ### chunk number 14: ################################################### summary(out.A) ################################################### ### chunk number 15: pairsA ################################################### pairs(out.A) ################################################### ### chunk number 16: ################################################### RT.A <- -1/out.A[["DomEig"]] summary(RT.A) ################################################### ### chunk number 17: ################################################### sum( RT.A > 150 ) / 2000 ################################################### ### chunk number 18: histA ################################################### A.fast <-RT.A[RT.A < 150] A.slow <-RT.A[RT.A >= 150] # Opens a new plotting device. On Macs use quartz() or x11() histA <- hist(A.fast, breaks=seq(0,150, by=5), main=NULL) ################################################### ### chunk number 19: ################################################### Eq = matrix(c( -1, 0, 0, -10, 0, -1, 0, -10, 0, 0, -1, -10, 0.1, 0.1, 0.1, 0), nrow=4, byrow=TRUE) ################################################### ### chunk number 20: ################################################### out.E <- pimmlawton(Eq, N=2000) summary(out.E) ################################################### ### chunk number 21: E1 ################################################### layout(matrix(1:2,nr=1)) plot(DomEig ~ IntraDD, data=out.E) plot(DomEig ~ I, data=out.E) ################################################### ### chunk number 22: histE ################################################### RT.E <- -1/out.E[["DomEig"]] E.fast <-RT.E[RT.E < 150] E.slow <-RT.E[RT.E >= 150] # Opens a new plotting device. On Macs use quartz() or x11() histE <- hist(E.fast, breaks=seq(0,150, by=5), main=NULL) ################################################### ### chunk number 23: ################################################### Bq = matrix(c( -1, -10, 0, 0, 0.1, 0, -10, -10, 0, 0.1, 0, -10, 0, 0.1, 0.1, 0), nrow=4, byrow=TRUE) ################################################### ### chunk number 24: ################################################### out.B <- pimmlawton(Bq, N=2000, omni.i=2, omni.j=4) summary(out.B) ################################################### ### chunk number 25: pairsB ################################################### pairs(out.B) ################################################### ### chunk number 26: histB ################################################### RT.B <- -1/out.B[["DomEig"]] B.fast <- RT.B[RT.B < 150 & RT.B>0 ] out.B.fast <- out.B[RT.B < 150 & RT.B>0 ,] out.B.stab <- out.B[RT.B>0 ,] histB <- hist(B.fast, breaks=seq(0,150, by=5), main=NULL) ################################################### ### chunk number 27: histsAB ################################################### hist(log(RT.A,10), probability=T, ylim=c(0,1), main=NULL, xlab=expression(log[10]("Return Time")) ) lines(density(log(RT.A,10)) ) lines(density(log(RT.B[RT.B>0],10)), lty=2, lwd=2 ) legend("topright", c("Chain A", "Chain B"), lty=1:2, lwd=1:2, bty='n')