#' @import foreach doRNG
#'
simHM.customEmigrRuleWeight <- function(x, network, sim.number, num.cores, fill.time){
if (fill.time == F){
parallelCustomMigr <- 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
state.var <- x$ssaObjet$state.var
emigrRule <- x$ssaObjet$emigrRule
probWeights <- x$ssaObjet$probWeights
# 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)
# eval emigrants
emigrants[, "emi"] <- floor(apply(emigrants, 1, function(x){
temp <- sim.result[tempo, paste(state.var, x[1], sep ='')]
colnames(temp) <- state.var
return(eval(parse(text=emigrRule), temp))}))
if (sum(emigrants$emi) != 0){
### sampling from nodes ###
# vector with emigrants tagged by its state, e.g., S or I
sampled <- apply(emigrants, 1,
function(x){
num <- sim.result[tempo, paste(state.var, x[1], sep = '')]
colnames(num) <- state.var
probW <- lapply(probWeights, function(x) eval(parse(text=x), num))
sampled <- sample(sample(rep(state.var, num), x[3], replace = F,
prob = rep(probW[state.var], num)))})
if (is.matrix(sampled) == T)
sampled <- as.list(data.frame(sampled, stringsAsFactors = F))
names(sampled) <- emigrants[,1]
# ----------- Randomly distributing individuals -------------
# connected.nodes is a 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)
connected.nodes[, state.var] <- 0
# distribution of individuals
for(donor.node in sample(emigrants[, from])){
first <- 1
last <- 0
d <- which(emigrants[, from] == donor.node)
for(reciever.node in which(connected.nodes[ , from] == donor.node)){
a <- ceiling(connected.nodes[reciever.node, arc] * emigrants$emi[d] / emigrants[d, arc])
last <- last + a
connected.nodes[reciever.node, state.var] <-
apply(as.matrix(state.var), 1,
function(x){
length(which(sampled[as.character(donor.node)][[1]][first:last] == x))})
first <- first + a
}
}
# balancing and # updating state variables
connected.emigrants <- stats::aggregate(connected.nodes[, state.var],
by = list(connected.nodes[, from]), sum)
connected.imigrants <- stats::aggregate(connected.nodes[, state.var],
by = list(connected.nodes[, to]), sum)
ssaObject$x0[as.vector(apply(as.matrix(state.var), 1, function(x)
paste(x, connected.emigrants[,'Group.1'], sep = '')))] <- as.vector(t(apply(connected.emigrants, 1, function(x){
ssaObject$x0[paste(state.var, x['Group.1'], sep = '')] - as.numeric(x[state.var])})))
ssaObject$x0[as.vector(apply(as.matrix(state.var), 1, function(x)
paste(x, connected.imigrants[,'Group.1'], sep = '')))] <- as.vector(t(apply(connected.imigrants, 1, function(x){
ssaObject$x0[paste(state.var, x['Group.1'], sep = '')] + as.numeric(x[state.var])})))
}
# 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 = ssaObject$ssa.method$method,
tau = ssaObject$ssa.method$tau, f = ssaObject$ssa.method$f,
epsilon = ssaObject$ssa.method$epsilon, nc = ssaObject$ssa.method$nc,
hor = ssaObject$ssa.method$hor, dtf = ssaObject$ssa.method$dtf,
nd = ssaObject$ssa.method$nd)$data
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
ssaObject$x0[names(out.sim)] <- out.sim
}
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% (parallelCustomMigr())
parallel::stopCluster(cl)
sim.result <- do.call(rbind.data.frame, sim.result)
return(sim.result)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.