diff --git a/DESCRIPTION b/DESCRIPTION
index 0ceb01a04f7f4c48222424c5129da7ef757bff36..65158a23e8947530467b1ab926ac0a637acee4b3 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -12,29 +12,30 @@ LazyData: true
Encoding: UTF-8
Depends: R (>= 4.0.0)
Imports:
- SimInf,
- data.table,
- magrittr,
- chron,
- terra,
- tidyterra,
- ggplot2,
- raster,
- sp,
- sf,
- geodata,
- zoo,
- stringdist,
- stringr,
- leaflet,
- RColorBrewer,
- imputeTS,
+ SimInf (>= 9.5.0),
+ data.table (>= 1.14.8),
+ magrittr (>= 2.0.3),
+ chron (>= 2.3-60),
+ terra (>= 1.7-18),
+ tidyterra (>= 0.4.0),
+ ggplot2 (>= 3.4.1),
+ raster (>= 3.6-20),
+ sp (>= 1.6-0),
+ sf (>= 1.0-12),
+ geodata (>= 0.5-3),
+ zoo (>= 1.8-11),
+ stringdist (>= 0.9.10),
+ stringr (>= 1.5.0),
+ leaflet (>= 2.1.2),
+ RColorBrewer (>= 1.1-3),
+ imputeTS (>= 3.3),
parallel,
- doParallel,
- foreach
+ doParallel (>= 1.0.17),
+ foreach (>= 1.5.2),
+ pbapply (>= 1.7-0)
Suggests:
- knitr,
- rmarkdown
+ knitr (>= 1.42),
+ rmarkdown (>= 2.20)
VignetteBuilder:
knitr
RoxygenNote: 7.2.3
diff --git a/NAMESPACE b/NAMESPACE
index 4bf2512649a66452a3b203aeb8f7be576729a696..857f0728af6b1a9db1c33541c2adaeef85b4a924 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,5 +7,15 @@ export(cleanmeteo)
export(create_ldata)
export(diapause)
export(iniState)
+export(runArboRisk)
import(data.table)
+importFrom(SimInf,mparse)
+importFrom(SimInf,run)
+importFrom(SimInf,trajectory)
importFrom(imputeTS,na_ma)
+importFrom(parallel,clusterEvalQ)
+importFrom(parallel,clusterExport)
+importFrom(parallel,detectCores)
+importFrom(parallel,makeCluster)
+importFrom(parallel,stopCluster)
+importFrom(pbapply,pbapply)
diff --git a/R/build_pts_fun.R b/R/build_pts_fun.R
index 333f804113244d0fcd3e2f4b275aa6f70dbf3bfe..14aa2157893f4d9e47a723f44ecd6f62ba2d20a0 100644
--- a/R/build_pts_fun.R
+++ b/R/build_pts_fun.R
@@ -7,7 +7,7 @@
#' @param u0 initial population stage for compartment driven by stochastic processes
#' @param v0 continuous variables. Contains post time step state of vector population as well as daily meteorological data.
#' @param gdata list of global data
-#' @param diapause_interv interval of days defining favourable period for mosquitoes
+#' @param diapause_interv interval of days defining favorable period for mosquitoes
#'
#' @return String containing C code used as pts_fun by SimInf package
#'
@@ -16,7 +16,7 @@
#'
#' @export
-build_pts_fun <- function(u0, v0, gdata, diapause_interv){
+build_pts_fun <- function(u0, v0, gdata, ldata, diapause_interv){
diapause_string <- ""
@@ -31,12 +31,32 @@ build_pts_fun <- function(u0, v0, gdata, diapause_interv){
diapause_string <- paste0(diapause_string, paste("(int)t >=",diapause_interv[i*2+1]," && (int)t < ",diapause_interv[i*2+2]))
}
+
+ # # infection of mosquitoes
+ # # weighted infected pop
+ # wIh <- c()
+ # for(patch in seq(0, nrow(u0)-1)){
+ # wIh %<>% append(.,
+ # paste0(c("ldata[",1 + round(ldata[1,1]) * 3 + patch, # proba that someone from patch i is in patch j
+ # "]*u_0[",2 + patch * ncol(u0),"]/(u_0[", # number of infectious ind. in patch i
+ # 0 + patch * ncol(u0),"] + u_0[", # number of susceptible ind. in patch i
+ # 1 + patch * ncol(u0),"] + u_0[", # number of exposed ind. in patch i
+ # 2 + patch * ncol(u0),"] + u_0[", # number of infectious ind. in patch i
+ # 3 + patch * ncol(u0),"])"), collapse = "")) # number of recovered ind. in patch i
+ # }
+ # wIh %<>% paste(., collapse= "+")
+
###########################################################
### Deterministic transitions for mosquitoes life cycle ###
###########################################################
options(scipen=999)
- pts_fun <- paste('
+ pts_fun <- paste0('
+
+const int nCpmt_u = (int)',ncol(u0), 'L;
+const int nCpmt_v = (int)',nrow(v0), 'L;
+const int * u_0 = &u[-node * nCpmt_u];
+const double * v_0 = &v[-node * nCpmt_v];
// update diapause
v_new[0] = 0;
@@ -50,6 +70,13 @@ build_pts_fun <- function(u0, v0, gdata, diapause_interv){
v_new[i] = ldata[1 + (int)ldata[0] * (i-1) + (int)t];
}
+ double temperature;
+ long kP;
+ long kL;
+ temperature = v_new[', which(rownames(v0) == "temperature") - 1,'];
+ kP = v_new[', which(rownames(v0) == "kP") - 1,'];
+ kL = v_new[', which(rownames(v0) == "kL") - 1,'];
+
// Development and mortality functions
// EGGS
@@ -59,75 +86,75 @@ build_pts_fun <- function(u0, v0, gdata, diapause_interv){
v_new[', which(rownames(v0) == "prevEggs") - 1,'] = u[', which(names(u0) == "Neggs") - 1,']/1000;
// eggs layed by susceptible mosquitoes
int neweggS;
- neweggS =', gdata["gammaAo"], ' * (v[10] * ', gdata["beta1"], ' + v[13] * ', gdata["beta2"], ');
+ neweggS =', gdata["gammaAo"], ' * (v[', which(rownames(v0) == "A1om") - 1,'] * ', gdata["beta1"], ' + v[', which(rownames(v0) == "A2om") - 1,'] * ', gdata["beta2"], ');
// eggs development into larvae and mortality
double fdevE;
- if(v_new[1] > ', gdata["TE"], '){
- fdevE = v_new[0] * ((v_new[1] - ', gdata["TE"], ') / ', gdata["TDDE"], ');
+ if(temperature > ', gdata["TE"], '){
+ fdevE = v_new[', which(rownames(v0) == "z") - 1,'] * ((temperature - ', gdata["TE"], ') / ', gdata["TDDE"], ');
}else{
fdevE = 0;
}
long varEm;
if((', gdata["muE"],' + fdevE) > 1){
- varEm = neweggI + neweggS - v[4];
+ varEm = neweggI + neweggS - v[', which(rownames(v0) == "Em") - 1,'];
} else{
- varEm = neweggI + neweggS - v[4] * (', gdata["muE"],' + fdevE);
+ varEm = neweggI + neweggS - v[', which(rownames(v0) == "Em") - 1,'] * (', gdata["muE"],' + fdevE);
}
// LARVAE
// larvae development into pupae
double fdevL;
- if((', gdata["q1L"], ' * v_new[1] * v_new[1] + ', gdata["q2L"], ' * v_new[1] + ', gdata["q3L"], ') < 0){
+ if((', gdata["q1L"], ' * temperature * temperature + ', gdata["q2L"], ' * temperature + ', gdata["q3L"], ') < 0){
fdevL = 0;
} else {
- fdevL = (', gdata["q1L"], ' * v_new[1] * v_new[1] + ', gdata["q2L"], ' * v_new[1] + ', gdata["q3L"], ');
+ fdevL = (', gdata["q1L"], ' * temperature * temperature + ', gdata["q2L"], ' * temperature + ', gdata["q3L"], ');
}
// larvae mortality
double fmL;
- fmL = (0.0007 * exp((v_new[1]-10)*0.1838) + 0.02) * (1 + v[5]/v_new[2]);
+ fmL = (0.0007 * exp((temperature-10)*0.1838) + 0.02) * (1 + v[', which(rownames(v0) == "Lm") - 1,']/kL);
long varLm;
if((fdevL + fmL) > 1){
- varLm = v[4] * fdevE - v[5] ;
+ varLm = v[', which(rownames(v0) == "Em") - 1,'] * fdevE - v[', which(rownames(v0) == "Lm") - 1,'] ;
} else {
- varLm = v[4] * fdevE - v[5] * (fdevL + fmL) ;
+ varLm = v[', which(rownames(v0) == "Em") - 1,'] * fdevE - v[', which(rownames(v0) == "Lm") - 1,'] * (fdevL + fmL) ;
}
// PUPAE
// pupae development into emerging adult
double fdevP;
- fdevP = (', gdata["q1P"], ' * v_new[1] * v_new[1] + ', gdata["q2P"], ' * v_new[1] + ', gdata["q3P"], ');
+ fdevP = (', gdata["q1P"], ' * temperature * temperature + ', gdata["q2P"], ' * temperature + ', gdata["q3P"], ');
if(fdevP < 0)
fdevP = 0;
// pupae mortality
double fmP;
- fmP = (0.0003 * exp((v_new[1]-10)*0.2228) + 0.02) * (1 + v[6]/v_new[3]);
+ fmP = (0.0003 * exp((temperature-10)*0.2228) + 0.02) * (1 + v[', which(rownames(v0) == "Pm") - 1,']/kP);
long varPm;
if((fdevP + fmP) > 1){
- varPm = v[5] * fdevL - v[6] ;
+ varPm = v[', which(rownames(v0) == "Lm") - 1,'] * fdevL - v[', which(rownames(v0) == "Pm") - 1,'] ;
} else {
- varPm = v[5] * fdevL - v[6] * (fdevP + fmP) ;
+ varPm = v[', which(rownames(v0) == "Lm") - 1,'] * fdevL - v[', which(rownames(v0) == "Pm") - 1,'] * (fdevP + fmP) ;
}
// ADULTS
// aldult mortality
double fmA;
- fmA = 0.0003 * exp((v_new[1] - 10)*0.1745) + 0.025;
+ fmA = 0.0003 * exp((temperature - 10)*0.1745) + 0.025;
if(fmA < ', gdata["muA"], ')
fmA = ', gdata["muA"], ';
// emerging aldult mortality
double fmAEm;
- fmAEm = exp(-', gdata["muEM"], ' * (1+v[6]/v_new[3]));
+ fmAEm = exp(-', gdata["muEM"], ' * (1+v[', which(rownames(v0) == "Pm") - 1,']/kP));
// aldult mortality in research activity
double fmAr;
@@ -135,96 +162,166 @@ build_pts_fun <- function(u0, v0, gdata, diapause_interv){
// aldult gorging activity
double fAg;
- if(v_new[1] > ', gdata["TAG"], '){
- fAg = (v_new[1] - ', gdata["TAG"], ')/', gdata["TDDAG"], ';
+ if(temperature > ', gdata["TAG"], '){
+ fAg = (temperature - ', gdata["TAG"], ')/', gdata["TDDAG"], ';
} else {
fAg = 0;
}
double varAemm;
- varAemm = v[6] * fdevP * ', gdata["sigma"],' * fmAEm - v[7] * (fmA + ', gdata["gammaAem"], ') ;
+ varAemm = v[', which(rownames(v0) == "Pm") - 1,'] * fdevP * ', gdata["sigma"],' * fmAEm - v[', which(rownames(v0) == "Aemm") - 1,'] * (fmA + ', gdata["gammaAem"], ') ;
double varA1hm;
- varA1hm = v[7] * ', gdata["gammaAem"], ' - v[8] * (fmAr + ', gdata["gammaAh"], ') ;
+ varA1hm = v[', which(rownames(v0) == "Aemm") - 1,'] * ', gdata["gammaAem"], ' - v[', which(rownames(v0) == "A1hm") - 1,'] * (fmAr + ', gdata["gammaAh"], ') ;
double varA1gm;
- varA1gm = v[8] * ', gdata["gammaAh"], ' - v[9] * (fmA + fAg) ;
+ varA1gm = v[', which(rownames(v0) == "A1hm") - 1,'] * ', gdata["gammaAh"], ' - v[', which(rownames(v0) == "A1gm") - 1,'] * (fmA + fAg) ;
double varA1om;
- varA1om = v[9] * fAg - v[10] * (fmAr + ', gdata["gammaAo"], ') ;
+ varA1om = v[', which(rownames(v0) == "A1gm") - 1,'] * fAg - v[', which(rownames(v0) == "A1om") - 1,'] * (fmAr + ', gdata["gammaAo"], ') ;
double varA2hm;
- varA2hm = (v[10] + v[13]) * ', gdata["gammaAo"], ' - v[11] * (fmAr + ', gdata["gammaAh"], ') ;
+ varA2hm = (v[', which(rownames(v0) == "A1om") - 1,'] + v[', which(rownames(v0) == "A2om") - 1,']) * ', gdata["gammaAo"], ' - v[', which(rownames(v0) == "A2hm") - 1,'] * (fmAr + ', gdata["gammaAh"], ') ;
double varA2gm;
- varA2gm = v[11] * ', gdata["gammaAh"], ' - v[12] * (fmA + fAg) ;
+ varA2gm = v[',which(rownames(v0) == "A2hm") - 1,'] * ', gdata["gammaAh"], ' - v[', which(rownames(v0) == "A2gm") - 1,'] * (fmA + fAg) ;
double varA2om;
- varA2om = v[12] * fAg - v[13] * (fmAr + ', gdata["gammaAo"], ');
+ varA2om = v[', which(rownames(v0) == "A2gm") - 1,'] * fAg - v[', which(rownames(v0) == "A2om") - 1,'] * (fmAr + ', gdata["gammaAo"], ');
// update mosquito population vector
// eggs
- v_new[4] = v[4] + varEm ;
+ v_new[', which(rownames(v0) == "Em") - 1,'] = v[', which(rownames(v0) == "Em") - 1,'] + varEm ;
// larvae
- v_new[5] = v[5] + varLm ;
+ v_new[', which(rownames(v0) == "Lm") - 1,'] = v[', which(rownames(v0) == "Lm") - 1,'] + varLm ;
// pupae
- v_new[6] = v[6] + varPm ;
+ v_new[', which(rownames(v0) == "Pm") - 1,'] = v[', which(rownames(v0) == "Pm") - 1,'] + varPm ;
// emerging adult
- v_new[7] = v[7] + varAemm ;
+ v_new[', which(rownames(v0) == "Aemm") - 1,'] = v[', which(rownames(v0) == "Aemm") - 1,'] + varAemm ;
// nulliparous host-seeking adult
- v_new[8] = v[8] + varA1hm ;
+ v_new[', which(rownames(v0) == "A1hm") - 1,'] = v[', which(rownames(v0) == "A1hm") - 1,'] + varA1hm ;
// nulliparous gorged adult
- v_new[9] = v[9] + varA1gm ;
+ v_new[', which(rownames(v0) == "A1gm") - 1,'] = v[', which(rownames(v0) == "A1gm") - 1,'] + varA1gm ;
// nulliparous oviposition-site-seeking adult
- v_new[10] = v[10] + varA1om ;
+ v_new[', which(rownames(v0) == "A1om") - 1,'] = v[', which(rownames(v0) == "A1om") - 1,'] + varA1om ;
// parous host-seeking adult
- v_new[11] = v[11] + varA2hm ;
+ v_new[', which(rownames(v0) == "A2hm") - 1,'] = v[', which(rownames(v0) == "A2hm") - 1,'] + varA2hm ;
// parous gorged adult
- v_new[12] = v[12] + varA2gm ;
+ v_new[', which(rownames(v0) == "A2gm") - 1,'] = v[', which(rownames(v0) == "A2gm") - 1,'] + varA2gm ;
// parous oviposition-site-seeking adult
- v_new[13] = v[13] + varA2om;
+ v_new[', which(rownames(v0) == "A2om") - 1,'] = v[', which(rownames(v0) == "A2om") - 1,'] + varA2om;
//Calcul of R0
int atot;
- atot = v[7] + v[8] + v[9] + v[10] + v[11] + v[12] + v[13];
+ atot = u[',which(colnames(u0) == "A1gmI") - 1,
+ '] + u[',which(colnames(u0) == "A1omI") - 1,
+ '] + u[',which(colnames(u0) == "A2hmI") - 1,
+ '] + u[',which(colnames(u0) == "A2gmI") - 1,
+ '] + u[',which(colnames(u0) == "A2omI") - 1,
+ ']+ v[', which(rownames(v0) == "Aemm") - 1,
+ ']+ v[', which(rownames(v0) == "A1hm") - 1,
+ ']+ v[', which(rownames(v0) == "A1gm") - 1,
+ ']+ v[', which(rownames(v0) == "A1om") - 1,
+ ']+ v[', which(rownames(v0) == "A2hm") - 1,
+ ']+ v[', which(rownames(v0) == "A2gm") - 1,
+ ']+ v[', which(rownames(v0) == "A2om") - 1,'];
double taux_survie;
- taux_survie = 1-fmA-((v[11]+v[13])/atot)*fmAr;
+ taux_survie = 1 - ((v[', which(rownames(v0) == "Aemm") - 1,'] + v[', which(rownames(v0) == "A1hm") - 1,'] + v[', which(rownames(v0) == "A1gm") - 1,'] + v[', which(rownames(v0) == "A1om") - 1,'] + v[', which(rownames(v0) == "A2gm") - 1,'])/atot) * fmA - ((v[', which(rownames(v0) == "A2hm") - 1,']+v[', which(rownames(v0) == "A2om") - 1,'])/atot) * fmAr;
double incub_extr;
- incub_extr = 0.11*(v_new[1]*v_new[1])-7.13*v_new[1]+121.17;
+ incub_extr = 0.11*(temperature*temperature)-7.13*temperature+121.17;
double comp_vect;
- if(((-0.0043)*(v_new[1]*v_new[1])+(0.2593*v_new[1])-3.2705) > 0){
- comp_vect = ((-0.0043)*(v_new[1]*v_new[1])+(0.2593*v_new[1])-3.2705);
+ if(((-0.0043)*(temperature*temperature)+(0.2593*temperature)-3.2705) > 0){
+ comp_vect = ((-0.0043)*(temperature*temperature)+(0.2593*temperature)-3.2705);
} else {
comp_vect = 0;
}
- int pop;
- pop = u[0] + u[1] + u[2] + u[3];
+ // Total human population
+ int TP;
+ TP = u[',which(colnames(u0) == "Sh") - 1,'] + u[',which(colnames(u0) == "Eh") - 1,'] + u[',which(colnames(u0) == "Ih") - 1,'] + u[',which(colnames(u0) == "Rh") - 1,'];
+
+ // Infected human population
+ int IP;
+ IP = u[',which(colnames(u0) == "Ih") - 1,'];
double capac_vect;
- if(pop > 0 && atot > 0){
- capac_vect = atot/pop*pow((v[8] + v[11])*', gdata["gammaAh"], '/atot,2)*pow(taux_survie,incub_extr)/-log(taux_survie);
+ if(TP > 0 && atot > 0){
+ capac_vect = atot/TP*pow((v[8] + v[11])*', gdata["gammaAh"], '/atot,2)*pow(taux_survie,incub_extr)/-log(taux_survie);
} else {
capac_vect = 0;
}
int R0;
- if(R0 > pop){
- v_new[15] = pop;
+ if(R0 > TP){
+ v_new[', which(rownames(v0) == "R0") - 1,'] = TP;
} else {
- v_new[15] = (comp_vect*capac_vect)/', gdata["rhoH"],';
+ v_new[', which(rownames(v0) == "R0") - 1,'] = (comp_vect*capac_vect)/', gdata["rhoH"],';
}
+ // proportion of infected mosquitoes weighted by the probability of being in contact (probability to be in the patch)
+ // initialize proportion
+ double wIm;
+ wIm = 0;
+
+ // total number of mosquitoes looking for host in the patch
+ int nMh;
+ // proportion of infected mosquitoes looking for host in the patch
+ double pIMh;
+
+ int j;
+ for(j = 0; j <= ',nrow(u0)-1,'; ++j){
+ nMh = u_0[',which(colnames(u0) == "A2hmI") - 1,' + j * ',ncol(u0),'] + v_0[', which(rownames(v0) == "A1hm") - 1,' + j * ',nrow(v0),'] + v_0[', which(rownames(v0) == "A2hm") - 1,' + j * ',nrow(v0),'];
+ if(nMh > 0){
+ pIMh = (double) u_0[',which(colnames(u0) == "A2hmI") - 1,' + j * ',ncol(u0),'] / nMh;
+ } else {
+ pIMh = (double) 0;
+ }
+
+ pIMh = pIMh * ldata[1 + (int)ldata[0] * 3 + j];
+ wIm = wIm + pIMh;
+ }
+
+ v_new[', which(rownames(v0) == "pIm") - 1,'] = (double) wIm;
+
+
+ // proportion of infected humans weighted by the probability of being in contact (probability to be in the patch)
+ // initialize proportion
+ double wIh;
+ wIh = 0;
+
+
+ // total number of mosquitoes looking for host in the patch
+ int nHtot;
+ // proportion of infected mosquitoes looking for host in the patch
+ double pIH;
+
+ int k;
+ for(k = 0; k <= ',nrow(u0)-1,'; ++k){
+ nHtot = u_0[',which(colnames(u0) == "Sh") - 1,' + k * ',ncol(u0),
+ '] + u_0[',which(colnames(u0) == "Eh") - 1,' + k * ',ncol(u0),
+ '] + u_0[',which(colnames(u0) == "Ih") - 1,' + k * ',ncol(u0),
+ '] + u_0[',which(colnames(u0) == "Rh") - 1,' + k * ',ncol(u0),'];
+ if(nHtot > 0){
+ pIH = (double) u_0[',which(colnames(u0) == "Ih") - 1,' + k * ',ncol(u0),'] / nHtot;
+ } else {
+ pIH = (double) 0;
+ }
+
+ pIH = pIH * ldata[1 + (int)ldata[0] * 3 + k];
+ wIh = wIh + pIH;
+ }
+
+
+ v_new[', which(rownames(v0) == "pIh") - 1,'] = (double) wIh;
+
+
return 1;
')
options(scipen=0)
return(pts_fun)
}
-
-
-
diff --git a/R/build_transitions.R b/R/build_transitions.R
index 3a32baf73b94988b56a04dd1010d52618b63c905..22efbe75bd2a65225ac9047f43b9c7cb1ffe228a 100644
--- a/R/build_transitions.R
+++ b/R/build_transitions.R
@@ -42,8 +42,8 @@ stoch_transitions <- c(
#### MOSQUITOS INFECTION ####
# We do not consider "exposed" stage and temperature-dependent extrinsic incubation period (EIP) -- all females directly become infectious // we do not include recovery
- "@ -> 0.0003 * exp((temperature-10)*0.1745) < 0 ? A1hm * (gammaAh - muA - muR) * Ih/(Sh+Eh+Ih+Rh) * bHM : A1hm * (gammaAh - ((0.0003 * exp((temperature-10) * 0.1745) + muA + muR))) * Ih/(Sh+Eh+Ih+Rh) * bHM -> A1gmI",
- "@ -> 0.0003 * exp((temperature-10)*0.1745) < 0 ? A2hm * (gammaAh - muA - muR) * Ih/(Sh+Eh+Ih+Rh) * bHM : A2hm * (gammaAh - ((0.0003 * exp((temperature-10) * 0.1745) + muA + muR))) * Ih/(Sh+Eh+Ih+Rh) * bHM -> A2gmI",
+ "@ -> 0.0003 * exp((temperature-10)*0.1745) < 0 ? A1hm * (gammaAh - muA - muR) * pIh * bHM : A1hm * (gammaAh - ((0.0003 * exp((temperature-10) * 0.1745) + muA + muR))) * pIh * bHM -> A1gmI",
+ "@ -> 0.0003 * exp((temperature-10)*0.1745) < 0 ? A2hm * (gammaAh - muA - muR) * pIh * bHM : A2hm * (gammaAh - ((0.0003 * exp((temperature-10) * 0.1745) + muA + muR))) * pIh * bHM -> A2gmI",
"A1gmI -> temperature > TAG ? A1gmI * (temperature - TAG) / TDDAG : 0 -> A1omI",
"A1omI -> A1omI * gammaAo -> A2hmI + 95 * Neggs", ## New eggs
diff --git a/R/create_ldata.R b/R/create_ldata.R
index 4c80bb1b7743c677aaa25ba8992cc0e398762887..4b0c56acc8d95c72ea45d832a2912cd5bfa2e722 100644
--- a/R/create_ldata.R
+++ b/R/create_ldata.R
@@ -1,65 +1,76 @@
#' Write local data matrix
#'
-#' @description Internal function to write the pts_fun code driving SimInf package to implement deterministic population dynamics of Aedes Albopictus.
+#' @description Function to define local data
#'
-#' @usage create_ldata(PARCELLE, meteo, gdata, time_max, nYR_init, nodeID)
+#' @usage create_ldata(PARCELLE, METEO, gdata, time_max, nYR_init, nodeID)
#'
#' @param PARCELLE data.frame or data.table
-#' @param meteo data.frame or data.table
+#' @param METEO data.frame or data.table
#' @param gdata list
#' @param time_max numerical value
#' @param nYR_init numerical value
#' @param nodeID string
#'
+#' @importFrom pbapply pbapply
+#' @importFrom parallel makeCluster
+#' @importFrom parallel detectCores
+#' @importFrom parallel clusterExport
+#' @importFrom parallel clusterEvalQ
+#' @importFrom parallel stopCluster
+#'
#' @return matrix
#'
-#' @keywords ldata meteo
+#' @keywords ldata METEO
#'
#' @export
-create_ldata <- function(PARCELLE, meteo, gdata, time_max, nYR_init = 2, nodeID = "ID"){
+create_ldata <- function(PARCELLE, METEO, gdata, mMov, time_max, nYR_init = 1){
## We have seven variables that we want to incorporate in each node.
## Create a matrix for each variable, where each column contains the data for one node.
## Format
+ PARCELLE %<>% copy
PARCELLE %>% setDT
- meteo %>% setDT
+
+ METEO %<>% copy
+ setDT(METEO)
## CHECKS
- try(if(NA %in% match(PARCELLE$STATION, meteo$POSTE))
- stop("Some stations are missing from the meteorological dataset"))
+ try(if(NA %in% match(PARCELLE$STATION, METEO$ID))
+ stop("Some stations are missing from the METEOrological dataset"))
- try(if(meteo[POSTE %in% PARCELLE$STATION, .N, by = POSTE] %>% .[,N] %>% min %>% is_less_than(365))
- stop("At least one year of meteorological data are required for initialization"))
+ try(if(METEO[ID %in% PARCELLE$STATION, .N, by = ID] %>% .[,N] %>% min %>% is_less_than(365))
+ stop("At least one year of METEOrological data are required for initialization"))
- try(if(meteo[POSTE %in% PARCELLE$STATION, .N, by = POSTE] %>% .[,N] %>% min %>% is_less_than(time_max))
- stop("All stations must be associated to a number of daily meteorological data equal or greater than time_max"))
+ try(if(METEO[ID %in% PARCELLE$STATION, .N, by = ID] %>% .[,N] %>% min %>% is_less_than(time_max))
+ stop("All stations must be associated to a number of daily METEOrological data equal or greater than time_max"))
## CALCULATION
- # Add diff_alt
- PARCELLE$diff_alt <- PARCELLE$ALTITUDE_UNIT - PARCELLE$ALTITUDE_STATION
-
## for each administrative unit
+ message("## Generating parameters for all patches can take time, please be patient. ##")
#make it parallel
+
clust <- detectCores() %>% makeCluster
#export variables
- clusterExport(clust, list("PARCELLE", "meteo", "time_max", "nYR_init"), envir = environment())
+ clusterExport(clust, list("PARCELLE", "METEO", "time_max", "nYR_init"), envir = environment())
#export required library
clusterEvalQ(clust, library("data.table"))
clusterEvalQ(clust, library("magrittr"))
clusterEvalQ(clust, library("zoo"))
+ clusterEvalQ(clust, library("pbapply"))
+
+ldata <- pbapply(PARCELLE, 1, function(x){
-ldata <- parApply(clust, PARCELLE, 1, function(x){
# Daily temperature corrected by the altitude
- temperature <- (meteo[POSTE == x["STATION"] %>% as.numeric, TP] - as.numeric(x["diff_alt"])*(4.2/1000)) %>% .[seq(time_max)]
+ temperature <- (METEO[ID == x["STATION"] %>% as.numeric, TP] - as.numeric(x["DIFF_ALT"])*(4.2/1000)) %>% .[seq(time_max)]
# Cumulative rainfall over 7 days
- raincumul7 <- meteo[POSTE == x["STATION"] %>% as.numeric, rollapply(c(rep(0,6), RR), 7, sum)][seq(time_max)]
+ raincumul7 <- METEO[ID == x["STATION"] %>% as.numeric, rollapply(c(rep(0,6), RR), 7, sum)][seq(time_max)]
kLfix <- x["KLfix"] %>% as.numeric
kLvar <- x["KLvar"] %>% as.numeric
@@ -98,13 +109,17 @@ ldata <- parApply(clust, PARCELLE, 1, function(x){
# Carrying capacities related to rainfall
kL_init, kL,
kP_init, kP
- )
+ ) %>% as.numeric
-})
+}, cl = clust)
stopCluster(clust)
-parse(text=(paste0("colnames(ldata) <- PARCELLE$",nodeID)))
+colnames(ldata) <- PARCELLE$ID
+
+ldata %<>% rbind(., mMov)
+
+rownames(ldata) <- NULL
return(ldata)
}
diff --git a/R/generate_events.R b/R/generate_events.R
index 5b97e5ae26a6a38dc6005420f8defef86621c898..21dd7239d5908111a3719361f6730f83530edd06 100644
--- a/R/generate_events.R
+++ b/R/generate_events.R
@@ -15,7 +15,7 @@
#'
#' @export
-build_E <- function(timetosample, nintro, compartments, nodeintro){
+build_E <- function(timetosample, nintro, nind, compartments, nodeintro){
# E Matrix for the select process
E <- structure(.Data = c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, #1
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, #2
@@ -45,7 +45,7 @@ events <- data.frame(
time = sample(timetosample, nintro, replace = T), ## The time that the event happens
node = sample(nodeintro, nintro, replace = T), ## In which node does the event occur
dest = rep(0, nintro), ## Which node is the destination node
- n = rep(1, nintro), ## How many individuals are moved
+ n = rep(nind, nintro), ## How many individuals are moved
proportion = rep(0, nintro), ## This is not used when n > 0
select = rep(2, nintro), ## Use the 4th column in the model select matrix
shift = rep(0, nintro) ## Not used in this example
diff --git a/R/iniState.R b/R/iniState.R
index 41291589aa80b3990953fe4fbd033a6ec8f27d3e..f25c931239817af6b0a38bf95fda61f552026852 100644
--- a/R/iniState.R
+++ b/R/iniState.R
@@ -1,11 +1,9 @@
#' Write initial state of the meta-population
#'
#' @description initialize
-#' @usage iniState(PARCELLE, nodeID, popID)
+#' @usage iniState(PARCELLE)
#'
#' @param PARCELLE data.frame or data.table
-#' @param nodeID string or data.table
-#' @param popID string
#' @param initMosq numerical value Initial number of eggs in each node.
#'
#' @return list
@@ -15,9 +13,9 @@
#' @export
-iniState <- function(PARCELLE, nodeID, popID, initMosq = 100000){
+iniState <- function(PARCELLE, initMosq = 100000){
- nh <- PARCELLE[, ..popID] %>% round
+ nh <- PARCELLE[, POP] %>% round
nadm <- nrow(PARCELLE)
u0 <- data.frame(
@@ -52,13 +50,13 @@ cont_var <- c(
"z", "temperature", "kL", "kP",
"Em", "Lm", "Pm",
"Aemm", "A1hm", "A1gm", "A1om", "A2hm", "A2gm", "A2om", "prevEggs",
- "R0")
+ "R0", "pIh", "pIm")
v0 <- matrix(0,
nrow = length(cont_var),
- ncol = PARCELLE[, ..nodeID] %>% nrow,
+ ncol = PARCELLE[, ID] %>% length,
dimnames=list(cont_var,
- unlist(PARCELLE[, ..nodeID])
+ unlist(PARCELLE[, ID])
))
## Random initial mosquito population size (eggs)
diff --git a/R/runArboRisk.R b/R/runArboRisk.R
new file mode 100644
index 0000000000000000000000000000000000000000..b26dc2338fb1c350d6370005d87ec4f0c5f834db
--- /dev/null
+++ b/R/runArboRisk.R
@@ -0,0 +1,266 @@
+#' Run a simulation
+#'
+#' @description function to run simulations
+#'
+#' @usage runArboRisk(PARCELLE, METEO, gdata, mMov = NULL, start_date, end_date, nYR_init = 2, nodeID = "ID")
+#'
+#' @param PARCELLE data.frame or data.table
+#' @param METEO data.frame or data.table
+#' @param gdata list
+#' @param mMov double matrix
+#' @param start_date date in "%Y-%m-%d" format
+#' @param end_date date in "%Y-%m-%d" format
+#' @param nYR_init numerical value. Number of years used to initialize the population dynamics (default is 1)
+#'
+#' @importFrom SimInf mparse
+#' @importFrom SimInf run
+#' @importFrom SimInf trajectory
+#'
+#' @return data.table
+#'
+#' @export
+
+
+runArboRisk <- function(PARCELLE,
+ METEO,
+ gdata = NULL,
+ vector_species = "Ae. albopictus",
+ mMov = NULL,
+ start_date = NULL,
+ end_date = NULL,
+ nYR_init = 1,
+ diapause = TRUE){
+
+
+ ####################
+ ### General data ###
+ ####################
+ # input data:
+ # PARCELLE : table where each row is a node and columns are attributes
+ # METEO: table with four columns nodeID, date, daily rainfall, daily mean temperature
+
+ if(!(is.data.frame(PARCELLE) | is.data.table(PARCELLE) | inherits(PARCELLE, "SpatVector")))
+ stop("PARCELLE must be a data.table or a data.frame")
+
+ if(inherits(PARCELLE, "SpatVector")){
+ shape <- PARCELLE
+ PARCELLE %<>% as.data.frame
+ }
+
+ PARCELLE <- copy(PARCELLE)
+ setDT(PARCELLE)
+
+ if(NA %in% match(c("ID", "POP", "SURF_HA", "KLfix", "KLvar", "KPfix", "KPvar", "STATION", "DIFF_ALT"), names(PARCELLE)))
+ stop('PARCELLE must contain at least 16 columns: "ID", "POP", "SURF_HA", "KLfix", "KLvar", "KPfix", "KPvar", "STATION", "DIFF_ALT"')
+
+ ## Inform user on number of PARCELLES, include possibility to use shape
+ ## check PARCELLE columns
+ # popID = "POP"
+
+ if(!(is.data.frame(METEO) | is.data.table(METEO)))
+ stop("METEO must be a data.table or a data.frame")
+
+ METEO <- copy(METEO)
+ setDT(METEO)
+
+ if(NA %in% match(c("ID", "DATE", "RR", "TP"), names(METEO)))
+ stop("METEO must contain at least 4 columns: 'ÃŒD', 'DATE', 'RR', 'TP' ")
+
+
+ ## FIX ME check if meteo contains daily records and print period
+ ## check METEO columns
+ # rename POSTE
+ # DATE "%Y-%m-%d"
+ # RR double normalized between 0 and 1
+ # TP double < 70
+
+ ############################
+ ### Model initialization ###
+ ############################
+
+ if(is.null(start_date))
+ start_date <- min(METEO$DATE)
+
+ if(is.null(end_date))
+ end_date <- max(METEO$DATE)
+
+ if(start_date < min(METEO$DATE) | end_date > max(METEO$DATE))
+ stop("start_date and end_date must be included in meteorological records")
+
+ #### Simulation duration ####
+ time_max <- seq(from = min(METEO$DATE),
+ to = max(METEO$DATE), by = "days") %>% length
+
+ ###########################################
+ ## Create u0, v0 and compartment objects ##
+ ###########################################
+
+ message("Initialization")
+
+ init <- iniState(PARCELLE, initMosq = 100000)
+
+ u0 <- init$u0
+ v0 <- init$v0
+
+ ###############################################################
+ ## Define the compartments, gdata and stochastic transitions ##
+ ###############################################################
+ ## Stochastic transitions are related to the epidemiological model and the life cycle of infected mosquitoes.
+
+ #### Global data (general parameters) ####
+
+ if(is.null(gdata)){
+ if(is.null(vector_species))
+ stop("If vector_species is NULL, a set of global parameters gdata must be provided")
+
+ if(vector_species == "Ae. albopictus"){
+ gdata <- c(
+ # humans
+ bMH = 0.5, # Infection probability from vector to host - 0.5
+ delta = 4.5,
+ muH = 1/5, # transition rate from exposed (E) to infected (I) for humans (1/days) (doi: 10.1038/s41598-019-53127-z : 1/2)
+ rhoH = 1/5, # Recovery rate of humans (1/days) (doi: 10.1038/s41598-019-53127-z : 1/4)
+ # eggs
+ muE = 0.05,
+ TE = 10, # 10.4 in Tran 2013
+ TDDE = 110,
+ # larva
+ # muL = 0.08,
+ q1L = -0.0007,
+ q2L = 0.0392,
+ q3L = -0.3911,
+ # pupa
+ # muP = 0.03,
+ muEM = 0.1,
+ q1P = 0.0008,
+ q2P = -0.0051,
+ q3P = 0.0319,
+ # emerging adults
+ muA = 0.025, # 0.02 dans le papier // 0.025 dans ocelet
+ gammaAem = 0.4,
+ sigma = 0.5, # proportion of females
+ # host-seeking adults
+ gammaAh = 0.3, # 0.2 dans Tran 2013
+ muR = 0.08,
+ # engorged adults
+ TAG = 10,
+ TDDAG = 77,
+ # oviposition-site-seeking adults
+ gammaAo = 0.28, # fao == gammaAo = 0.2 in Tran 2013 ?
+ beta1 = 95,
+ beta2 = 75,
+ # mosquitoes infection
+ bHM = 0.31,
+ startFav = "10/03/2020",
+ endFav = "30/09/2020"
+ )
+ } else
+ stop(paste("Parameters for", vector_species, "haven't been implemented, please provide a set of gdata parmaeters or choose 'Ae. albopictus' vector_species"))
+ }
+
+ if(NA %in% match(c("bMH","delta","muH","rhoH","muE","TE","TDDE","q1L","q2L","q3L","muEM","q1P","q2P","q3P","muA","gammaAem","sigma","gammaAh","muR","TAG","TDDAG","gammaAo","beta1","beta2","bHM"), names(gdata)))
+ stop('gdata must contain at least the following parameters:
+ "bMH","delta","muH","rhoH","muE","TE","TDDE","q1L","q2L","q3L","muEM","q1P","q2P","q3P","muA","gammaAem","sigma","gammaAh","muR","TAG","TDDAG","gammaAo","beta1","beta2","bHM"')
+
+ ##############
+ ## Diapause ##
+ ##############
+
+ # Calculates days for diapause over all the simulated period and build pts_fun
+ if(isTRUE(diapause)){
+ diapause_interv <- diapause(startFav = gdata["startFav"],
+ endFav = gdata["endFav"],
+ time_max,
+ nYR_init)
+ gdata %<>% .[-which(names(gdata) %in% c("startFav","endFav"))]
+ }
+
+ mode(gdata) <- "numeric"
+
+
+ message("Build stochastic transitions")
+ #### Stochastic transitions ####
+ stoch_transitions <- build_transitions()
+
+ ##################
+ ## Define ldata ##
+ ##################
+ ## ldata is a vector of local data: daily temperature, kL and kP and number of non infected mosquitoes in each stage.
+ ## the first item (index = 0) is the number of simulated days, it will be used to daily update v0 with ldata content.
+ ## Now improve this and use time to index data in ldata (local data).
+
+ message("Format local data")
+
+ if(is.null(mMov)){
+ mMov <- diag(nrow(PARCELLE))
+ rownames(mMov) <- PARCELLE[, ID]
+ }
+
+ ldata <- create_ldata(PARCELLE, METEO, gdata, mMov, time_max, nYR_init)
+
+ ###########################################################
+ ### Deterministic transitions for mosquitoes life cycle ###
+ ###########################################################
+ # pts_fun is a post time step function in C_code updating v0 at each time step.
+ # It is used for meterological and evironemental data, diapause trigger and
+ # vector population stages driven by deterministic events.
+
+ pts_fun <- build_pts_fun(u0, v0, gdata, ldata, diapause_interv)
+
+ ###########################################
+ ### Introduction of exposed individuals ###
+ ###########################################
+
+ # mov <- build_E(timetosample = diapause_interv,
+ # nintro = 10,
+ # nind = 100,
+ # u0 %>% colnames,
+ # nodeintro = 1)
+ #
+ # E <- mov$E
+ # N <- mov$N
+ # events <- mov$events
+
+ ###################
+ ### Simulations ###
+ ###################
+
+ # u0[1, 2] <- 1000
+ # u0[, 7] <- 1000
+
+ # start_time <- Sys.time()
+
+ message("Build and run the model")
+ message(Sys.time())
+
+ model <- mparse(
+ transitions = stoch_transitions,
+ compartments = u0 %>% colnames,
+ gdata = gdata,
+ ldata = ldata,
+ u0 = u0,
+ v0 = v0,
+ tspan = seq(from = 0, to = (time_max + nYR_init*365) - 1, by = 1), # number of days simulated. Start at 0 to get the indexing correct in the post-time-step function.
+ pts_fun = pts_fun#,
+ # introduction of exposed individuals
+ # events = events,
+ # E = E,
+ # N = N
+ )
+
+ result <- run(model = model)
+
+ traj <- trajectory(result)
+ # end_time <- Sys.time()
+ # end_time - start_time
+
+ setDT(traj)
+
+ traj[, DATE := seq.Date(from = (start_date - (365*nYR_init)),
+ to = end_date,
+ by ="day"), by = node]
+
+ traj[, ID := PARCELLE$ID[node]]
+
+ return(traj)
+}
diff --git a/data/.gitignore b/data/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..b1e54437876c9b8433e283a46992e87c816bf1b5
--- /dev/null
+++ b/data/.gitignore
@@ -0,0 +1,2 @@
+PARCELLE.rda
+meteo.rda
diff --git a/data/PARCELLE.rda b/data/PARCELLE.rda
index 918e03c432f1c0d106dd86583f4b0e844594f030..96997417e8f6863f471a223f80e44bd4cdee0dec 100644
Binary files a/data/PARCELLE.rda and b/data/PARCELLE.rda differ
diff --git a/data/meteo.rda b/data/meteo.rda
index c4fbc5e444e4148c3e523955a13b16b560fe4645..7d2941e55968a2dad41da7a3688a67fbc9b12699 100644
Binary files a/data/meteo.rda and b/data/meteo.rda differ
diff --git a/man/create_ldata.Rd b/man/create_ldata.Rd
index fa3b327043c78e4fb58aac36ed20dceac5f1ba97..b99cd1149a02d6634d4fe811922320ae83e9a480 100644
--- a/man/create_ldata.Rd
+++ b/man/create_ldata.Rd
@@ -4,12 +4,12 @@
\alias{create_ldata}
\title{Write local data matrix}
\usage{
-create_ldata(PARCELLE, meteo, gdata, time_max, nYR_init, nodeID)
+create_ldata(PARCELLE, METEO, gdata, time_max, nYR_init, nodeID)
}
\arguments{
\item{PARCELLE}{data.frame or data.table}
-\item{meteo}{data.frame or data.table}
+\item{METEO}{data.frame or data.table}
\item{gdata}{list}
@@ -23,7 +23,7 @@ create_ldata(PARCELLE, meteo, gdata, time_max, nYR_init, nodeID)
matrix
}
\description{
-Internal function to write the pts_fun code driving SimInf package to implement deterministic population dynamics of Aedes Albopictus.
+Function to define local data
}
+\keyword{METEO}
\keyword{ldata}
-\keyword{meteo}
diff --git a/man/iniState.Rd b/man/iniState.Rd
index 770f61052952d348764197d4c6f0b6a5db915b65..f852b9a56a33927bdae477747ab0e043cf7e4533 100644
--- a/man/iniState.Rd
+++ b/man/iniState.Rd
@@ -4,15 +4,11 @@
\alias{iniState}
\title{Write initial state of the meta-population}
\usage{
-iniState(PARCELLE, nodeID, popID)
+iniState(PARCELLE)
}
\arguments{
\item{PARCELLE}{data.frame or data.table}
-\item{nodeID}{string or data.table}
-
-\item{popID}{string}
-
\item{initMosq}{numerical value Initial number of eggs in each node.}
}
\value{
diff --git a/man/runArboRisk.Rd b/man/runArboRisk.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..1a2282ba5d8735952f9395cc07fa72cc1c22adb4
--- /dev/null
+++ b/man/runArboRisk.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/runArboRisk.R
+\name{runArboRisk}
+\alias{runArboRisk}
+\title{Run a simulation}
+\usage{
+runArboRisk(PARCELLE, METEO, gdata, mMov = NULL, start_date, end_date, nYR_init = 2, nodeID = "ID")
+}
+\arguments{
+\item{PARCELLE}{data.frame or data.table}
+
+\item{METEO}{data.frame or data.table}
+
+\item{gdata}{list}
+
+\item{mMov}{double matrix}
+
+\item{start_date}{date in "%Y-%m-%d" format}
+
+\item{end_date}{date in "%Y-%m-%d" format}
+
+\item{nYR_init}{numerical value. Number of years used to initialize the population dynamics (default is 1)}
+}
+\value{
+data.table
+}
+\description{
+function to run simulations
+}
diff --git a/work_in_progress/.gitignore b/work_in_progress/.gitignore
index d7f2e836aeaf086129fd7b6c1ceee6e0bf75e439..6e55a39499fd61f6087c9fc754eed2cc2c708889 100644
--- a/work_in_progress/.gitignore
+++ b/work_in_progress/.gitignore
@@ -8,3 +8,4 @@ u0.Rda
v0.Rda
data
PARCELLE.rda
+METEO.rda
diff --git a/work_in_progress/contact_method.R b/work_in_progress/contact_method.R
new file mode 100644
index 0000000000000000000000000000000000000000..eb7ab9f63f000410d5155704b966f9b025d65942
--- /dev/null
+++ b/work_in_progress/contact_method.R
@@ -0,0 +1,80 @@
+# Create contact
+
+## Levels of contact
+
+## Probabilities of contact
+# Very high within patch (patch level)
+Hp <- 0.95
+# Medium (eg. between patches of a city ; city level)
+Mp <- 0.25
+# Low (eg. between patches of a department ; department level)
+Lp <- 0.05
+# Very low (eg. between patches of different departments ; global level)
+Gp <- 0.01
+
+library(ohenery)
+normalize(c(Hp, Mp, Lp, Gp))
+
+# High within patch (patch level)
+Hp <- 0.59375
+# Medium between patches of a city (city level)
+Mp <- 0.28125
+# Low between patches of a department (department level)
+Lp <- 0.09375
+# Very low between patches of different departments (national level)
+Gp <- 0.03125
+
+PARCELLE[, SL2 := NOM_COM %>% factor %>% as.numeric]
+PARCELLE[, SL3 := INSEE_COM %>% substr(., 1, 2) %>% factor %>% as.numeric]
+save(PARCELLE, file = "D:/Pachka_tmp/1-Modeling/dengue-risk-assessment/data/PARCELLE.rda")
+save(PARCELLE, file = "D:/Pachka_tmp/1-Modeling/dengue-risk-assessment/work_in_progress/PARCELLE.rda")
+
+
+int IP;
+IP = u[3];
+
+int TP;
+TP = pop;
+
+int IC;
+int TC;
+if(ldata[1] == ){
+ IC = u[...] + u[...] + ...;
+ TC = u[...] + u[...] + ...;
+}
+
+
+
+ISL2 <- c()
+for(SpUnit in PARCELLE$SL2 %>% unique){
+ ISL2 %<>% append(., paste0(c("x = if(ldata[1] == ", SpUnit ,"){ IC = ", paste0("u[", 2 + which(PARCELLE$SL2 == SpUnit) - 1,"]", collapse = " + "), "; TC = ", paste0("u[", 1 + which(PARCELLE$SL2 == SpUnit) - 1,"] + u[", 2 + which(PARCELLE$SL2 == SpUnit) - 1,"] + u[",3 + which(PARCELLE$SL2 == SpUnit) - 1,"] + u[", 4 + which(PARCELLE$SL2 == SpUnit) - 1,"]", collapse = " + "), ";}"), collapse = ""))
+ }
+
+
+ISL3 <- c()
+for(SpUnit in PARCELLE$SL3 %>% unique){
+ ISL3 %<>% append(., paste0(c("x = if(ldata[2] == ", SpUnit ,"){ IC = ", paste0("u[", 2 + which(PARCELLE$SL3 == SpUnit) - 1,"]", collapse = " + "), "; TC = ", paste0("u[", 1 + which(PARCELLE$SL3 == SpUnit) - 1,"] + u[", 2 + which(PARCELLE$SL3 == SpUnit) - 1,"] + u[",3 + which(PARCELLE$SL3 == SpUnit) - 1,"] + u[", 4 + which(PARCELLE$SL3 == SpUnit) - 1,"]", collapse = " + "), ";}"), collapse = ""))
+}
+
+int ID;
+int TD;
+if(ldata[2] == ){
+ ID = u[...] + u[...] + ... ;
+ TD = u[...] + u[...] + ... ;
+}
+
+double pI;
+pI = PLp * (IP/TP) +
+ CLp * (sum(IC) - IP)/(sum(TC)-TP) +
+ DLp * (sum(ID) - IC)/sum(TD))-TC) +
+ NLp * (sum(IN) - ID)/sum(TN)-TD)
+
+
+## >>> cummulative probabilities for the patch
+
+# Medium between patches of a city (city level)
+CLp <- 0.65
+# Low between patches of a department (department level)
+DLp <- 0.35
+# Very low between patches of different departments (national level)
+NLp <- 0.05
diff --git a/work_in_progress/main.R b/work_in_progress/main.R
index 7abc7233757852989675660ffeaaf9aa96e80217..a31aa8e35fd1d1efcaf7b633a92282418f5da6f9 100644
--- a/work_in_progress/main.R
+++ b/work_in_progress/main.R
@@ -14,6 +14,7 @@ gc() # free up memory and report the memory usage.
###################################
### Load packages and functions ###
###################################
+devtools::document()
devtools::install(quick = T)
library(ArboRisk)
@@ -21,18 +22,6 @@ library(ArboRisk)
library(data.table)
library(magrittr)
-library(ArboRisk)
-library(SimInf)
-
-# convert time as date in the output
-library(lubridate)
-library(hydroTSM)
-
-# Parallel processing
-library(parallel)
-library(doParallel)
-library(foreach)
-
####################
### General data ###
####################
@@ -41,116 +30,9 @@ library(foreach)
# meteo: table with four columns nodeID, date, daily rainfall, daily mean temperature
PARCELLE
-PARCELLE %<>% .[startsWith(PARCELLE$INSEE_COM, "34"),]
-meteo
-
-############################
-### Model initialization ###
-############################
-
-#### Simulation duration ####
-time_max <- seq(from = as.Date(as.character(min(meteo$DATE)), format="%Y%m%d"),
- to = as.Date(as.character(max(meteo$DATE)), format="%Y%m%d"), by = "days") %>% length
-# time_max <- 400 # number of days simulated
-nYR_init <- 2 # number of years used to initialize the population dynamics
-
-###########################################
-## Create u0, v0 and compartment objects ##
-###########################################
-
-init <- iniState(PARCELLE, nodeID = "INSEE_COM", popID = "POP", initMosq = 100000)
-
-u0 <- init$u0
-v0 <- init$v0
-
-###############################################################
-## Define the compartments, gdata and stochastic transitions ##
-###############################################################
-## Stochastic transitions are related to the epidemiological model and the life cycle of infected mosquitoes.
-
-#### Global data (general parameters) ####
-gdata <- c(
- # humans
- bMH = 0.5, # Infection probability from vector to host - 0.5
- delta = 4.5,
- muH = 1/5, # transition rate from exposed (E) to infected (I) for humans (1/days) (doi: 10.1038/s41598-019-53127-z : 1/2)
- rhoH = 1/5, # Recovery rate of humans (1/days) (doi: 10.1038/s41598-019-53127-z : 1/4)
- # eggs
- muE = 0.05,
- TE = 10, # 10.4 in Tran 2013
- TDDE = 110,
- # larva
- # muL = 0.08,
- q1L = -0.0007,
- q2L = 0.0392,
- q3L = -0.3911,
- # pupa
- # muP = 0.03,
- muEM = 0.1,
- q1P = 0.0008,
- q2P = -0.0051,
- q3P = 0.0319,
- # emerging adults
- muA = 0.025, # 0.02 dans le papier // 0.025 dans ocelet
- gammaAem = 0.4,
- sigma = 0.5, # proportion of females
- # host-seeking adults
- gammaAh = 0.3, # 0.2 dans Tran 2013
- muR = 0.08,
- # engorged adults
- TAG = 10,
- TDDAG = 77,
- # oviposition-site-seeking adults
- gammaAo = 0.28, # fao == gammaAo = 0.2 in Tran 2013 ?
- beta1 = 95,
- beta2 = 75,
- # mosquitoes infection
- bHM = 0.31
-)
-
-#### Stochastic transitions ####
-stoch_transitions <- build_transitions()
-
-##################
-## Define ldata ##
-##################
-## ldata is a vector of local data: daily temperature, kL and kP and number of non infected mosquitoes in each stage.
-## the first item (index = 0) is the number of simulated days, it will be used to daily update v0 with ldata content.
-## Now improve this and use time to index data in ldata (local data).
-
-ldata <- create_ldata(PARCELLE, meteo, gdata, time_max, nYR_init = 2, nodeID = "INSEE_COM")
-
-##############
-## Diapause ##
-##############
-
-# Calculates days for diapause over all the simulated period and build pts_fun
-diapause_interv <- diapause(startFav = "10/03/2020",
- endFav = "30/09/2020",
- time_max,
- nYR_init)
-
-###########################################################
-### Deterministic transitions for mosquitoes life cycle ###
-###########################################################
-# pts_fun is a post time step function in C_code updating v0 at each time step.
-# It is used for meterological and evironemental data, diapause trigger and
-# vector population stages driven by deterministic events.
-
-pts_fun <- build_pts_fun(u0, v0, gdata, diapause_interv)
-
-###########################################
-### Introduction of exposed individuals ###
-###########################################
-
-mov <- build_E(timetosample = diapause_interv,
- nintro = 10,
- u0 %>% colnames,
- nodeintro = nrow(PARCELLE))
-
-E <- mov$E
-N <- mov$N
-events <- mov$events
+PARCELLE %<>% .[startsWith(PARCELLE$ID, "34"),]
+# PARCELLE %<>% .[1:100,]
+METEO
###################
### Simulations ###
@@ -158,24 +40,15 @@ events <- mov$events
start_time <- Sys.time()
-model <- mparse(
- transitions = stoch_transitions,
- compartments = u0 %>% colnames,
- gdata = gdata,
- ldata = ldata,
- u0 = u0,
- v0 = v0,
- tspan = seq(from = 0, to = (time_max + nYR_init*365) - 1, by = 1), # number of days simulated. Start at 0 to get the indexing correct in the post-time-step function.
- pts_fun = pts_fun#,
- # introduction of exposed individuals
- # events = events,
- # E = E,
- # N = N
-)
-
-result <- run(model = model)
-
-traj <- trajectory(result)
+traj <- runArboRisk(PARCELLE,
+ METEO,
+ gdata = NULL,
+ vector_species = "Ae. albopictus",
+ mMov = NULL,
+ start_date = NULL,
+ end_date = NULL,
+ nYR_init = 1,
+ diapause = TRUE)
end_time <- Sys.time()
end_time - start_time
@@ -184,7 +57,6 @@ end_time - start_time
##### Explore output
#####
-traj %>% setDT
traj
t2d <- rev(ymd('2022/11/30') - days(0:max(traj$time)))
diff --git a/work_in_progress/processInput.R b/work_in_progress/processInput.R
index f27f6549b7b005a70739b490c5531c78d555ad98..be5f63f5a9de12cd90cbde2b55e0743c6c18cf79 100644
--- a/work_in_progress/processInput.R
+++ b/work_in_progress/processInput.R
@@ -2,21 +2,21 @@ library(terra)
library(geodata)
#### Meteorological data ####
-meteo <- read.csv2("test/data/meteo_20210101_20221130.txt")
-meteo %>% setDT
-meteo %<>% cleanmeteo
+METEO <- read.csv2("work_in_progress/data/METEO_20210101_20221130.txt")
+METEO %>% setDT
+METEO %<>% cleanMETEO
-meteo[, `:=`(TN = as.numeric(TN),
+METEO[, `:=`(TN = as.numeric(TN),
TX = as.numeric(TX))]
-meteo[, TP := (TN+TX)/2]
-meteo[, `:=`(TN = NULL,
+METEO[, TP := (TN+TX)/2]
+METEO[, `:=`(TN = NULL,
TX = NULL)]
### List stations ###
-filename <- "test/data/stations.txt"
+filename <- "work_in_progress/data/stations.txt"
STATIONS <- read.csv2(filename, sep = ";", dec = ".")
setDT(STATIONS)
-STATIONS %<>% .[INSEE %in% meteo$POSTE] # Keep stations for which we have meteorological data
+STATIONS %<>% .[INSEE %in% METEO$ID] # Keep stations for which we have METEOrological data
#remove unsed columns
STATIONS[, `:=`(TYPE = NULL,
NOM = NULL,
@@ -25,7 +25,7 @@ STATIONS[, `:=`(TYPE = NULL,
STATIONS[, ALTITUDE := substr(ALTITUDE, 1, 4) %>% as.numeric]
#### List of administrative unit ####
-filename <- "test/data/IRIS-GE/IRIS_GE.SHP"
+filename <- "work_in_progress/data/IRIS-GE/IRIS_GE.SHP"
PARCELLE <- vect(filename)
terra::crs(PARCELLE, describe = T)
### Test limit the number of parcelles
@@ -47,7 +47,7 @@ PARCELLE$KPvar <- 50 * PARCELLE$SURF_HA
## Average altitude ##
# download raster
-elevation <- elevation_30s(country="FRA", path = "test/data/")
+elevation <- elevation_30s(country="FRA", path = "work_in_progress/data/")
# project raster
elevation %<>% project(., "epsg:2154")
# compute average
@@ -62,7 +62,7 @@ PARCELLE[is.na(PARCELLE$ALTITUDE_UNIT)]$ALTITUDE_UNIT <- 0
## Population ##
# download raster
-pop <- geodata::population(2020, 0.5, path = "test/data/")
+pop <- geodata::population(2020, 0.5, path = "work_in_progress/data/")
# project raster
pop %<>% terra::project(., elevation)
# compute average
@@ -72,7 +72,7 @@ PARCELLE$POP %>% is.na %>% sum %>% equals(0)
# Turn density into number of persons
PARCELLE$POP %<>% multiply_by(PARCELLE$SURF_HA/100) %>% round
-## Unique meteorological station ##
+## Unique METEOrological station ##
# two options: only spatial considerations or the weighted by the human population
# look at the function nearest from terra package
@@ -85,9 +85,11 @@ STATIONS_pts %<>% terra::project(., elevation)
PARCELLE$STATION <- STATIONS[nearest(PARCELLE, STATIONS_pts) %>% .$to_id, INSEE]
PARCELLE$ALTITUDE_STATION <- STATIONS[nearest(PARCELLE, STATIONS_pts) %>% .$to_id, ALTITUDE]
-## add the difference of altitude between the average altitude of the administrative unit and the altitude of the meteorological station (diff_alt = ALTITUDE_UNIT - ALTITUDE_STATION)
+## add the difference of altitude between the average altitude of the administrative unit and the altitude of the METEOrological station (diff_alt = ALTITUDE_UNIT - ALTITUDE_STATION)
PARCELLE$diff_alt <- PARCELLE$ALTITUDE_UNIT - PARCELLE$ALTITUDE_STATION
+names(PARCELLE)[4] <- "ID"
+
# turn PARCELLE into table
PARCELLE %<>% as.data.frame
setDT(PARCELLE)