Aedes Aegypti ABM

clima-1516.csv

dengue-ABM.R

rm()
set.seed(1000)
setwd("~/Analysis/EmAnalise/Mauricio/Aedes")
# setwd("~/Documents/Aedes")
#dadosClima <- read.table("dados_clima_v2.txt", sep=",", head = T )
#dadosClima <- read.table("dados_clima_v3.txt",  head = T )
dadosClima <- read.table("dados_climaPF2015_2016.csv",  head = T, sep = "," )
dadosClima$Temp.Comp.Media <- (dadosClima$TMAX + dadosClima$TMIN) / 2

colnames(dadosClima) <- c('DATE','SRAD','TMAX','TMIN','Precipitacao', 'Temp.Comp.Media')
#===============================================
# (1) parametrosGlobais
#===============================================
# numero máximo de registros de clima do arquivo de dados
#dataFim <- 629 VEJAM MAIS ABAIXO
dataFim <- nrow(dadosClima)

timeStep <- 1 # d
steps <- 60

taxaMortalidadeOvo <- 0.15
taxaMortalidadeLarva <- 0.05
taxaMortalidadePupa <- 0.05
taxaMortalidadeAdulto <- 0.05

taxaCrescimentoOvo <- 1
taxaCrescimentoLarva <- 0.99
taxaCrescimentoPupa <- 0.6
taxaCrescimentoAdulto <- 0

#temperatura <- 27

#morte <- 0.02

#limite de idade para mudança de fase
idadeLimiteOvo <- 2
idadeLimiteLarva <- 4
idadeLimitePupa <- 6

idadeAdulto <- 8
idadeMorte <- 38

# quantidade de ovos por dia (newaedes)
ovosPorDia <- 10

# quantidade de ovos por dia (newaedes)
# quantidade de ovos para temperatura < 18 graus e > 34 (temperaturas extremas)
ovos_tempExt <- 10
# quantidade de ovos para temperatura > 18 graus e <= 21
ovos_temp18g <- 20
# quantidade de ovos para temperatura >= 22 graus e <= 25
ovos_temp22g <- 30
# quantidade de ovos para temperatura == 26
ovos_temp26g <- 40
# quantidade de ovos para temperatura > 26 e < 30
ovos_temp30g <- 50
# quantidade de ovos para temperatura > 30 graus e <= 34
ovos_temp34g <- 30


#
# template para cada individuo
## fase 1 = ovo 2 = larva 3 = pupa 4 = adulto
newaedes <- data.frame(idade = 0,
                         fase = 1,
                       procuraHumano = 0)

#===============================================
# (2) Trabalhadinha nos Dados do arquivo
#===============================================
library(plyr)
novosDadosClima <- dadosClima
#novosDadosClima <- ddply(dadosClima, .(Data), summarise, Precipitacao = sum(Precipitacao), Temp.Comp.Media = sum(Temp.Comp.Media))
#novosDadosClima$Data <- as.Date(novosDadosClima$Data, "%d/%m/%Y")
#novosDadosClima <- novosDadosClima[order(as.Date(novosDadosClima$Data, format="%d/%m/%Y")),]

#===============================================
# (3) Nova Soma Chuva
#===============================================


#retira 10% de agua do primeiro dia que passou, 20 do segundo, 30 do terceiro e assim por diante.. (10 dias)
somaChuva <- function(timestp){
  
  fatorEvap <- seq(0,1,0.1)
  i <- dataFim-steps
  
  sumTotal <- sum( novosDadosClima$Precipitacao[((i+timestp)-10):((i+timestp))] * fatorEvap )
  #sumTotal <- sum( dadosClima$Precipitacao[i:(i+(timestp))])
  sumTotal
}

calcTemp <- function(timestp){
  # considerando o num de steps estipulados ? verificado
  # em que passo do arquivo de dados climáticos deve começar os cálculos
  # obs: é realizado os steps x 2, pois no arquivo há 2 registros de dados por dia
  dt_atual <- dataFim - steps + timestp
  # a partir da data inicial de verificação é somado quantos dias já se passou 
  # nos timestemp
  somaTM <- sum(novosDadosClima$Temp.Comp.Media[(dt_atual - (3)):dt_atual])
  TMdia <- somaTM / 3
  TMdia
}
#===============================================
# (3) m?todos do ciclo de vida dos individuos
#===============================================
#Soma Chuva
#somaChuva <- function(timestp){
  # considerando o num de steps estipulados ? verificado
  # em que passo do arquivo de dados climáticos deve começar os cálculos
  # obs: ? realizado os steps x 2, pois no arquivo h? 2 registros de dados por dia
  #i <- dataFim-(steps*2) 
  #dataFim <- dataFim+1
  
  # a partir da data inicial de verifica??o ? somado quantos dias j? se passou 
  # nos timestemp
  
  # A soma dos dias anteriores - 10%;
  
  # se i esta no dia atual 549 entao soma a 10% da precipitacao do dia
  #incrementa a porcentagem
  # quando chegar no decimo dia a porcentagem volta ser 10%
  #soma da chuva sera sempre incrementada
  #sumTotal <- sumTotal - sumTotal*0.1

  #sumTotal <- sum( dadosClima$Precipitacao[i:(i+(timestp*2))]*0.9 )
  #sumTotal
#}



## funcao crescer envelhece em 1 o individuo e muda a fase conforme as idades limite
crescer <- function(inds, numdias){
  ninds <- length(inds$idade)
  inds$idade <- inds$idade + timeStep
  
  newinds <- NULL
  
  #VERIFICA FASE
  print(paste(" Numero de Inds",ninds))
  for (i in 1: ninds) {
    # Evolu??o do Ovo
    
    if (inds$fase[i] == 1) {
      
      #Eclosao do ovo
      
      if (inds$idade[i] >= idadeLimiteOvo && somaChuva(numdias) >= 30){
        inds$fase[i] = inds$fase[i] + 1;
      }
      
    }
    # Evolucao da Larva
    if (inds$fase[i] == 2){
      if (inds$idade[i] >= idadeLimiteLarva){
        inds$fase[i] = inds$fase[i] + 1;
      }
    }
    
    # Evolução da Pupa
    if (inds$fase[i] == 3){
      if (inds$idade[i] >= idadeLimitePupa){
        inds$fase[i] = inds$fase[i] + 1;
      }
    }
    
    
    # Adulto - Mosquito
    if (inds$fase[i] == 4){
      #ProcuraHumano vai verificar se o mosquito se alimentou, considerando
      # que todos os dias há alimento disponível
      #assim a cada 3 dias ele ir? colocar ovo
      if (inds$procuraHumano[i] < 3){
        inds$procuraHumano[i] = inds$procuraHumano[i] + 1
      }
      if (inds$procuraHumano[i] == 3){
        tempMedia <- calcTemp(numdias)
        #Analise da temperatura dos últimos dias, referente ao tempo da alimentação
        if (tempMedia < 18 || tempMedia > 34){
          for (j in 1:ovos_tempExt) {
            ## a lista de novos individuos recebe um novo aedes
            newinds <- rbind(newinds,
                             newaedes)
          }
        }
        if (tempMedia >= 18 && tempMedia <= 21){
          for (j in 1:ovos_temp18g) {
            ## a lista de novos individuos recebe um novo aedes
            newinds <- rbind(newinds,
                             newaedes)
          }
        }
        if (tempMedia >= 22 && tempMedia <= 25){
          for (j in 1:ovos_temp22g) {
            ## a lista de novos individuos recebe um novo aedes
            newinds <- rbind(newinds,
                             newaedes)
          }
        }
        if (tempMedia == 26){
          for (j in 1:ovos_temp26g) {
            ## a lista de novos individuos recebe um novo aedes
            newinds <- rbind(newinds,
                             newaedes)
          }
        }
        if (tempMedia > 26 && tempMedia < 30){
          for (j in 1:ovos_temp30g) {
            ## a lista de novos individuos recebe um novo aedes
            newinds <- rbind(newinds,
                             newaedes)
          }
        }
        if (tempMedia >= 30 && tempMedia <= 34){
          for (j in 1:ovos_temp34g ) {
            ## a lista de novos individuos recebe um novo aedes
            newinds <- rbind(newinds,
                             newaedes)
          }
        }
        inds$procuraHumano[i] <- 0
      }
    }
    
  }
  rbind(inds, newinds)
  
}



##
## Os métodos de morrer mantem somente os individuos que passam nos teste
## Subset receber dois paramentros, o primeiro é a lista de elementos
## O segundo é o parametro para saber quais elementos serao mantidos
##
## os metodos morrer abaixo usam as taxas de mortalidade para saber quem ? mantido
## fase 1 = ovo 2 = larva 3 = pupa 4 = adulto
## o método abaixo recebe os individuos, a taxa de mortalidade e a fase que será aplicada a mortalidade
morrer <- function(inds,taxa,fase) {
  # captura a taxa de mortalidade assim podemos alterar conforme precipitacao ou temperatura e não perder o valor original
  taxaMortalidade <- taxa
  # separa os individuos que não ser testados para morrer
  indsOutros <- subset(inds, inds$fase != fase)
  # separa os individuos que devem ser testados para morrer
  indsTestados   <- subset(inds, inds$fase == fase)
  # aplica a taxa de mortalidade
  if (fase == 1){
    indsTemp <- subset(indsTestados, runif(indsTestados$fase) > taxaMortalidade) 
  } else {
    indsTemp <- subset(indsTestados, runif(indsTestados$fase) > taxaMortalidade & indsTestados$idade <= idadeMorte)
  }
  # agrupa os que não morreram com os que não foram testados para morrer
  rbind(indsOutros,  indsTemp)
}

#===============================================
# (4) Inicializa 10 individuos 
#===============================================
inds <- data.frame(idade = c(7, 14, 1, 11, 8,
                           27, 2, 20, 7, 20),
                   fase = c(1, 2, 3, 4,
                            3, 1, 2, 4,
                            2, 3),
                   procuraHumano = c(0, 0, 0, 3,
                                     0, 0, 0, 1,
                                     0, 0))
#===============================================
# (5) life loop
#===============================================
total.ovos <- NULL
total.larvas <- NULL
total.pupas <- NULL
total.adultos <- NULL
total.chuva <- NULL
for (k in 1:steps) {
  print(paste("timestep",k))
  # var K ? necessário na função crescer para que a função somaChuva 
  # verifique em que timestemp esta para acumular os dias que passaram
  
  # Quanto de chuva acumulada o dia?
  print(paste(" Chuva",somaChuva(k)))
  
  #tche.. se nao entrar aqui.. acabou-se a vida.
  if(dim(inds)[1] > 0){
    
    inds <- crescer(inds, k)  
    ## morrer ovos
    inds <- morrer(inds,taxaMortalidadeOvo,1)
    ## morrer larva
    inds <- morrer(inds,taxaMortalidadeLarva,2)  
    ## morrer pupa
    inds <- morrer(inds,taxaMortalidadePupa,3)   
    ## morrer adulto
    inds <- morrer(inds,taxaMortalidadeAdulto,4)    
    # guardando os totais da população
    total.chuva <- c(total.chuva,somaChuva(k))
    total.ovos <- c(total.ovos,length(subset(inds,inds$fase == 1)$fase))
    total.larvas <- c(total.larvas,length(subset(inds,inds$fase == 2)$fase))
    total.pupas <- c(total.pupas,length(subset(inds,inds$fase == 3)$fase))
    total.adultos <- c(total.adultos,length(subset(inds,inds$fase == 4)$fase))
  }
}
#===============================================
# (6) results and graphics
#===============================================

par(mfrow=c(1,1))
simulation <- data.frame(total.ovos,total.larvas,total.pupas,total.adultos,total.chuva)

# Trocado os steps pelo numero de rodadas até que morreram todos. Se sobreviverem ate o numero de steps.. BLZ!
#dt <- seq(1,steps,1)

dt <- seq(1, nrow(simulation) ,1)

colnames(simulation) <- c('Ovos','Larvas','Pupas','Adultos','Chuva')
cores <- c('pink','orange','green','red','blue')
matplot(dt, simulation, type = "l", xlab = "Tempo", ylab = "Quantidade",
        main = "População", lwd = 1, lty = 1, bty = "l", col = cores)
legend("topleft", c("Ovos", "Larvas", "Pupas","Adultos","Chuva"), bty="n" , col = cores, fill=cores)