# Hunter and Caswell (2005) - Ecological Modelling 188:15-21.
# works with falcon example
require(MASS)

# Set key features of population
p=2 # Number of patches
s=2 # Number of stages

# Build the vec-permutation matrix that is
#  used in many calculations below
P=matrix(0,s*p,s*p)
for (i in 1:s){
  for (j in 1:p){
    E=matrix(0,s,p)
    E[i,j]=1
    E
    P=P+kronecker(E,t(E))
}}
P

# Build the movement matrix
d=.27 # same for young in northern & southern pops
M0=matrix(0,p,p) # empty blocks
M1=matrix(c(1-d,d, d,1-d),p,p,byrow=T) # young
M2=diag(1,p,p)   # adults
M=rbind(cbind(M1,M0),cbind(M0,M2))
M

# Build the demography matrix
fN=.2556; fS=.1908 # fecundity in northern & southern pops
SjN=.72; SjS=.72   # juvenile survival in N. & S. pops
SaN=.77; SaS=.77   # adult survival in N. & S. pops
B0=matrix(0,p,p)   # empty blocks
B1=matrix(c(0,fN, SjN, SaN),p,p,byrow=T) # demography in N.
B2=matrix(c(0,fS, SjS, SaS),p,p,byrow=T) # demography in S.
B=rbind(cbind(B1,B0),cbind(B0,B2))
B

# Build the projection matrix for 'Dispersal before Demography'
# By Patches
A=B%*%t(P)%*%M%*%P
A

# Obtain asymptotic lambda for B1, B2, A
(LambdaB1<-max(Re(eigen(B1)$values)))
(LambdaB2<-max(Re(eigen(B2)$values)))
(LambdaA<-max(Re(eigen(A)$values)))

# Obtain the Sensitivity & Elasticity matrices
eig<-eigen(A)
EigVecs<-eig$vectors # used later
Lambdas<-Re(eig$values)
pos<-which.max(Lambdas)
V <- Conj(ginv(EigVecs))
w <- Re( t(t(EigVecs[,pos])))
v <- Re(t(t(V[pos,])))

# Sensitivity matrix for A
SensA <- v%*%t(w) # t(w) = transpose of 'w'
round(SensA,4)

# Elasticity matrix for A
tmp <- SensA*A
ElastA <- tmp/LambdaA
round(ElastA,4)

# Sensitivity matrix for B
SensB=SensA%*%t(P)%*%t(M)%*%P
round(SensB,4)

# Elasticity matrix for B
ElastB=B*SensB/LambdaA
round(ElastB,4)

#Sum diagonal blocks of Elast. B
ElastB1=sum(ElastB[1:2,1:2])
round(ElastB1,2)
ElastB2=sum(ElastB[3:4,3:4])
round(ElastB2,2)

# Sensitivity matrix for M
SensM=P%*%t(B)%*%SensA%*%t(P)
round(SensM,4)

# Elasticity matrix for M
ElastM=M*SensM/LambdaA
round(ElastM,4)


# Alternative representation of 'A',
#  which is like what's seen in Lebreton 1996
# Here, the matrix is organized by stages
#  upper left block  = F for juveniles
#  upper right block = F for adults
#  lower left block  = surv. & moves for juv.
#  lower right block = surv. & moves for ads.
AltA=P%*%B%*%t(P)%*%M
round(AltA,3)
round(A,3)

#########################################################
##############  Population Projection ###################
#########################################################

# Set up variables for projecting the population forward
# Set initial values for age classes
#  & store in matrix N in which
#  rows=stages & columns=patches
N=matrix(data=c(40,200, 10,50),nrow=s,ncol=p,byrow=F)

# create vectors that have n for patches & stages
# stored together by (a) patch=np or (b) stage=ns
np=matrix(N,nrow=p*s,ncol=1)
ns=matrix(t(N),nrow=p*s,ncol=1)

# Project the population forward 3 & 4 years and
#  convert the output to a block matrix
#  that is called either NfinalA or NfinalAltA &
#  whose rows  = stages &
#     columns  = patches
(NAyr4=A%*%A%*%A%*%np)              # data arranged by patches
(NAyr5=A%*%A%*%A%*%A%*%np)          # data arranged by patches
(NfinalA=matrix(NAyr5,2,2,byrow=F)) # put data into matrix

(NAlt4=AltA%*%AltA%*%AltA%*%ns)               # data arranged by stages
(NAlt5=AltA%*%AltA%*%AltA%*%AltA%*%ns)        # data arranged by stages
(NfinalAltA=matrix(NAlt5,2,2,byrow=T)) # put data into matrix

sum(NAlt5)/sum(NAlt4)
LambdaA