###############################################
###### Exemplo de mutualismo -  linear ########
###############################################
rm(list=ls())
Mut=matrix(c(1/2,  5/6,
             12/5, 3/2), ncol=2,nrow=2, byrow=T)


EV.mut=eigen(Mut)
max(abs(EV.mut$values[1]))
razao.EV.Sp1=EV.mut$vectors[1,1]/EV.mut$vectors[2,1]
razao.EV.Sp2=EV.mut$vectors[2,1]/EV.mut$vectors[2,1]
razao.EV.Sp1
razao.EV.Sp2

#########

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

#definicao dos vetores de cada categoria e pop. total para cada instante de tempo
Sp1=matrix(0,T.max,1)
Sp2=matrix(0,T.max,1)

#Primeiro valor dos vetores das spp
Sp1[1]<-4
Sp2[1]<-2

#Calculando para cada  instante de tempo
for (n in 1:T.max){
  Sp1[n+1]<-Mut[1,1]*Sp1[n]+Mut[1,2]*Sp2[n]
  Sp2[n+1]<-Mut[2,1]*Sp1[n]+Mut[2,2]*Sp2[n]
}

#Juntando os dados
tempo<-seq(0,T.max)
dados=data.frame(tempo,Sp1,Sp2)
head(dados)
tail(dados)
razao=dados
razao$Sp1=dados$Sp1/dados$Sp2
razao$Sp2=dados$Sp2/dados$Sp2

tail(razao)
razao.EV.Sp1
razao.EV.Sp2

#Plot numero de indivíduos variando no tempo
plot(dados$tempo, dados$Sp1, type="b", pch=19,xlim=c(0,10), ylim=c(0,8000),col="red", 
     xlab="tempo", ylab="Numero de Individuos") 
lines(dados$tempo, dados$Sp2, pch=18, col="blue", type="b", lty=2)
legend(5, 7000, legend=c("Sp1", "Sp2"),
       col=c("red", "blue"), lty=1:2, cex=0.8)



###################################################
###### Exemplo de predacao fraca -  linear ########
##################################################
rm(list=ls())
Pred=matrix(c(0.5,  0.25,
             -0.5, 1.5), ncol=2,nrow=2, byrow=T)


EV.pred=eigen(Pred)
EV.pred
max(abs(EV.pred$values[1]))
razao.EV.Sp1=EV.pred$vectors[1,1]/EV.pred$vectors[2,1]
razao.EV.Sp2=EV.pred$vectors[2,1]/EV.pred$vectors[2,1]
razao.EV.Sp1
razao.EV.Sp2

#########

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

#definicao dos vetores de cada categoria e pop. total para cada instante de tempo
Sp1=matrix(0,T.max,1)
Sp2=matrix(0,T.max,1)

#Primeiro valor dos vetores das spp
Sp1[1]<-2
Sp2[1]<-4

#Calculando para cada  instante de tempo
for (n in 1:T.max){
  Sp1[n+1]<-Pred[1,1]*Sp1[n]+Pred[1,2]*Sp2[n]
  Sp2[n+1]<-Pred[2,1]*Sp1[n]+Pred[2,2]*Sp2[n]
}

#Juntando os dados
tempo<-seq(0,T.max)
dados=data.frame(tempo,Sp1,Sp2)
head(dados)
tail(dados)
razao=dados
razao$Sp1=dados$Sp1/dados$Sp2
razao$Sp2=dados$Sp2/dados$Sp2

tail(razao)
razao.EV.Sp1
razao.EV.Sp2

#Plot numero de indivíduos variando no tempo
plot(dados$tempo, dados$Sp1, type="b", pch=19,xlim=c(0,10), ylim=c(0,60),col="red", 
     xlab="tempo", ylab="Numero de Individuos") 
lines(dados$tempo, dados$Sp2, pch=18, col="blue", type="b", lty=2)
legend(5, 50, legend=c("Sp1", "Sp2"),
       col=c("red", "blue"), lty=1:2, cex=0.8)




###################################################
###### Exemplo de predacao forte -  linear ########
##################################################
rm(list=ls())
Pred.intensa=matrix(c(0.5,  1,
                     -2.5, 1.5), ncol=2,nrow=2, byrow=T)


EV.predint=eigen(Pred.intensa)
EV.predint
max(abs(EV.predint$values[1]))

#########


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

#definicao dos vetores de cada categoria e pop. total para cada instante de tempo
Sp1=matrix(0,T.max,1)
Sp2=matrix(0,T.max,1)

#Primeiro valor dos vetores das spp
Sp1[1]<-5
Sp2[1]<-1

#Calculando para cada  instante de tempo
for (n in 1:T.max){
  Sp1[n+1]<-Pred.intensa[1,1]*Sp1[n]+Pred.intensa[1,2]*Sp2[n]
  Sp2[n+1]<-Pred.intensa[2,1]*Sp1[n]+Pred.intensa[2,2]*Sp2[n]
}

#Juntando os dados
tempo<-seq(0,T.max)
dados=data.frame(tempo,Sp1,Sp2)
head(dados)
tail(dados)


#Plot numero de indivíduos variando no tempo
plot(dados$tempo, dados$Sp1, type="b", pch=19,xlim=c(0,10), ylim=c(-2000,2000),col="red", 
     xlab="tempo", ylab="Numero de Individuos") 
lines(dados$tempo, dados$Sp2, pch=18, col="blue", type="b", lty=2)
legend(5, 50, legend=c("Sp1", "Sp2"),
       col=c("red", "blue"), lty=1:2, cex=0.8)




###################################
###### 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 para cada  instante de tempo
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
plot(dados$tempo, dados$NF, type="b", pch=19,xlim=c(0,100),  ylim=c(0,5),col="red", 
     xlab="tempo", ylab="Numero de Individuos") 
lines(dados$tempo, dados$NI, pch=18, col="blue", type="b", lty=2)
lines(dados$tempo, dados$NM, pch=17, col="grey", type="b", lty=2)
lines(dados$tempo, dados$NR, pch=16, col="orange", type="b", lty=2)
legend(70, 5, legend=c("Filhotes", "Imaturos", "Maduros", "Reprodutivos"),
       col=c("red", "blue", "grey", "orange"), lty=1:2, cex=0.8)


# 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
plot(proporcao.df$tempo, proporcao.df$NF, type="b", pch=19, xlim=c(0,100), ylim=c(0,1),col="red", 
     xlab="tempo", ylab="Proporção de Individuos (%)")
lines(proporcao.df$tempo, proporcao.df$NI, pch=18, col="blue", type="b", lty=2)
lines(proporcao.df$tempo, proporcao.df$NM, pch=17, col="grey", type="b", lty=2)
lines(proporcao.df$tempo, proporcao.df$NR, pch=16, col="orange", type="b", lty=2)
legend(70, 1, legend=c("Filhotes", "Imaturos", "Maduros", "Reprodutivos"),
       col=c("red", "blue", "grey", "orange"), lty=1:2, cex=0.8)



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


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]
}


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 Crescimento de cada classe de indivíduos variando no tempo
plot(crescimento$tempo, crescimento$NF, type="b", pch=19, xlim=c(0,50), ylim=c(0,2),col="red", 
     xlab="tempo", ylab="Taxa de aumento da Populaçao")
lines(crescimento$tempo, crescimento$NI, pch=18, col="blue", type="b", lty=2)
lines(crescimento$tempo, crescimento$NM, pch=17, col="grey", type="b", lty=2)
lines(crescimento$tempo, crescimento$NR, pch=16, col="orange", type="b", lty=2)
legend(35, 2, legend=c("Filhotes", "Imaturos", "Maduros", "Reprodutivos"),
       col=c("red", "blue", "grey", "orange"), lty=1:2, cex=0.8)


#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)

#autovalores
tail(crescimento) # pelos dados calculados na mao

EV$values

#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
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

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

tail(proporcao) # pelos dados calculados na mao