###############################################################################

Model <- function(  # all available actions with a model

         A,         # projection matrix
         tmax=20    # length of projection
         ){

  MP <- MatProp(A)
  print(MP)
  Pr <- Projection(A,tmax=tmax)
  plotPCensus(Pr)
  plotPStructure(Pr)
  plotEigen(MP)
  plotSensitivity(A,MP$S,MP$E)
}


###############################################################################
###############################################################################

Square <- function(A){ # Test whether matrix A is a square one

    s <- sqrt(floor(length(A)))
    {if (!(length(A) == s^2)){
        print("Warning! The matrix is not a square one.",quote=FALSE)
        R <- FALSE}
     else{R <- TRUE}}
    R}


###############################################################################

MatProp <- function(A){ # Check properties of matrix A.

if (Square(A)){
   s <- sqrt(floor(length(A)))
   {if (s > 1){
        Ap <- 1*(A>0)
        for (i in 1:s){Ap[i,i] <- 1}
        App <- Ap
        if (s > 2){for (i in 1:(s-2)){App <- 1*((App %*% Ap)>0)}}
        reducible <- (sum(App) < s*s)

        {if (!reducible){
           Ap <- 1*(A>0)
           App <- Ap
           for (i in 1:((s-1)^2)){App <- 1*((App %*% Ap)>0)}
           primitive <- !(sum(App) < s*s)
           }
           else{primitive <- FALSE}
        }
      }
      else{primitive <- (sum(A) > 0)
           reducible <- !primitive}
   }

   EV <- eigen(t(A),symmetric=FALSE)
   EVval <- EV$values
   EVvec <- EV$vectors
   j <- 1
   if (!(primitive)){
     repeat{if ((Im(EVval[j]) > pi/s) | (Re(EVval[j]) < 0)){j <- j+1}
              else{break}}
     }
   v <- abs(Re(EVvec[,j]))
   norm <- sum(v)
   v <- v/norm                         # right eigenvector (reproductive value)

   EV <- eigen(A,symmetric=FALSE)
   EVval <- EV$values
   EVvec <- EV$vectors
   j <- 1
   if (!(primitive)){
     repeat{if ((Im(EVval[j]) > pi/s) | (Re(EVval[j]) < 0)){j <- j+1}
              else{break}}
     }

   lambda <- abs(EVval[j])             # dominant eigenvalue
   w <- abs(Re(EVvec[,j]))
   norm <- sum(w)
   w <- w/norm                         # left eigenvector (stable structure)

   d <- 1
   repeat{
     if (d == s) {break}
       else{if (abs(1-abs(EVval[d+1])/lambda) < 1e-15){d <- d+1}
              else{break}}
     }                                  # period of matrix

   rho <- Inf                           # damping ratio
   if (d < s){rho <- lambda/abs(EVval[d+1])}

   S <- (v %*% t(w))/sum(w * v)         # sensitivity matrix
   E <- S*A/lambda                      # elasticity matrix

   MatProp <- list(EVval = EVval,       # eigenvalues
                   EVvec = EVvec,       # eigenvectors
                   S = S,               # sensitivity matrix (lambda1 on A)
                   E = E,               # elasticity matrix
                   reducible = reducible,
                   primitive = primitive,
                   d = d,              # period
                   w = w,              # stable structure
                   v =v,               # reproductive value
                   rho=rho,            # damping ratio
                   lambda1 = lambda    # dominant eigenvalue
                                      )
   }
}

###############################################################################

Projection <- function( # calculate projection

     A,                 # projection matrix
     n0 = NULL,         # initial value
     tmax = 10          # length of projection
     ){

if (Square(A)){
   s <- sqrt(floor(length(A)))
   if (length(n0) == 0){n0 <- array(0,s); n0[1] <- 1}

   nt <- array(NA,c(s,tmax+1))
   nt[,1] <- n0
   for (i in 1:tmax){nt[,i+1] <- A %*% nt[,i]}
   Projection <- nt
   }
}


###############################################################################

plotPCensus <- function(nt){# plot projected census

  s <- length(nt[,1])
  tmax <- length(nt[1,])-1
  cols <- rainbow(s)
  census <- array(NA,tmax+1)
  for (i in 1:(tmax+1)){census[i] <- sum(nt[,i])}

  st <- array(0,c(s,tmax+1))
  st[1,] <- nt[1,]
  for (i in 2:s){st[i,] <- st[i-1,]+nt[i,]}
  pt <- array(0,2*tmax+1)
  tp <- array(0,2*tmax+1)
  tp[1:(tmax+1)] <- seq(0,tmax)
  for (j in 1:tmax){tp[2*tmax+2-j] <- j}

  X11()
  plot(c(0,tmax),c(0,max(census)),type="n",xlab="Time",ylab="Numbers",bty="l",
       main="Census with stage or age structure")
  polygon(c(0,seq(0,tmax),tmax),c(0,nt[1,],0),density=20,
            angle=45,col=cols[1],lty=3)
  for (i in 2:s){
    for (j in 1:tmax){pt[j] <- st[i,j]
                      pt[2*tmax+2-j] <- st[i-1,j+1]}
    pt[tmax+1] <- st[i,tmax+1]
    polygon(tp,pt,density=20,angle=45+(i-1)*180/s,col=cols[i],lty=3)
    }
  points(seq(0,tmax),census,type="l")
  for (i in 1:(s-1)){points(seq(0,tmax),st[i,],type="l")}
}


###############################################################################

plotPStructure <- function(nt){# plot projected age or stage structure

  s <- length(nt[,1])
  tmax <- length(nt[1,])-1
  cols <- rainbow(s)
  census <- array(NA,tmax+1)
  for (i in 1:(tmax+1)){census[i] <- sum(nt[,i])}


  X11()
  split.screen(c(2,1))
  screen(1)
  plot(c(0,1.02*tmax),c(0,max(nt)),type="n",xlab="Time",ylab="Numbers",
       bty="l",main="Structure")
  for (i in 1:s){points(seq(0,tmax),nt[i,],type="l",col=cols[i])
                 text(1.02*tmax,nt[i,tmax+1],i,col=cols[i])
                 }
  screen(2)
  ft <- array(NA,c(s,tmax+1))
  for (i in 1:s){ft[i,] <- nt[i,]/census}
  plot(c(0,1.02*tmax),c(0,1.15*max(ft)),type="n",
       xlab="Time",ylab="Frequencies",bty="l")
  for (i in 1:s){points(seq(0,tmax),nt[i,]/census,type="l",col=cols[i])
                 text(1.02*tmax,nt[i,tmax+1]/census[tmax+1],i,col=cols[i])
                 }
  close.screen(all=TRUE)
}


###############################################################################

plotEigen <- function(  # plot eigenvalues and dominant eigenvectors
             MP  # list of matrix properties (output of function MatProp)
             ){
  s <- length(MP$w)
  cols <- rainbow(s)

  X11()
  split.screen(c(2,1))

  screen(1)
  xl <- max(1,MP$lambda1)
  plot(51*c(-xl,xl)/17,c(-xl,xl),type="n",main="Eigenvalues",
      xlab="Real part",ylab="Imaginary part",bty="n")
  xp <- seq(-1,1,0.002)
  points(xp, sqrt(1-xp^2),type="l",lty=3,lwd=0.5,col="red")
  points(xp,-sqrt(1-xp^2),type="l",lty=3,lwd=0.5,col="red")
  xp <- seq(-MP$lambda1,MP$lambda1,2*MP$lambda1/1000)
  points(xp, sqrt(MP$lambda1^2-xp^2),type="l",lty=3,lwd=0.5)
  points(xp,-sqrt(MP$lambda1^2-xp^2),type="l",lty=3,lwd=0.5)
  points(0,0,pch="|",col="blue")
  points(Re(MP$EVval),Im(MP$EVval),pch=19,col="blue")
  for(i in 1:length(MP$EVval)){
    points(c(0,Re(MP$EVval[i])),c(0,Im(MP$EVval[i])),type="l",col="blue")}

  split.screen(c(1,2),screen=2)
  screen(3)
  plot(c(0,s),c(0,max(MP$w)),type="n",bty="n",xaxt="n",
       main="Right eigenvector\n(stable stage distribution)",xlab="",ylab="")
  for (i in 1:s){
    polygon(c(i-1,i-1,i,i),c(0,MP$w[i],MP$w[i],0),density=20,
            angle=45+(i-1)*180/s,col=cols[i],lty=3)
    points(c(i-1,i-1),c(0,MP$w[i]),type="l")
    points(c(i-1,i),c(MP$w[i],MP$w[i]),type="l")
    points(c(i,i),c(MP$w[i],0),type="l")
    text(i-0.5,0,i)
    }
  screen(4)
  plot(c(0,s),c(0,max(MP$v)),type="n",bty="n",xaxt="n",
       main="Left eigenvector\n(reproductive values)",xlab="",ylab="")
  for (i in 1:s){
    polygon(c(i-1,i-1,i,i),c(0,MP$v[i],MP$v[i],0),density=20,
            angle=45+(i-1)*180/s,col=cols[i],lty=3)
    points(c(i-1,i-1),c(0,MP$v[i]),type="l")
    points(c(i-1,i),c(MP$v[i],MP$v[i]),type="l")
    points(c(i,i),c(MP$v[i],0),type="l")
    text(i-0.5,0,i)
    }

  close.screen(all=TRUE)
}


###############################################################################

plotSensitivity <- function(# plot image of sensitivity and elasticity matrices

        A,                  # projection matrix
        S=NULL,             # matrix of sensitivities
        E=Null              # matrix of elasticities
        ){

  if (length(S) == 0){MP <- MatProp(A)
                      S <- MP$S
                      E <- MP$E}

  s <- floor(sqrt(length(A)))
  k <- 1+floor((s-1)/10)

  plotImageAndGauge <- function(Mat,cols,Ttl){
    NoCol <- length(cols)
    plot(c(0,51*(s+k)/17),c(-k,s),type="n",bty="n",
         xlab="",ylab="",xaxt="n",yaxt="n",
         main=paste(Ttl,"of the growth rate to the projection matrix entries"))

    mx <- max(Mat)
    L <- floor(log10(mx))
    {if (L > 2){f <- 0; r <- L+1}
        else{f <- 2-L
             {if (L < 0){r <- f+2}
                 else{r <- 5}
                 }
        }
    }
    form <- paste("%",r,".",f,"f",sep="")    

    text((2.55+0.06*r)*s,0,sprintf(form,0),adj=c(1,0.5))
    for (i in 1:NoCol){
      polygon(c(2.1*s,2.5*s,2.5*s,2.1*s),
              c((i-1)*s/NoCol,(i-1)*s/NoCol,i*s/NoCol,i*s/NoCol),
              col=cols[i])
      text((2.55+0.06*r)*s,i*s/NoCol,sprintf(form,i*mx/NoCol),adj=c(1,0.5))
      }
    nk <- floor(s/k)
    for (i in 1:nk){
      text(1-s*0.05,s-i*k+0.5,i*k,adj=c(1,0.5))
      text(i*k+0.5,-s*0.05,i*k,adj=c(0.5,1))
      }
    for (i in 1:s){
      for (j in 1:s){
         colNo <- 1+floor((NoCol-1)*Mat[i,j]*1.01/mx)
         {if (A[i,j] == 0){den <- 20} else{den <- NULL}}
         polygon(c(j,j+1,j+1,j),c(s-i,s-i,s+1-i,s+1-i),
                 col=cols[colNo],density=den,border="black")
         }
      }
  }

  X11()
  split.screen(c(2,1))
  screen(1)
  plotImageAndGauge(S,topo.colors(6),"Sensitivity")
  screen(2)
  plotImageAndGauge(E,heat.colors(5)[seq(5,1,-1)],"Elasticity")
  close.screen(all=TRUE)
}


###############################################################################

# some examples - functions return a projection matrix

KillerWhale <- function(){
A <- array(NA,c(4,4))
A[1,] <- c(0,      0.0043, 0.1132, 0     )
A[2,] <- c(0.9775, 0.9111, 0,      0     )
A[3,] <- c(0,      0.0736, 0.9534, 0     )
A[4,] <- c(0,      0,      0.0452, 0.9804)
A}


USApopulation <- function(){
A <- array(0,c(10,10))
A[1,]<-c(0,     .00102,.08515,.30574,.40002,.28061,.15260,.06420,.01483,.00089)
A[2,1] <-.99670
A[3,2] <-       .99837
A[4,3] <-              .99780
A[5,4] <-                     .99672
A[6,5] <-                            .99607
A[7,6] <-                                   .99472
A[8,7] <-                                          .99240
A[9,8] <-                                                 .98867
A[10,9]<-                                                        .98274
A}

Teasel <- function(){
A <- array(NA,c(6,6))
A[1,] <- c(0,     0,     0,     0,     0,   322.38 )
A[2,] <- c(0.966, 0,     0,     0,     0,     0    )
A[3,] <- c(0.013, 0.010, 0.125, 0,     0,     3.448)
A[4,] <- c(0.007, 0,     0.125, 0.238, 0,    30.170)
A[5,] <- c(0.001, 0,     0.036, 0.245, 0.167, 0.863)
A[6,] <- c(0,     0,     0,     0.023, 0.750, 0    )
A}

TwoStage <- function(sigma1,sigma2,gamma,phi){
A <- array(NA,c(2,2))
A[1,] <- c(sigma1*(1-gamma),    phi)
A[2,] <- c(    sigma1*gamma, sigma2)
A}