====== Chuva de Propágulos ======
{{:roteiro:chuvaprop.jpg?150 | Ênio e Beto andando na chuva http://desenhos.kids.sapo.pt/a-chuva.htm}}

Neste exercício, iremos simular a dinâmica de uma metapopulação onde as probabilidades de extinção e de colonização de manchas são constantes (modelo de chuva de propágulos). Neste modelo há sempre uma fonte constante de imigrantes que podem colonizar qualquer mancha vazia. Relembrando, a variação da fração de manchas ocupadas no tempo é descrito pela seguinte fórmula geral:

$$\frac{df}{dt}=I - E $$

 onde $I$ é a taxa de entrada de migrantes, e $E$ a taxa de saída. A partir dele podemos definir um modelo simples para a dinâmica de ocupação de manchas que formam uma metapopulação:



$$\frac{df}{dt}=p_i(1 - f)-p_e f $$
 

onde $p_i$ é a taxa de imigração ou colonização, $ p_e$ é a taxa de extinção e $ f $ é a fração de manchas ocupadas (número de manchas ocupadas / número total de manchas).

Em primeiro lugar, vamos estabelecer a probabilidade de colonização de manchas vazias (pi), a probabilidade de extinção em manchas ocupadas (pe) e a fração inicial de manchas ocupadas (fi) como 30%, 15% e 40%, respectivamente.
<code>
pi=0.3
pe=0.15
fi=0.4
</code>

Agora que já temos os parâmetros do nosso modelo (pi, pe e fi) vamos criar nossa paisagem virtual. Essa paisagem deve ser constituída por manchas de habitat que podem estar ocupadas ou não. Para tanto, criaremos uma matriz de 10 linhas e 10 colunas, sendo que cada célula dessa matriz representará uma mancha. Serão, portanto, 100 manchas.
<code>
paisag=array(0,dim=c(10,10,1))
paisag
</code>

Muito bem, mas temos somente a matriz no tempo inicial. Para acompanhar a dinâmica de ocupação de manchas é preciso criarmos uma terceira dimensão: o tempo. Por enquanto vamos criar 10 passos no tempo.
<code>
paisag=array(0,dim=c(10,10,10))
paisag
</code>
Ótimo! Agora temos nossa paisagem em 10 tempos diferentes.

Finalmente, antes de começarmos a brincadeira, precisamos definir quais manchas estarão ocupadas no tempo inicial. Vamos combinar que quando a mancha está ocupada ela recebe o valor 1 e quando está vazia recebe o valor 0. Do jeito que a matriz está, repleta de zeros, nenhuma mancha está ocupada. Então vamos preencher ao acaso algumas manchas, usando um fi de 40%:
<code>
fi
1-fi
rep(1, fi*100)
ocor<-rep(1,fi*100)
nocor<-rep(0,(1-fi)*100)
nocor
c(ocor,nocor)
sample(c(ocor,nocor))
paisag[,,1]<-sample(c(ocor,nocor))
paisag[,,1]
image(paisag[,,1], col=c("white", "green"))
</code>
Agora sim, estamos prontos para a simulação. 

O passo seguinte é fazer as coisas acontecerem! Mas para isso temos que ter em mente todas as possibilidades:
  * **Manchas ocupadas**
    * podem permanecer ocupadas (1-pe = 0,85)
    * podem sofrer extinção local (pe = 0,15)
  * **Manchas vazias**
    * podem permanecer vazias (1-pi = 0,7)
    * podem ser colonizadas (pi = 0,3) 

Calma, não se assuste. O monstrinho abaixo vai tratar apenas das manchas que estavam ocupadas no tempo inicial.
<code>
sum(paisag[,,1]) # numero de manchas ocupadas no tempo 1
paisag[,,2][paisag[,,1]==1]<-sample(c(0,1),sum(paisag[,,1]),replace=T,prob=c(pe,1-pe))
paisag[,,2]
</code>

Se você observar atentamente, algumas manchas que estavam ocupadas (1) ficaram vazias (0):

<code>
par(mfrow=c(1,2))
image(paisag[,,1], col=c("white", "green"))
image(paisag[,,2], col=c("white", "green"))
</code>
Estes foram eventos de extinção local. Note que o número de manchas ocupadas diminuiu.
<code>
sum(paisag[,,1]) # total de manchas ocupadas no tempo inicial
sum(paisag[,,2]) # total de manchas ocupadas no tempo 2
</code>

Mas essa é apenas metade da história, não? Estamos esquecendo que as manchas vazias podem ser ocupadas. Vamos simular a colonização com um outro monstrinho:
<code>
length(paisag[,,1]) ## esse é o tamanho da nossa paisagem ou o número de manchas!
nmanchas=length(paisag[,,1]) ## vamos guardar esse valor que não muda durante toda a simulação
paisag[,,2][paisag[,,1]==0]<-sample(c(0,1),nmanchas-sum(paisag[,,1]),replace=T,prob=c(1-pi,pi))
paisag[,,2]
</code>

Compare a paisagem no tempo inicial e no tempo 2 e veja se no final o número de manchas ocupadas aumentou ou  diminuiu.
<code>
paisag[,,1] 
paisag[,,2]

sum(paisag[,,1]) # total de manchas ocupadas no tempo inicial
sum(paisag[,,2]) # total de manchas ocupadas no tempo 2
</code>
Note que algumas manchas que estavam ocupadas no tempo 1 ficaram vazias no tempo 2 e vice-versa. Note também que algumas manchas continuaram vazias e outras continuaram ocupadas. O próximo passo é calcular a fração de manchas ocupadas inicialmente (f1) e no tempo 2 (f2) e depois a diferença entre elas:
<code>
nmanchas
f1=sum(paisag[,,1])/nmanchas# fração das manchas ocupadas no tempo 1
f2=sum(paisag[,,2])/nmanchas # fração das manchas ocupadas no tempo 2

f1
f2
f2-f1
</code>
Eis a dinâmica de metapopulações! ;-) Mas até aqui vimos apenas a variação do f do tempo 1 para o tempo 2. Restam ainda 8 passos para chegarmos no tempo 10, só que para não termos que repetir todo o trabalho 8 vezes usaremos um comando chamado **for**, que fará o serviço por nós:
<code>
paisag=array(0,dim=c(10,10,10))
nmanchas=length(paisag[,,1])
paisag[,,1]=sample(c(rep(0,((1-fi)*nmanchas)),rep(1,fi*nmanchas)))
resultado=numeric()
for(t in 2:10)
               {
	       paisag[,,t][paisag[,,(t-1)]==1]<-sample(c(0,1),sum(paisag[,,(t-1)]) ,replace=T, prob=c(pe,1-pe))
	       paisag[,,t][paisag[,,(t-1)]==0]<-sample(c(0,1),nmanchas - sum(paisag[,,(t-1)]), replace=T,prob=c(1-pi,pi))
	       resultado[t-1]=sum(paisag[,,t])/nmanchas
	      }
resultado
</code>
Vôa lá! Agora temos as frações de manchas ocupadas (f) ao longo do tempo. Só para deixar nossos resultados mais bonitos, vamos colocá-los em uma tabela: 
<code>
resultado=data.frame(t=1:10,f=c(fi,resultado))
resultado
</code>

Agora, para ficar mais bonito ainda, vamos criar um gráfico com esses resultados:
<code>
par(mfrow=c(1,1))
plot(1:10,resultado$f,type="l",xlab="Tempo",ylab="Fração de manchas ocupadas",
ylim=c(0,1),main="Dinâmica de ocupação de manchas",font.lab=2,lwd=2)
</code>

Como nossa filosofia de vida é (ou pelo menos deveria ser) melhorar sempre, vamos acrescentar a esse gráfico uma informação muito importante: a fração de manchas ocupadas no equilíbrio (F). O F é calculado da seguinte forma:

$$F=\frac{p_i}{p_i+p_e}$$

Portanto...
<code>
F=pi/(pi+pe)
F
</code>
Agora que conhecemos o F vamos colocá-lo no gráfico na forma de uma linha horizontal:
<code>
plot(1:10,resultado$f,type="l",xlab="Tempo",ylab="Fração de manchas ocupadas",
ylim=c(0,1),main="Dinâmica de ocupação de manchas",font.lab=2,lwd=2)
abline(h=F,col=2,lwd=2,lty=2)
</code>

O serviço sujo está feito, agora é hora da diversão! Como somos bonzinhos, resolvemos criar uma função para facilitar a vida de vocês. Com essa função vocês poderão variar os parâmetros do nosso modelo à vontade, sem medo de ser feliz. Por favor, retribuam a gentileza e testem vários valores para cada parâmetro e vejam o que acontece com nossas metapopulações virtuais. Abaixo a função:
<code>
metapop=function(tf,c,l,fi,pi,pe){
	paisag=array(0,dim=c(l,c,tf))
        nmanchas=c*l
	paisag[,,1]=matrix(sample(c(1,0),nmanchas,prob=c(fi,1-fi), replace=T),l,c)
	resultado=numeric()
	for(t in 2:tf){
	       paisag[,,t][paisag[,,(t-1)]==1]<-sample(c(0,1),sum(paisag[,,(t-1)]), replace=T, prob=c(pe,1-pe))
	       paisag[,,t][paisag[,,(t-1)]==0]<-sample(c(0,1),c*l-sum(paisag[,,(t-1)]), replace=T, prob=c(1-pi,pi))
	       resultado[t-1]=sum(paisag[,,t])/(c*l)
	      }

	F=pi/(pi+pe)

	plot(1:tf,c(fi,resultado),type="l",xlab="Tempo",ylab="Fração de manchas ocupadas",
	ylim=c(0,1),main=paste("Chuva de Propágulos","\n c=",c," l=",l," fi=",fi," pi=",pi," pe=",pe),font.lab=2,lwd=2)
	abline(h=F,col=2,lwd=2,lty=2)
	
      return(paisag)
	}
</code>

Maravilha! A função acima tem 6 argumentos (tf, c, l, fi, pi e pe), mas para que a função funcione é preciso atribuir um valor para cada argumento, como no exemplo abaixo:
<code>
metapop(tf=100,c=10,l=10,fi=0.5,pi=0.3,pe=0.15)
</code>
<box 50%|Parâmetros da função metapop>
  * tf - número total de passos no tempo
  * c - número de colunas da matriz //paisag//
  * l - número de linhas da matriz //paisag//
  * fi - fração de manchas ocupadas inicial
  * pi - probabilidade de colonização
  * pe - probabilidade de extinção
</box>

<box 80% red |Exercício>
Agora é com vocês! Tentem responder às seguintes perguntas:
  * O que acontece você rodar mais de uma vez a função com os mesmos parâmetros ? Por quê?((Você deve entender isto para responder várias das questões seguintes.))
  * O que acontece com a variação na fração de manchas ocupadas quando:
      *  o número de manchas é muito grande? E quando é muito pequeno?
      * e se as probabilidades de extinção e colonização são muito altas?
      * e se forem muito baixas?
  * E quando a fração inicial (fi) é muito diferente da fração em equilíbrio (F)?
  * Quais condições levam à extinção da população na paisagem neste modelo?
    
</box> 
== DICA ==

Para responder às perguntas acima você precisa comparar resultados de diferentes simulações.  Abrir espaço para mais de um gráfico na janela gráfica do R vai ajudar bastante. Para ter 4 gráficos digite na linha de comando do R:

<code>
par(mfrow=c(2,2))
</code>
   
Os dois números indicam o número de linhas e colunas desejadas, na janela gráfica. No caso, pedimos duas linhas e duas colunas, ou seja, espaço para quatro gráficos. Para voltar a um gráfico apenas na janela digite:

<code>
par(mfrow=c(1,1))
</code>
 


Para vocês que chegaram vivos até aqui, uma recompensa. Rode os comandos abaixo e descubra:
<code>

anima=function(dados){
        x11()
	for(i in 1:dim(dados)[3])
        {
	image(dados[,,i], main=("Ocupação de manchas"),sub=paste("simulação no.= ", i), col=c("white","red"), bty="n",xaxt='n',yaxt='n')
	grid(dim(dados)[1],dim(dados)[2])
	Sys.sleep(.2)
	}
}


simula1<- metapop(20,10,10,1.0, 0.4,0.2)
anima(simula1)


</code>

Para fazer rodar a animação você precisa apenas salvar o resultado da função //metapop//  e chamar esse objeto na função anima. Para modificar parmâmetros do modelo é só rodar outra simulação.
<code>
simula2=metapop(tf=25,c=100,l=100,fi=.01,pi=0.2,pe=0.5)
anima(simula2)
</code>