R/sim_SI_wo_demogr_migr.R

In hybridModels: An R Package for the Stochastic Simulation of Disease Spreading in Dynamic Networks

#' @import foreach doRNG
#' 
simHM.siWoDemogrMigr <- function(x, network, sim.number, num.cores){
  
  siWoDemogrMigr <- function(){
    
    # making it readable
    sim.result <- x$results
    ssaObject <- x$ssaObjet
    from <- x$ssaObjet$var.names$from
    arc <- x$ssaObjet$var.names$arc
    to <- x$ssaObjet$var.names$to
    Time <- x$ssaObjet$var.names$Time
    
    # starting
    sim.result$sim <- sims
    
    for(tempo in 1:length(ssaObject[['mov.dates']])){
      
      # Extracting the total number of migrants per node per tempo
      # emigrants
      emigrants <- 
        stats::aggregate(network[which(network[, Time] == ssaObject[['mov.dates']][tempo]),
                                  c(from,arc)][, arc],
                  by = list(network[which(network[, Time] == ssaObject[['mov.dates']][tempo]),
                                            c(from,arc)][, from]),
                  FUN = sum)
      colnames(emigrants) <- c(from, arc)
      
      # imigrants
      imigrants <-
        stats::aggregate(network[which(network[, Time] == ssaObject[['mov.dates']][tempo]),
                                  c(to,arc)][, arc],
                  by = list(network[which(network[, Time] == ssaObject[['mov.dates']][tempo]),
                                            c(to,arc)][, to]),
                  FUN = sum)
      colnames(imigrants) <- c(to, arc)
      
      # Creating a Vector with current infected/Suscetibles/Efetivo/Traded animals
      s.individuals <- as.vector(sim.result[which(sim.result[, Time] == ssaObject[['mov.dates']][tempo]),
                                            ssaObject$IDs$S.ID], mode = 'integer')
      names(s.individuals) <- ssaObject$IDs$nodes.ID
      i.individuals <- as.vector(sim.result[which(sim.result[, Time] == ssaObject[['mov.dates']][tempo]),
                                            ssaObject$IDs$I.ID], mode = 'integer')
      names(i.individuals) <- ssaObject$IDs$nodes.ID
      n.individuals <- as.vector(sim.result[which(sim.result[, Time] == ssaObject[['mov.dates']][tempo]),
                                            ssaObject$IDs$N.ID], mode = 'integer')
      names(n.individuals) <- ssaObject$IDs$nodes.ID
      
      
      # ---------------- Hypergeom Wiki--------------------
      # mean da hypergeom  <- n*(m/N}, N is poputiation size
      # m, number of successes states in the population
      # n, number of draws
      # --------------------------------------------------
      infected.emigrants <- apply(cbind(i.individuals[as.character(emigrants[, from])],
                                        s.individuals[as.character(emigrants[, from])],
                                        emigrants[, arc]),
                                  1, function(x) {stats::rhyper(1,x[[1]],x[[2]],x[[3]])})
      
      susceptible.emigrants <- emigrants[, arc] - infected.emigrants
      
      # -------------        Expensive Option        --------------
      # ------------- Randomly distributing infected --------------
      # connected.nodes is a vector data frame with the connected nodes in the time tempo
      connected.nodes <-
        stats::aggregate(network[which(network[, Time] == ssaObject[['mov.dates']][tempo]),
                                  c(from,arc)][, arc],
                  by = list(network[which(network[, Time] == ssaObject[['mov.dates']][tempo]),
                                            c(from,arc)][, from],
                            network[which(network[, Time] == ssaObject[['mov.dates']][tempo]),
                                            c(to,arc)][, to]),
                  FUN = sum)
      colnames(connected.nodes) <- c(from, to, arc)
      
      # individual.emigrants is a list of traded individuals (susceptible = zero and infected = 1 in a random 
      # order in each node)
      individual.emigrants <- apply(cbind(rbind(susceptible.emigrants, infected.emigrants,
                                                emigrants = emigrants[, arc]), foo = 0), 2,
                                    function(x){sample(c(rep.int(0, x[1]), rep.int(1, x[2])), x[3],
                                                       replace=FALSE)})
      
      connected.nodes <- cbind.data.frame(connected.nodes, i.in = integer(length(connected.nodes[, to])))
      
      # distribution of infected and susceptible
      for(donor.node in emigrants[, from]){
        
        first <- 1
        last <- 0
        
        for(reciever.node in which(connected.nodes[ , from] == donor.node)){
          
          last <- last + connected.nodes[reciever.node, arc]
          
          connected.nodes[reciever.node, 'i.in'] <- sum(individual.emigrants[[as.character(donor.node)]][first:last])
          
          first <- first + connected.nodes[reciever.node, arc]
          
        }
      }
      connected.nodes <- cbind.data.frame(connected.nodes,
                                          s.in = connected.nodes[, arc] - connected.nodes[, 'i.in'])
      # ----------------------------------------------------------
      
      # first balance
      emigrants <- cbind.data.frame(emigrants, i.out = sapply(individual.emigrants,sum)[-length(individual.emigrants)])
      emigrants <- cbind.data.frame(emigrants, s.out = emigrants[, arc] - emigrants$i.out)
      
      imigrants <- cbind.data.frame(imigrants, 
                                    stats::aggregate(connected.nodes[, c('i.in', 's.in')],
                                              by = list(connected.nodes[, to]), sum))
      imigrants <- imigrants[, -which(names(imigrants) == 'Group.1')]
      
      # second balance
      s.individuals[as.character(emigrants[, from])] <- s.individuals[as.character(emigrants[, from])] - emigrants$s.out
      s.individuals[as.character(imigrants[, to])] <- s.individuals[as.character(imigrants[, to])] + imigrants$s.in
      
      i.individuals[as.character(emigrants[, from])] <- i.individuals[as.character(emigrants[, from])] - emigrants$i.out
      i.individuals[as.character(imigrants[, to])] <- i.individuals[as.character(imigrants[, to])] + imigrants$i.in
      
      n.individuals[as.character(emigrants[, from])] <- n.individuals[as.character(emigrants[, from])] - emigrants[, arc]
      n.individuals[as.character(imigrants[, to])] <- n.individuals[as.character(imigrants[, to])] + imigrants[, arc]
      
      if(x$ssaObjet$pop.correc == TRUE){
        
        ## Population correction +nout-nin #####
        correc <- integer(ssaObject$number.nodes)
        names(correc) <- names(n.individuals)
        correc[as.character(imigrants[, to])] <- imigrants[, arc]
        correc[as.character(emigrants[, from])] <- - emigrants[, arc]
        
        phi <- ceiling(i.individuals * (correc / (s.individuals + i.individuals)))
        
        names(i.individuals) <- paste('I', names(i.individuals), sep = '')
        names(n.individuals) <- paste('N', names(n.individuals), sep = '')
        
        # Alocating N and I
        ssaObject$x0 <- i.individuals[names(ssaObject$x0)] + phi[gsub("^.","",names(ssaObject$x0))]
        ssaObject$parms[names(n.individuals)] <- as.vector(sim.result[1, names(n.individuals)], mode = 'numeric')
        
      } else{
        
        # Alocating N and I
        names(i.individuals) <- paste('I', names(i.individuals), sep = '')
        names(n.individuals) <- paste('N', names(n.individuals), sep = '')
        
        ssaObject$x0 <- i.individuals[names(ssaObject$x0)]
        ssaObject$parms[names(n.individuals)] <- n.individuals
        
      }
      
      # checking whether will be another trade
      out.sim <- GillespieSSA::ssa(x0 = ssaObject$x0, a = ssaObject$propFunction, nu = ssaObject$sCMatrix,
                                   parms = ssaObject$parms, tf = ssaObject$time.diff[tempo],
                                   censusInterval = 0, verbose = FALSE, method = "D")$data
      
      #     out.sim <- ssa(x0=ssaObject$x0, a=ssaObject$propFunction, nu=ssaObject$sCMatrix,
      #                    parms=ssaObject$parms, tf = ssaObject$tFinal, method="BTL",f=100000);
      #     
      #     
      #     out.sim <- ssa(x0=ssaObject$x0, a=ssaObject$propFunction, nu=ssaObject$sCMatrix,
      #                    parms=ssaObject$parms, tf = ssaObject$tFinal, method="OTL",
      #                    nc=50,epsilon=0.02,dtf=10,nd=100);
      #     out.sim <- ssa1(x0=ssaObject$x0, a=ssaObject$propFunction, nu=ssaObject$sCMatrix,
      #                     parms=ssaObject$parms, tf = ssaObject$tFinal, method="OTL",censusInterval=Inf,
      #                     nc=30,epsilon=0.005,dtf=10,nd=100);
      
      out.sim <- out.sim[(length(out.sim[,1])-2), -1]
      
      sim.result[which(sim.result[, Time] == ssaObject[['mov.dates']][tempo]) + 1,
                 names(out.sim)] <- out.sim
      
      sim.result[which(sim.result[, Time] == ssaObject[['mov.dates']][tempo]) + 1,
                 paste('S', names(s.individuals), sep = '')] <-
        ssaObject$parms[gsub('I','N',names(out.sim))] - out.sim
      
      sim.result[which(sim.result[, Time] == ssaObject[['mov.dates']][tempo]) + 1,
                 names(n.individuals)] <- ssaObject$parms[names(n.individuals)]
    }
    
    return(sim.result)
  }
  
  if(num.cores == 'max') {
    num.cores <- parallel::detectCores() 
  }
  
  cl <- parallel::makeCluster(num.cores, type = "PSOCK")
  doParallel::registerDoParallel(cl)
  sims <- NULL
  sim.result <- foreach(sims = 1:sim.number, .verbose=FALSE, .inorder=FALSE,
                        .packages = 'GillespieSSA') %dorng% (siWoDemogrMigr())
  
  parallel::stopCluster(cl)
  
  sim.result <- do.call(rbind.data.frame, sim.result)
  
  return(sim.result)
  
}