# For large matrices it is faster to calculate the dominant eigenvalue and eigenvectors associated with it via iteration rather than using eigen. get.eigen.stuff <- function(mat){ max.it=0 sz <- dim(mat)[1] t.now <- runif(sz) t.now <- t.now/sum(t.now) t.next <- mat%*%t.now t.next <- t.next/sum(t.next) i <- 0 while ((sum(abs(t.next-t.now))>0.0000001)&(i<2000)){ i <- i+1 #print(i) t.now <- t.next t.next <- mat%*%t.now lambda <- sum(t.next)/sum(t.now) if(sum(t.next)!=0){ t.next <- t.next/sum(t.next)} } if (i==2000){lambda=NA t.next=rep(NA,nn*4) r.next=rep(NA,nn*4) max.it=1 print('max iteration for lambda') }else{ r.now <- runif(sz) r.now <- r.now/sum(r.now) r.next <- r.now%*%mat r.next <- r.next/sum(r.next) while (sum(abs(r.next-r.now))>0.0000001){ r.now <- r.next r.next <- r.now%*%mat r.next <- r.next/sum(r.next)} } return(list(lambda,t.next,r.next,max.it)) } #Standardize a matrix such that each column sum to 1 normali<- function(mat){ nr=dim(mat)[1] nc=dim(mat)[2] mat<- ifelse(is.nan(mat),0,mat) mat<-mat/matrix(as.vector(apply(mat,2,sum)),nrow=nr,ncol=nc,byrow=TRUE) mat<- ifelse(is.nan(mat),0,mat) return(mat) } #Calculates the mean and variance of a quantity, typically a continuous character means.and.vars <- function(z,n){ s=ifelse(length(dim(n))==0,1,dim(n)[2]) me=va=rep(0,s) if(s==1){ me <- sum(z*n) / sum(n) va <- sum(z*z*n) / sum(n) - me^2 }else{ for (i in 1:s){ me[i] <- sum(z[,i]*n[,i]) / sum(n[,i]) va[i] <- sum(z[,i]*z[,i]*n[,i]) / sum(n[,i]) - me[i]^2 }} return(matrix(c(me,va),c(2,dim(n)[2]),byrow=T)) } # Compute the kernel component functions from the fitted models # Survival function sx<-function(xbm,a,bbm) { u<-exp(-exp(a+bbm*xbm)) return(u) } # growth function gxy<-function(y,x,K0,L0,d,dbm) { sigmax <- exp(d+dbm*x) sigmax2g <-sigmax^2 mux <- exp(-K0)*(x+(L0)*(1/exp(-K0)-1)) fac1 <- sqrt(2*pi)*sigmax fac2 <- ((y-mux)^2)/(2*sigmax2g) return(exp(-fac2)/fac1) } # inheritance function ixy<-function(y,x,a,bm,d,dbm) { sigmax <- exp(d+dbm*x) sigmax2g <-sigmax^2 mux <- a+bm*x fac1 <- sqrt(2*pi)*sigmax fac2 <- ((y-mux)^2)/(2*sigmax2g) return(exp(-fac2)/fac1) } ##### matrices for each season 1 and 2 #survival S1=diag(sx(meshpoints,sparam1[1],sparam1[2])^12) S2=diag(sx(meshpoints,sparam2[1],sparam2[2])^12) #growth G1 <- outer(meshpoints,meshpoints,gxy,gparam1[1],gparam1[2],gparam1[3],gparam1[4]) G1 <- normali(G1) #probability to produce eggs A=rep(1,nn) A[meshpoints<36]=0 #reproduction R1<- diag((rparam1[1]+rparam1[2]*meshpoints)*A*eggsur*juvsur) R2<- diag((rparam2[1]+rparam2[2]*meshpoints)*A*eggsur*juvsur) #inheritance I <- outer(meshpoints,meshpoints,ixy,inhparam[1],0,inhparam[2],0) I <- normali(I)