=====  Dinâmica de Metapopulações =====
==== Modelo Chuva de Propágulos ====


Nesta aula, iremos simular a dinâmica de uma metapopulação onde as taxas de extinção e de colonização são constantes (modelo de chuva de propágulos). Neste caso há sempre uma fonte constante de imigrantes que podem colonizar qualquer mancha vazia.\\

Primeiro vamos conhecer um novo comando do R que cria uma matriz de dados tridimensional. O comando é **array**. Teste a seguinte linha para ver um array de 3 dimensões.\\

  array(NA,dim=c(4,4,3))
  
Veja que o R criou 3 matrizes de 4 linhas por 4 colunas. Essa matriz assemelha-se a um bolo com 3 camadas, porém quando o objeto é chamado o R mostra as camadas separadamente. No exercício abaixo a primeira dessas 3 matrizes representará o tempo t = 1, e cada célula da matriz 
representará uma mancha no espaço, que poderá ou não estar ocupada.\\

  * Crie uma função chamada **metapop** com os seguintes argumentos: 
      * nl: número de linhas;
      * nc: número de colunas;
      * tmax: número máximo de vezes que você deseja simular a dinâmica de metapopulações;
      * psi0: proporção de manchas ocupadas no início da simulação;
      * cln: probabilidade de colonização de manchas vazias;
      * ext: probabilidade de extinção local em manchas ocupadas
  * O número de linhas nl e o número de colunas nc determinarão quantas manchas potencialmente ocupáveis haverá na metapopulação (nc*nl = número de manchas). Use os seguintes comandos para começar a sua função:

<code>
  ## crie esses objetos para ir testando as linhas da sua função enquanto cria
tmax = 20; nc = 10; nl = 10; psi0 = 0.5; cln = 0.3; ext = 0.15  
  ## Para criar a função.
metapop<-function(nl,nc,tmax,psi0,cln,ext){   
  ## Isso vai criar tmax matrizes com nl linhas e nc colunas.
resulta<-array(0,dim=c(nl,nc,tmax))  
  ## como vamos usar muito o nl*nc (número de manchas), vamos criá-lo

</code>

Agora, para cada mancha você deverá jogar uma moeda para decidir se a mancha está ou não ocupada no tempo t =0. A probabilidade de estar ocupada deve ser **psi0** (a probabilidade de ocorrência inicial). Para isso use a função ''sample()''.\\ Usando ''sample(c(0,1),1)'' obtemos uma amostra de um valor que pode ser 0 ou 1, com igual probabilidade. Mas nós precisamos amostrar **nl*nc** valores de 0 ou 1, pois temos esse número de manchas.\\

Então vamos usar:
  sample(c(0,1),nl*nc,replace=T )

onde ''replace=T'' significa amostragem com reposição.\\

Porém, ainda falta informarmos a probabilidade de ocorrência. Se amostrarmos com igual probabilidade de obter 0 (mancha desocupada) ou 1 (mancha ocupada) a probabilidade de ocorrência será sempre 0.5, mas não é isso que pretendemos. Queremos ter a opção de alterar a probabilidade de ocorrência. Para isso usamos o argumento  "prob=c(1-psi0,psi0)" que especifica a probabilidade de obter os valores indicados dentro da função ''sample''. Ou seja, ''1-psi0'' será a probabilidade de um 0 ser amostrado e ''psi0'' será a probabilidade de um 1 ser amostrado.\\ 

Então vamos simular a ocupação no tempo t=0 de acordo com psi0:\\

<code>
resulta[ , ,1]<-sample(c(0,1),nl*nc,replace=T,prob=c(1-psi0,psi0))
</code>

A probabilidade de ocorrência **psi0** irá variar nos próximos passos, então crie um objeto chamado de resultados "res.psi" e salve  o **psi0** na posição 1 (**psi0** é o valor inicial de fs). No Gotelli **f** é a fração de manchas ocupadas (número de manchas ocupadas / número de manchas).\\

Use um ''for()'' para simular o que acontece com a metapopulação em cada tempo subseqüente, até que **tmax** seja atingido. Repare que já preencheu o array de resultado na primeira posição relativa ao tempo inicial \\

Agora precisamos criar uma condição para saber se uma mancha que estava ocupada continuará ocupada e se uma mancha que estava desocupada continuará desocupada, ou se será ocupada. 

Veja o que apareçe no objeto resulta[,,1]: apenas valores 0 e 1, significando ausência ou ocupação de mancha no tempo t=1. No tempo t+1 ou no caso t=2 as possibilidades são:
  * Manchas desocupadas (valores 0) podem ser colonizados e passar a 1, ou  continuar 0.
  * Manchas ocupadas (valores 1)podem ter extinção e passsar a 0, ou continuar 1
Para o primeiro caso a probabilidade de 0 passar a 1 é  //cln// e a probabilidade de 1 passar a 0 é //ext//.
Uma forma de fazer isso é:
<code>
  num1=sum(resulta[,,t-1])
  num0=(nl*nc)-num1
  resulta[,,2][resulta[,,1]==1]<-sample(c(0,1),num1,replace=T,prob=c(ext,1-ext))
  ...
</code>
Esse comando pode ser lido da seguinte forma: nas posições do resulta t2 onde resulta t1 são iguais a 1, coloque 0 ou 1 com probabilidade ext e 1-ext respectivamente.

Faça o mesmo para  as posições em t1 iguais a 0, pensando na colonização.
Pronto, agora e só mudar o valores de tempo pelo seu indexador para que em cada ciclo um novo tempo seja calculado. 
Para finalizar o ciclo só falta calcular o res.psi de cada tempo.
  res.psi[t]=sum(resulta[,,t])/(nl*nc)
  

 

<box 100% |Dica avançada>
Sugerimos que use a solução anterior por ser mais intuitiva. Caso queira se aprofundar nos comandos do R, veja a solução abaixo para o mesmo problema:  substituir os dados de ocupação de manchas de um tempo para outro, utilizando as probabilidades de extinção e colonização:

É possível fazer o mesmo com uma única linha de comando utilizando a função 
''ifelse()''. Ela funciona assim: 
 
  ifelse(condição, "faça isso", "se não faça aquilo")


Então tente algo assim:\\

   resulta[,,t]<-ifelse(Passo 1, Passo 2, Passo3) 

Veja abaixo as instruções para os Passos 1, 2 e 3.''\\


**Instruções para os passos 1, 2 e 3 acima.**\\

Na função ''ifelse'' acima você deverá, no Passo 1, criar uma condição para saber se uma determinada mancha estava ou não ocupada no tempo t-1.\\ 
Caso a mancha estivesse desocupada você deverá, no Passo 2, dizer o que acontecerá com ela no tempo t. Para isso gere um número aleatório de uma distribuição uniforme e veja se esse valor é menor que a taxa de colonização, se o valor for menor a mancha deverá ser ocupada, se for maior ela permanece desocupada. Vamos ver como fazer isso:\\
	
**Atenção**, as duas linhas abaixo são para você testar no R e entender os comandos, estas linhas não devem estar no seu script final.\\

  runif(1)<cln  
Veja que o resultado será TRUE ou FALSE, mas nós queremos que seja 0 ou 1''\\

Precisamos transformar o TRUE e o FALSE em valores numéricos, usando a função ''as.numeric()''\\

  as.numeric(runif(1)<cln) 	
veja que o resultado agora será 0 ou 1.''\\

Agora você já sabe como fazer o Passo 2, mas lembre que no ''runif(1'') você não deverá ter apenas 1 valor, e sim nc*nl valores, ou seja, um valor pra cada mancha.\\

Passo 3: Aprenda com o Passo 2.  :-)\\

Fim das instruções para os passos 1, 2 e 3.\\
</box>

===== Representação gráfica =====

Agora faça um gráfico, **dentro da sua função**, que mostre a trajetória da proporção de manchas ocupadas ao longo do tempo. Acrescente uma reta mostrando o valor da fração de manchas ocupadas no equilíbrio **f.equilibrio**.\\
Equação 4.4 do Gotelli: f=Pi/(Pi+Pe) 

Primeiro calcule o f.equilibrio

  f.equilibrio<- ??

<code>
plot(1:tmax),fs ,type="l",xlab="tempo",ylab="fs",ylim=c(0,1),main="Metapopulação ")
abline(h=f.equilibrio, col=2)
</code>
===== Animação de um conjunto de manchas no tempo =====

Agora vamos criar uma animação para olhar nossos resultados ao longo do tempo. Este mostra uma animação da dinâmica da metapopulação ao longo do tempo.\\
Primeiro salve o resultado de uma simulação da função metapop e crie um gráfico do tempo 1.
 <code>
   image(dados[,,i], main=("Ocupação de manchas"),col=c("white","red"))
   grid(nl,nc)
 </code>  

Agora é só criar o ciclo com //for// para que os gráficos de cada tempo sejam sobrepostos no mesmo dispositivo gráfico.
Para isso vamos criar uma função //anima.meta//:
  - O único argumento é o dado que será igual ao resultado da função metapop.
  - crie um ciclo de 1 até tmax que é igual ao tamanho da dimensão três do seus dados
  - faça com que em cada ciclo seja feito um gráfico usando //image// dos dados do tempo t
  - crie um lapso de tempo para que os gráficos possam ser visualizados antes de trocar, colocando o comando abaixo:
  Sys.sleep(.5)

  * pronto!

 
Use sua função e altere os parâmetros iniciais para ver se a metapopulação persiste, dependendo das condições iniciais.  Se tiver tempo, procure alterar a sua função para implementar os outros modelos descritos na aula téorica.


<box 90% blue| Uma animação mais informativa>
Uma forma de representar a dinâmica espacial da metapopulação é incluir no mapa não só os espaços ocupados e vagos mas entre esses aqueles que foram recentemente colonizados ou extintos. Ou seja, teremos quatro tipos de informações em cada momento: ocupado antigo, ocupado recente, vago, extinção recente.
<code>

meta.anima2=function(dados)
{
nsim=dim(dados)[3]
image(dados[,,1], main=("Ocupacao de manchas"),col=c("white","green"),sub="tempo=1")
	for(i in 2:nsim)
	{
	conta12=dados[,,(i-1)]+ (2*dados[,,i])
	image(conta12, main=("Ocupacao de manchas"),col=c("white", "red","green",
"darkgreen"),sub=paste("tempo= ",i))
	Sys.sleep(.5)
	}
}


</code>

</box> 


--------------------

==== Modelo Colonização Interna ====

Podemos eliminar do modelo anterior o pressuposto de uma chuva de propágulos constante e fazer com que a colonização seja uma função do número de lugares ocupados. Em uma formulação simples desse modelo, a fonte de propágulos é unicamente interna (sistema fechado) e a probabilidade de colonização varia de forma linear à proporção de lugares ocupados. 
Nossa função começa então diferente com relação aos argumentos necessários:
  * sai o argumento //cln//;
  * entra o argumento //ef//: quanto a proporção de áreas ocupadas afeta a colonização.  

<code>
meta.inter<-function(nl,nc,tmax,psi0,ext, ef)
{   
...
</code>

  *  calcule a probabilidade de colonização no inicio de cada ciclo:
<code>
...
  for(t in 2: tmax)
  {
   cln=ef*res.psi[t-1]
   ...
   ...
  } 
</code>
 
  * modifique o valor da proporção de áreas ocupadas no equilíbrio para:
<code>
...
f.equilibrio=1-(ext/ef)
...   
</code>

  * com essas mudanças a novo modelo deve funcionar!