#########################################
###### Exemplo Juvenis e Adultos ########
#########################################
rm(list=ls())

#Tamanho inicial da Populacao
J0<-0
A0<-1

#Parametros
b<-10
p<-0.1
s<-0.2

matriz.JA=matrix(c(0, b, p, s), nrow=2, ncol=2, byrow=TRUE) #matriz com os paramentros
matriz.JA

#Calculando os autovalores e respectivos autovetores da matriz
eigen(matriz.JA)
eigenva.JA=eigen(matriz.JA)$values[1] # pegando o maior autovalor em modulo
eigenvc.JA=eigen(matriz.JA)$vectors[,1]  # pegando o espectivo autovetor do maior autovalor em modulo

#Proporcao de individuos:
#Juvenis:
prop.eigenJ=eigenvc.JA[1]/(eigenvc.JA[1]+eigenvc.JA[2])
prop.eigenJ
#Adultos:
prop.eigenA=eigenvc.JA[2]/(eigenvc.JA[1]+eigenvc.JA[2])
prop.eigenA

#Relacao A/J:
relAJ.eigen=eigenvc.JA[2]/eigenvc.JA[1]
relAJ.eigen


########################################################### 
#### Calculando o numero de individuos iterativamente: ####
###########################################################

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

#Definicao dos vetores de cada categoria e a populacao total para cada instante de tempo
J=rep(0,(T.max+1))
A=rep(0,(T.max+1))
Ntotal=J+A

#Primeiro valor dos vetores de categorias
J[1]<-J0
A[1]<-A0
Ntotal[1]<-J[1]+A[1]

#Calculando para cada instante de tempo (iterativamente):
for (n in 1:T.max){
  J[n+1]<-A[n]*b
  A[n+1]<-J[n]*p+A[n]*s  
  Ntotal[n+1]<-A[n+1]+J[n+1]
}

#Juntando os dados
tempo=seq(0,(T.max))
dados=data.frame(tempo,J,A, Ntotal)
head(dados)
tail(dados)

#Grafico do numero de indivíduos variando no tempo
plot(dados$tempo, dados$J, type="b", pch=19,ylim=c(0,100),xlim=c(0,50),
     col="red", xlab="tempo", ylab="Numero de Individuos") 
lines(tempo, A, pch=18, col="blue", type="b", lty=2)
legend(1, 100, legend=c("Juvenis", "Adultos"),col=c("red", "blue"), lty=1:2, cex=0.8)

plot(dados$tempo, dados$J, type="b", pch=19,ylim=c(0,10000),xlim=c(0,100),
     col="red", xlab="tempo", ylab="Numero de Individuos") 
lines(dados$tempo, dados$A, pch=18, col="blue", type="b", lty=2)
legend(1, 10000, legend=c("Juvenis", "Adultos"),col=c("red", "blue"), lty=1:2, cex=0.8)


############################################# 
#### Calculando crescimento da populacao ####
#############################################

# Definicao dos vetores para gaurdar o crescimento a cada
# instante de tempo em cada categoria e a populacao total:
crescimentoJ=rep(0,(T.max-1)) 
crescimentoA=rep(0,(T.max-1))
crescimento=rep(0,(T.max-1))

for (n in 1:(T.max)){ #calculando iterativamente o tamanho em um tempo seguinte em relacao ao tempo anterior
  crescimentoJ[n]<-J[n+1]/J[n]
  crescimentoA[n]<-A[n+1]/A[n]
  crescimento[n]<-Ntotal[n+1]/Ntotal[n]
}

tempo.c=seq(0,(T.max-1)) #definicao do tempo (decresce de 1 agora)
dados.c=data.frame(tempo.c,crescimentoJ,crescimentoA, crescimento) #juntando os dados
head(dados.c)
dados.c[35:55,]
tail(dados.c)

eigenva.JA #calculado pelo maior autovalor em modulo la no inicio

#Grafico do Crescimento da populacao pelo tempo
plot(dados.c$crescimentoJ ~ dados.c$tempo.c, type="b", pch=19,  ylim=c(0,5),
     xlim=c(0,100),col="red", xlab="Tempo", ylab="Taxa de aumento da Populaçao")
lines(dados.c$crescimentoA ~ dados.c$tempo.c, pch=18, col="blue", type="b", lty=2)
legend(80, 5, legend=c("Juvenis", "Adultos"),col=c("red", "blue"), lty=1:2, cex=0.8)
abline(h=eigenva.JA)

# Checando o valor de crescimento da populacao pela regressao linear nos tempos finais (quase estavel)
regressao.cresc=lm(crescimento[50:T.max]~tempo.c[50:(T.max)]) 
regressao.cresc
dados.c[100,3] #comparando com os dados de crescimento calculados acima
eigenva.JA  #comparando com o maior autovalor em modulo la no inicio


#################################################
#### Calculando proporcao A e J na populacao ####
#################################################
# Atribuindo as proporcoes baseado na populacao de 
# J e A dividido pelo Ntotal a cada instante de tempo:
proporcaoJ=J/Ntotal 
tail(proporcaoJ)
prop.eigenJ #checando com o valor calculado a partir do autovetor
proporcaoA=A/Ntotal
tail(proporcaoA)
prop.eigenA  #checando com o valor calculado a partir do autovetor

#Grafico da Porporcao de indivíduos variando no tempo
plot(tempo, proporcaoJ, type="b", pch=19, xlim=c(0,100), ylim=c(0,1),col="red",  xlab="tempo", ylab="Proporção de Individuos (%)")
abline(h=prop.eigenJ,col="red4", lwd=2, lty=2)
lines(tempo,proporcaoA, pch=18, col="blue", type="b", lty=2)
abline(h=prop.eigenA, col="blue4",lwd=2, lty=2)
legend(80, 0.8, legend=c("Juvenis", "Adultos"), col=c("red", "blue"), lty=1:2, cex=0.8)

#Calculando pela regressao
regressao.propJ=lm(proporcaoJ[50:T.max]~tempo.c[50:(T.max)]) 
regressao.propJ
prop.eigenJ  #comparando com o valor calculado a partir do autovetor

#Calculando pela regressao
regressao.propA=lm(proporcaoA[50:T.max]~tempo.c[50:(T.max)]) 
regressao.propA
prop.eigenA  #comparando com o valor calculado a partir do autovetor


############################################# 
#### Calculando relacao A/J na populacao ####
#############################################

# Calculando as relacoes A/J a cada instante de tempo:
relacaoAJ=A/J
head(relacaoAJ)
relacaoAJ[35:55]
tail(relacaoAJ)
relAJ.eigen #checando com o valor calculado a partir do autovetor

#Grafico de A x J. O coeficiente algular deve ser igual a relacao A/J
plot(A~J , type="b", pch=19,  ylim=c(0,1100),xlim=c(0,10000),xlab="Juvenis",ylab="Adultos")
abline(1,relAJ.eigen, col="red", lwd=2)
abline(1, relacaoAJ[100], lty=2, lwd=2)

#Zoom:
plot(A~J , type="b", pch=19,  ylim=c(0,12),xlim=c(0,100),xlab="Juvenis",ylab="Adultos")
abline(1,relAJ.eigen, col="red",lwd=2)
abline(1, relacaoAJ[100], lty=2, lwd=2)

#Calculando pela regressao
regressao.relAJ=lm(A[80:T.max]~J[80:T.max])
regressao.relAJ
relacaoAJ[100] # comparando com dados de relacaoAJ calculados acima
relAJ.eigen # comparando com o valor calculado a partir do autovetor