###################################
###### Exemplo das baleias ########
###################################
rm(list=ls())

# Tamanho inicial da Populacao
nf.0<-0
ni.0<-4
nm.0<-0
nr.0<-0

# Parametros
s.if<-0.92
s.ii<-0.86
s.mi<-0.08
s.ri<-0.02
s.mm<-0.80
s.mr<-0.88
s.rm<-0.19
s.rr<-0
b<-0.3

# Tempo total a ser analisado
T.max<-200

# Definicao dos vetores de cada categoria e pop. total para cada instante de tempo
NF=matrix(0,T.max,1)
NI=matrix(0,T.max,1)
NM=matrix(0,T.max,1)
NR=matrix(0,T.max,1)
Ntotal=NF+NI+NM+NR


# Primeiro valor dos vetores de categorias
NF[1]<-nf.0
NI[1]<-ni.0
NM[1]<-nm.0
NR[1]<-nr.0
Ntotal[1]<-NF[1]+NI[1]+NM[1]+NR[1]

# Calculando o tamanho de cada classe para cada instante de tempo (iterativamente):
for (t in 1:T.max){
  NF[t+1]<-NR[t]*b
  NI[t+1]<-NI[t]*s.ii+NF[t]*s.if
  NM[t+1]<-NM[t]*s.mm+NI[t]*s.mi+NR[t]*s.mr
  NR[t+1]<-NR[t]*s.rr+NI[t]*s.ri+NM[t]*s.rm
  Ntotal[t+1]<-NF[t+1]+NI[t+1]+NM[t+1]+NR[t+1]
}

# Juntando os dados
tempo<-seq(0,T.max)
dados=data.frame(tempo,NF,NI,NM,NR, Ntotal=Ntotal)
head(dados)
tail(dados)

# Plot numero de indivíduos variando no tempo
pdf("/Users/flaviamarquitti/Documents/LivroModEcoEvo/templateBonito/fig1Baleias.pdf", height=6, width=6)
plot(dados$tempo, dados$NF, type="o", pch=15, xlim=c(0,50),  ylim=c(0,5),col="red", 
     xlab="tempo", ylab="Número de Individuos", lty=2, cex.lab=1.5,cex=0.8) 
lines(dados$tempo, dados$NI, pch=18, col="blue", type="o", lty=2,cex=0.8)
lines(dados$tempo, dados$NM, pch=17, col="grey", type="o", lty=2,cex=0.8)
lines(dados$tempo, dados$NR, pch=16, col="orange", type="o", lty=2,cex=0.8)
legend(35, 5, legend=c("Filhotes", "Imaturos", "Maduros", "Reprodutivos"),
       col=c("red", "blue", "grey", "orange"), lty=1, lwd=2, cex=0.8)
dev.off()

# Proporcao em cada unidade de tempo
proporcao= matrix(0, ncol=4, nrow=(T.max+1))
proporcao[,1]=dados$NF/dados$Ntotal
proporcao[,2]=dados$NI/dados$Ntotal
proporcao[,3]=dados$NM/dados$Ntotal
proporcao[,4]=dados$NR/dados$Ntotal
tempo<-seq(0,T.max)

proporcao.1=cbind(tempo,proporcao)
proporcao.df=as.data.frame(proporcao.1)
colnames(proporcao.df)=c("tempo","NF","NI","NM","NR")

tail(proporcao.df)

#Plot Porporcao de indivíduos variando no tempo
?pdf()
pdf("/Users/flaviamarquitti/Documents/LivroModEcoEvo/templateBonito/fig2Baleias.pdf", height=6, width=6)
plot(proporcao.df$tempo, proporcao.df$NF, type="o", pch=15, xlim=c(0,50), ylim=c(0,1),col="red", 
     xlab="tempo", ylab="Proporção de Indivíduos (%)", cex.lab=1.5, cex=0.8, lty=2)
lines(proporcao.df$tempo, proporcao.df$NI, pch=18, col="blue", type="o", lty=2,cex=0.8)
lines(proporcao.df$tempo, proporcao.df$NM, pch=17, col="grey", type="o", lty=2,cex=0.8)
lines(proporcao.df$tempo, proporcao.df$NR, pch=16, col="orange", type="o", lty=2,cex=0.8)
legend(35, 1, legend=c("Filhotes", "Imaturos", "Maduros", "Reprodutivos"),
       col=c("red", "blue", "grey", "orange"), lty=1, lwd=2, cex=0.8)
dev.off()


# Calculando o crescimento da populacao total
crescimento=matrix(0,nrow=(T.max),ncol=5)

# Calculando o crescimento iterativamente:
for (i in 1:(T.max)){
  crescimento[i,1]=dados$NF[i+1]/dados$NF[i]
  crescimento[i,2]=dados$NI[i+1]/dados$NI[i]
  crescimento[i,3]=dados$NM[i+1]/dados$NM[i]
  crescimento[i,4]=dados$NR[i+1]/dados$NR[i]
  crescimento[i,5]<-dados$Ntotal[i+1]/dados$Ntotal[i]
}

# Compilando os dados de crescimento:
tempo<-seq(0,(T.max-1))
crescimento.1=cbind(tempo,crescimento)
crescimento=as.data.frame(crescimento.1)
colnames(crescimento)=c("tempo","NF","NI","NM","NR", "Ntotal")


# Plot Crescimento da populacao
plot(crescimento$tempo, crescimento$Ntotal, type="b", pch=18,  ylim=c(0.96,1.005),xlim=c(0,100), 
     xlab="Tempo", ylab="Taxa de aumento da Populaçao")

# Plot do Crescimento de cada classe de indivíduos variando no tempo
pdf("/Users/flaviamarquitti/Documents/LivroModEcoEvo/templateBonito/fig3Baleias.pdf", height=6, width=6)
par(bg="gray")
plot(crescimento$tempo, crescimento$NF, type="o", pch=15, xlim=c(0,50), ylim=c(0,2),col="red", 
     xlab="tempo", ylab="Taxa de aumento da Populaçao", cex=0.8, cex.lab=1.5, lty=2)
lines(crescimento$tempo, crescimento$NI, pch=18, col="blue", type="o", lty=2,cex=0.8)
lines(crescimento$tempo, crescimento$NM, pch=17, col="grey", type="o", lty=2,cex=0.8)
lines(crescimento$tempo, crescimento$NR, pch=16, col="orange", type="o", lty=2,cex=0.8)
lines(crescimento$tempo, crescimento$Ntotal, pch=19, col="black", type="o", lty=1, cex=0.3)
legend(35, 2, legend=c("Filhotes", "Imaturos", "Maduros", "Reprodutivos", "Total"),
       col=c("red", "blue", "grey", "orange", "black"), lty=1, lwd=2,cex=0.8)
dev.off()

# Matriz de projeção usada no exemplo acima:
Mat = matrix(c(0.0, 0.0, 0.0, b,
               s.if,s.ii,0.0,0.0,
               0.0,s.mi,s.mm,s.mr,
               0.0,s.ri,s.rm,s.rr),nrow=4,ncol=4,byrow = TRUE)
Mat

# Calculando os autovalores e autovetores da matriz de transição Mat
EV=eigen(Mat)
EV

# Autovalores
EV$values

# Compare o maior autovalor em valor absoluto com o crescimento calculado iterativamente:
tail(crescimento) # pelos dados calculados na mao

# Autovetores associados aos autovalores (os autovetores estao por coluna)
EV$vectors

# Calculando a razão do número filhotes, imaturos, e reprodutivos em relação ao número
# maduros através das entradas do autovetor associado ao maior autovalor (neste caso eh o primeiro):
razaoF=EV$vectors[1,1]/EV$vectors[3,1]
razaoI=EV$vectors[2,1]/EV$vectors[3,1]
razaoM=EV$vectors[3,1]/EV$vectors[3,1]
razaoR=EV$vectors[4,1]/EV$vectors[3,1]

razaoF
razaoI
razaoM
razaoR

# Calculando a proporcao (porcentagem) do número filhotes, imaturos, e reprodutivos em relação ao número ao total 
Propocao.EV.F=EV$vectors[1,1]/(EV$vectors[1,1]+EV$vectors[2,1]+EV$vectors[3,1]+EV$vectors[4,1])
Propocao.EV.I=EV$vectors[2,1]/(EV$vectors[1,1]+EV$vectors[2,1]+EV$vectors[3,1]+EV$vectors[4,1])
Propocao.EV.M=EV$vectors[3,1]/(EV$vectors[1,1]+EV$vectors[2,1]+EV$vectors[3,1]+EV$vectors[4,1])
Propocao.EV.R=EV$vectors[4,1]/(EV$vectors[1,1]+EV$vectors[2,1]+EV$vectors[3,1]+EV$vectors[4,1])

Propocao.EV.F
Propocao.EV.I
Propocao.EV.M
Propocao.EV.R

# Compare com a proporcao de individuos em cada classe calculado iterativamente:
tail(proporcao.df)




#########################################################################################
###### Calculo Autovalores e Autovetores da tabela de vida de insetos hipoteticos #######
#########################################################################################
rm(list=ls())

# Matriz de transição de uma populaçao hipotetica de insetos montada a partir de dados da tabela de vida
M = matrix(c(1.6, 1.5, 0.25, 0.0,
             0.8, 0.0, 0.0 , 0.0,
             0.0, 0.5, 0.0 , 0.0,
             0.0, 0.0, 0.25, 0.0), nrow=4, ncol=4, byrow = TRUE)
M

# Calculo dos autovalores
EV=eigen(M)
EV

# Obtendo o maximo dos autovalores em valor absoluto:
max(abs(EV$values))

#Calculo dos estrutura etaria estavel
stable.age=EV$vectors[,which(abs(EV$values)==max(abs(EV$values)))]/sum(EV$vectors[,which(abs(EV$values)==max(abs(EV$values)))]) 

stable.age

#tempo total a ser analisado
T.max<-300

#definicao da matriz com cada categoria/classe nas colunas
x=matrix(0,T.max,4)

#x[1,]=c(500,0,0,0) # numero de individuos em cada classe no instante inicial 1

x[1,]=c(0,40,20,10) # numero de individuos em cada classe no instante inicial 2

# Calculando iterativamente:
for(n in 1:(T.max-1)){
  x[n+1,] = M %*% x[n,] # multiplicacao matricial
  
}

# Montando a tabela com os dados:
tempo<-seq(0,(T.max-1))
dados=cbind(tempo,x)
dados=as.data.frame(dados)
colnames(dados)=c("tempo", "x.0","x.1","x.2","x.3")

# Plot Numero de individuos em cada classe etaria variando no tempo:
plot(dados$tempo, dados$x.0, type="b", pch=19, xlim=c(0,30), ylim=c(0,100000),col="red", 
     xlab="tempo", ylab="Tamanho de cada classe da Populacao")
lines(dados$tempo, dados$x.1, pch=18, col="blue", type="b", lty=2)
lines(dados$tempo, dados$x.2, pch=17, col="grey", type="b", lty=2)
lines(dados$tempo, dados$x.3, pch=16, col="orange", type="b", lty=2)
legend(20, 100000, legend=c("x0", "x1", "x2", "x3"),
       col=c("red", "blue", "grey", "orange"), lty=1:2, cex=0.8)


# Plot Log Numero de individuos em cada classe etaria variando no tempo:
plot(dados$tempo, log(dados$x.0), type="b", pch=19, xlim=c(0,50), ylim=c(0,60),col="red", 
     xlab="tempo", ylab="Log(tamanho de cada classe da Populacao)")
lines(dados$tempo, log(dados$x.1), pch=18, col="blue", type="b", lty=2)
lines(dados$tempo, log(dados$x.2), pch=17, col="grey", type="b", lty=2)
lines(dados$tempo, log(dados$x.3), pch=16, col="orange", type="b", lty=2)
legend(40, 60, legend=c("x0", "x1", "x2", "x3"),
       col=c("red", "blue", "grey", "orange"), lty=1:2, cex=0.8)


# Mostrando que as retas do grafico acima estao parealelas, portanto a pop cresce a uma mesma taxa
# em todas as classes. Note que o coeficiente angular dessas retas eh exatamente o Ln do 
# maior autovalor em modulo encontrado pela matriz de transicao "M" acima. Note que independente
# da condicao inicial, para um tempo longo, este coeficiente angular eh o mesmo (assim como eh o autovalor dominante)
lm(log(dados$x.0[100:300])~dados$tempo[100:300]) 
lm(log(dados$x.1[100:300])~dados$tempo[100:300]) 
lm(log(dados$x.2[100:300])~dados$tempo[100:300]) 
lm(log(dados$x.3[100:300])~dados$tempo[100:300]) 


# Curva de sobrevivencia da populaçao hipotetica de insetos
l_i=c(1,0.8,0.4,0.1,0)
i=c(0,1,2,3,4)
plot(log(lx)~x, pch=19, xlab="idade (i)", ylab="ln(l_i)", type="b", lty=2)



# Fazendo a estimativa de crescimento da pop pela tabela de vida (usada em aula):
tabela.de.vida=matrix(c(0, 500, 0, 1.00, NaN,  0.80,
                        1, 400, 2, 0.80, 0.80, 0.50,
                        2, 200, 3, 0.40, 0.50, 0.25,
                        3, 50 , 1, 0.10, 0.25, 0.00,
                        4, 0  , 0, 0.00, 0.00, NaN), nrow=5, ncol=6, byrow=T)

tabela.de.vida=as.data.frame(tabela.de.vida)
colnames(tabela.de.vida)=c("i", "S.i", "b.i", "l.i", "p.i", "sigma.i")

# Veja a tabela de vida:
tabela.de.vida

# Calculando o valor de R0:
R.0=sum(tabela.de.vida$l.i*tabela.de.vida$b.i)
R.0

# Calculando o valor de tempo de geracao:
G=sum(tabela.de.vida$i*tabela.de.vida$l.i*tabela.de.vida$b.i)/R.0
G

# Usando a forma da eq N.G=N.0*(Lambda^G) para estimar o valor de Lambda (ou R). N.G/N.0 eh 
# aproximadamente R.0. Desta forma, pela eq, temos que log(R.0) = log(Lambda^G) = G*log(Lambda).
# Portando Lambda = exp(log(lambda))
Log.Lambda=log(R.0)/G
Log.Lambda

lambda=exp(Log.Lambda)
lambda

# Compare o lambda acima com o maior autovalor encontrado pela matriz de trasicao: 
max(abs(EV$values))