# Constant matrix

A <- matrix(c(0,0.3,0,  1,0,0.5,  5,0,0),c(3,3))

T <- 50    # Time of projection

n <- array(NA,c(3,T+1))
N <- array(NA,T+1)

n[,1] <- c(1,0,0)   # Initial condition
N[1] <- sum(n[,1])
for (i in 1:T){
    n[,i+1] <- A %*% n[,i]
    N[i+1] <- sum(n[,i+1])
}

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

screen(1)
plot(NA,NA,bty="l",xlim=c(0,T+1),ylim=c(0.001,max(n)),xlab="t",ylab="n",log="y")
points(seq(0,T),n[1,],type="l",lwd=2)
points(seq(0,T),n[2,],type="l",lwd=2,col="red")
points(seq(0,T),n[3,],type="l",lwd=2,col="blue")
text(T+1.5,n[1,T+1],expression(n[1]))
text(T+1.5,n[2,T+1],expression(n[2]),col="red")
text(T+1.5,n[3,T+1],expression(n[3]),col="blue")

screen(2)
plot(NA,NA,bty="l",xlim=c(0,T+1),ylim=c(0,1),xlab="t",ylab="n/N")
points(seq(0,T),n[1,]/N,type="l",lwd=2)
points(seq(0,T),n[2,]/N,type="l",lwd=2,col="red")
points(seq(0,T),n[3,]/N,type="l",lwd=2,col="blue")

close.screen(all=TRUE)

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

# Impact of initial conditions

NoSim <- 200

n0 <- array(NA,c(3,NoSim+3))
for (j in 1:NoSim){
  n0[,j] <- runif(3,0,10)
  N0 <- sum(n0[,j])
  for (i in 1:3){n0[i,j] <- n0[i,j]/N0}
  }
n0[,NoSim+1] <- c(1,0,0)
n0[,NoSim+2] <- c(0,1,0)
n0[,NoSim+3] <- c(0,0,1)

result <- array(NA,c(3,T+1,NoSim+3))
N <- array(NA,c(T+1,NoSim+3))

for (j in (1:(NoSim+3))){
   result[,1,j] <- n0[,j]
   N[1,j] <- 1
   for (i in 2:(T+1)){
     result[,i,j] <- A %*% result[,i-1,j]
     N[i,j] <- sum(result[,i,j])}
   }

x11()
par(mfcol=c(1,2),mfg=c(1,1,1,2))
plot(c(0,T),c(min(N),max(N)),type="n",bty="l",xlab="t",ylab="N",
     log="y")
for (j in 1:NoSim){points(seq(0,T),N[,j],type="l")}
points(seq(0,T),N[,NoSim+1],type="l",col="green",lwd=2)
points(seq(0,T),N[,NoSim+2],type="l",col="violet",lwd=2)
points(seq(0,T),N[,NoSim+3],type="l",col="cyan",lwd=2)

par(mfcol=c(3,2),mfg=c(1,1,3,2))
ylabP <- c(expression(n[1]/N),expression(n[2]/N),expression(n[3]/N))
mfgP <- array(NA,c(3,2))
mfgP[1,] <- c(1,2)
mfgP[2,] <- c(2,2)
mfgP[3,] <- c(3,2)
for (k in 1:3){
   par(mfg=mfgP[k,])
   plot(c(0,T),c(0,1),type="n",bty="l",xlab="t",ylab=ylabP[k])
   for (j in 1:NoSim){
     points(seq(0,T),result[k,,j]/N[,j],type="l")}
   points(seq(0,T),result[k,,NoSim+1]/N[,NoSim+1],type="l",col="green")
   points(seq(0,T),result[k,,NoSim+2]/N[,NoSim+2],type="l",col="violet")
   points(seq(0,T),result[k,,NoSim+3]/N[,NoSim+3],type="l",col="cyan")
   }
close.screen(all=TRUE)

################################################################################
# Impact of matrix entries variation

T <- 80    # Time of projection


n0 <- matrix(c(1,0,0),c(3,1))

result <- array(NA,c(3,T+1))
N <- array(NA,T+1)

result[,1] <- n0
for (i in 2:(T+1)){result[,i] <- A %*% result[,i-1]}
for (i in 1:(T+1)){N[i] <- sum(result[,i])}


x11()
plot(seq(0,T),N,type="l",bty="l",xlab="t",ylab="N",ylim=c(0,max(N)),
     xlim=c(0,1.1*T))

Ap <- A
Ap[1,2] <- 0.9*Ap[1,2]
for (i in 2:(T+1)){result[,i] <- Ap %*% result[,i-1]}
for (i in 1:(T+1)){N[i] <- sum(result[,i])}
points(seq(0,T),N,type="l",col="red")
text(1.03*T,N[T+1]*1.05,expression(F[2]==0.9),col="red")

Ap <- A
Ap[1,3] <- 0.9*Ap[1,3]
for (i in 2:(T+1)){result[,i] <- Ap %*% result[,i-1]}
for (i in 1:(T+1)){N[i] <- sum(result[,i])}
points(seq(0,T),N,type="l",col="orange")
text(1.03*T,N[T+1]*1.2,expression(F[3]==4.5),col="orange")

Ap <- A
Ap[2,1] <- 0.9*Ap[2,1]
for (i in 2:(T+1)){result[,i] <- Ap %*% result[,i-1]}
for (i in 1:(T+1)){N[i] <- sum(result[,i])}
points(seq(0,T),N,type="l",col="blue")
text(1.03*T,N[T+1]*0.65,expression(P[1]==0.27),col="blue")

Ap <- A
Ap[3,2] <- 0.9*Ap[3,2]
for (i in 2:(T+1)){result[,i] <- Ap %*% result[,i-1]}
for (i in 1:(T+1)){N[i] <- sum(result[,i])}
points(seq(0,T),N,type="l",col="green")
text(1.03*T,N[T+1]*0.8,expression(P[2]==0.45),col="green")


################################################################################
# External variability

T <- 50

h <- runif(T,-10,10)
for (i in 1:(T)){
  if (h[i] > 0){h[i] <- 2}
    else{if (h[i] < 1){h[i] <- 0.5}
           else{h[i] <- 1}}
  }
result <- array(NA,c(3,T+1))
N <- array(NA,T+1)
result[,1] <- n0
N[1] <- sum(n0)
for (i in 1:T){
   Ap <- A
   Ap[1,2] <- A[1,2]*h[i]
   Ap[1,3] <- A[1,3]*h[i]
   result[,i+1] <- Ap %*% result[,i]
   N[i+1] <- sum(result[,i+1])
   }


x11()
par(mfcol=c(1,2),mfg=c(1,1,1,2))
plot(seq(0,T),N,type="l",bty="l",xlab="t",ylab="N",log="y",
     xlim=c(0,1.03*T))

par(mfg=c(1,2,1,2))
plot(seq(0,T),result[1,]/N,type="l",bty="l",xlab="t",ylab="Proportions",
     ylim=c(0,1),xlim=c(0,1.03*T))
points(seq(0,T),result[2,]/N,type="l",col="blue")
points(seq(0,T),result[3,]/N,type="l",col="red")
text(1.035*T,mean(result[1,]/N),expression(n[1]/N))
text(1.035*T,mean(result[2,]/N),expression(n[2]/N),col="blue")
text(1.035*T,mean(result[3,]/N),expression(n[3]/N),col="red")

close.screen(all=TRUE)

### Multiple simulation

NoSim <- 150

result <- array(NA,c(3,T+1,NoSim))
N <- array(NA,c(T+1,NoSim))

for (j in 1:NoSim){
  h <- runif(T,-100,100)
  for (i in 1:(T)){if (h[i] > 0){h[i] <- 2}
                             else{if (h[i] < 1){h[i] <- 0.5}
                                    else{h[i] <- 1}}}
  result[,1,j] <- n0
  N[1,j] <- sum(n0)
  for (i in 1:T){
    Ap <- A
    Ap[1,2] <- A[1,2]*h[i]
    Ap[1,3] <- A[1,3]*h[i]
    result[,i+1,j] <- Ap %*% result[,i,j]
    N[i+1,j] <- sum(result[,i+1,j])
    }
  }


x11()
par(mfcol=c(3,1),mfg=c(1,1,3,1))
plot(c(0,T),c(min(N),max(N)),type="n",bty="l",xlab="t",ylab="N",
     log="y",xlim=c(0,T*1.05))
for (j in 1:NoSim){points(seq(0,T),N[,j],type="l")}

par(mfg=c(2,1,3,1))
M <- array(NA,T+1)
V <- array(NA,T+1)
for (i in 1:(T+1)){M[i] <- mean(N[i,])
                           V[i] <- var(N[i,])}
plot(c(0,T),c(0.01,max(c(M,V))),log="y",type="n",bty="l",
     xlab="Time",ylab="Mean and Variance of N",xlim=c(0,T*1.05))
points(seq(0,T),M,type="l",col="violet",lwd=3)
points(seq(0,T),V,type="l",col="cyan",lwd=3)
text(T*1.005,V[T+1],"variance",adj=c(0,NULL))
text(T*1.005,M[T+1],"mean",adj=c(0,NULL))

P <- array(NA,c(T+1,3))
P10 <- array(NA,c(T+1,3))
P90 <- array(NA,c(T+1,3))
for (i in 1:(T+1)){
  for (k in 1:3){
    P[i,k] <- mean(result[k,i,]/N[i,])
    Q <- quantile(result[k,i,]/N[i,],c(0.1,0.9),names=FALSE)
    P10[i,k] <- Q[1]
    P90[i,k] <- Q[2]
    }
  }

par(mfcol=c(3,3),mfg=c(3,1,3,3))
plot(c(0,T),c(0,1),type="n",bty="l",xlab="t",
     ylab=expression(n[1]/N))
points(seq(0,T),P[,1],type="l",lwd=2)
points(seq(0,T),P10[,1],type="l")
points(seq(0,T),P90[,1],type="l")

par(mfg=c(3,2,3,3))
plot(c(0,T),c(0,1),type="n",bty="l",xlab="t",
     ylab=expression(n[2]/N))
points(seq(0,T),P[,2],type="l",col="blue",lwd=2)
points(seq(0,T),P10[,2],type="l",col="blue")
points(seq(0,T),P90[,2],type="l",col="blue")

par(mfg=c(3,3,3,3))
plot(c(0,T),c(0,1),type="n",bty="l",xlab="t",
     ylab=expression(n[3]/N))
points(seq(0,T),P[,3],type="l",col="red",lwd=2)
points(seq(0,T),P10[,3],type="l",col="red")
points(seq(0,T),P90[,3],type="l",col="red")

close.screen(all=TRUE)

################################################################################
# Internal variability

T <- 300

result <- array(NA,c(3,T+1))
N <- array(NA,T+1)

g <- function(N,b,R){R*exp(-b*N)}


R <- 1
b <- 0.005
result[,1] <- n0
N[1] <- sum(result[,1])
for (i in 1:T){
  An <- A
  An[1,2] <- A[1,2]*g(N[i],b,R)
  An[1,3] <- A[1,3]*g(N[i],b,R)
  result[,i+1] <- An %*% result[,i]
  N[i+1] <- sum(result[,i+1])}


x11()
par(mfcol=c(2,2),mfg=c(1,1,2,2))
plot(seq(0,T),result[1,],type="l",bty="l",xlab="t",
     ylab=expression(list(n[1],n[2],n[3])))
points(seq(0,T),result[2,],type="l",col="blue")
points(seq(0,T),result[3,],type="l",col="red")
text(T,result[1,T+1]-0.2,expression(n[1]))
text(T,result[2,T+1]-0.2,expression(n[2]),col="blue")
text(T,result[3,T+1]-0.2,expression(n[3]),col="red")

par(mfg=c(1,2,2,2))
plot(seq(0,T),N,type="l",bty="l",xlab="t",ylab="N")

NoSim <- 40

par(mfcol=c(2,1),mfg=c(2,1,2,1))
plot(seq(0,T),N,type="l",bty="l",xlab="t",ylab="N",ylim=c(0.2,50),
     log="y")

for (j in 1:floor(NoSim/2)){
  result[,1] <- runif(3,0,1)
  N[1] <- sum(result[,1])
  for (i in 1:T){
    An <- A
    An[1,2] <- A[1,2]*g(N[i],b,R)
    An[1,3] <- A[1,3]*g(N[i],b,R)
    result[,i+1] <- An %*% result[,i]
    N[i+1] <- sum(result[,i+1])}
  points(seq(0,T),N,type="l",col="green")
  }
for (j in 1:floor(NoSim/2)){
  result[,1] <- runif(3,0,5)
  N[1] <- sum(result[,1])
  for (i in 1:T){
    An <- A
    An[1,2] <- A[1,2]*g(N[i],b,R)
    An[1,3] <- A[1,3]*g(N[i],b,R)
    result[,i+1] <- An %*% result[,i]
    N[i+1] <- sum(result[,i+1])}
  points(seq(0,T),N,type="l",col="green")
}
close.screen(all=TRUE)
######################

T <- 80

result <- array(NA,c(3,T+1))
N <- array(NA,T+1)

R <- 45
B <- 0.005
result[,1] <- n0
N[1] <- sum(result[,1])
for (i in 1:T){
  An <- A
  An[1,2] <- A[1,2]*g(N[i],b,R)
  An[1,3] <- A[1,3]*g(N[i],b,R)
  result[,i+1] <- An %*% result[,i]
  N[i+1] <- sum(result[,i+1])}


x11()
plot(seq(0,T),N,type="l",bty="l",xlab="t",ylab="N",lwd=2)
close.screen(all=TRUE)