diff --git a/DESCRIPTION b/DESCRIPTION
index 0ceb01a04f7f4c48222424c5129da7ef757bff36..216a1544620352a92d83143de5588aed1a8400cd 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -12,29 +12,34 @@ 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),
+ igraph,
+ fields
Suggests:
- knitr,
- rmarkdown
+ knitr (>= 1.42),
+ rmarkdown (>= 2.20),
+ testthat (>= 3.0.0)
VignetteBuilder:
knitr
RoxygenNote: 7.2.3
+Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
index 4bf2512649a66452a3b203aeb8f7be576729a696..0e970517d4c247f54b96dcb0623e7741dc921f33 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,11 +1,30 @@
# Generated by roxygen2: do not edit by hand
-export(build_E)
+export(build_E_random)
+export(build_gdata)
+export(build_ldata)
+export(build_mMov)
export(build_pts_fun)
export(build_transitions)
export(cleanmeteo)
export(create_ldata)
export(diapause)
export(iniState)
+export(runArboRisk)
import(data.table)
+import(igraph)
+importFrom(SimInf,mparse)
+importFrom(SimInf,run)
+importFrom(SimInf,trajectory)
+importFrom(fields,image.plot)
+importFrom(igraph,as_adjacency_matrix)
+importFrom(igraph,erdos.renyi.game)
importFrom(imputeTS,na_ma)
+importFrom(magrittr,"%>%")
+importFrom(parallel,clusterEvalQ)
+importFrom(parallel,clusterExport)
+importFrom(parallel,detectCores)
+importFrom(parallel,makeCluster)
+importFrom(parallel,stopCluster)
+importFrom(pbapply,pbapply)
+importFrom(pbapply,pblapply)
diff --git a/R/build_E_random.R b/R/build_E_random.R
new file mode 100644
index 0000000000000000000000000000000000000000..5d34d8f2df7dc0ee7526a377c2d9bf7652a91a39
--- /dev/null
+++ b/R/build_E_random.R
@@ -0,0 +1,46 @@
+#' Generate the list of introduction event
+#'
+#' @description Function used to generate the event matrix and associated objects required by SimInf package
+#'
+#' @usage build_E_random(period_start, period_end, n_ind = 1, n_events = 1, stage = "Eh", loc)
+#'
+#' @param period_start date in '\%Y-\%m-\%d' format. Define the beginning of the period for introduction.
+#' @param period_end date in '\%Y-\%m-\%d' format. Define the end of the period for introduction.
+#' @param n_ind Numerical vector. Range of number of individuals introduced.
+#' @param n_events String vector. Number of introduction generated.
+#' @param stage String vector. Epidemiological stage of introduced individuals.
+#' @param loc String or numerical vector. List of patch to seed random introductions.
+#'
+#'
+#' @return events matrix with scheduled introductions.
+#'
+#' @keywords events
+#'
+#' @examples
+#' build_E_random(period_start = as.Date("2020/03/10"), period_end = as.Date("2020/09/30"), n_ind = 1, n_events = 10, stage = "Eh", loc = LETTERS[1:3])
+#'
+#' @export
+
+build_E_random <- function( period_start,
+ period_end,
+ n_ind = 1,
+ n_events = 1,
+ stage = "Eh",
+ loc){
+
+ period <- seq(period_start, period_end, by = "day")
+
+ if(length(period) < n_events) stop("The number of events cannot exceed the number of days in the period")
+
+ sampling <- expand.grid(period, loc, n_ind)
+ sampling <- sampling[sample(seq(nrow(sampling)), n_events),]
+
+events <- data.frame(
+ time = sampling$Var1, ## The time that the event happens
+ node = sampling$Var2, ## In which node does the event occur
+ n = sampling$Var3, ## How many individuals are moved
+ select = stage ## Use the 4th column in the model select matrix
+ )
+
+return(events)
+}
diff --git a/R/build_gdata.R b/R/build_gdata.R
new file mode 100644
index 0000000000000000000000000000000000000000..655873c535b14c9e6575cb52f540303ed1987cbb
--- /dev/null
+++ b/R/build_gdata.R
@@ -0,0 +1,180 @@
+#' Generate global parameters dataset
+#'
+#' @description
+#' Then the functions formats all global parameters required for the transitions as a named list readable by SimInf package.
+#'
+#' @usage build_gdata()
+#'
+#' @importFrom magrittr %>%
+#'
+#' @param vector_species
+#'
+#' @return Named list of parameters
+#'
+#' @examples build_gdata()
+#'
+#' @export
+
+build_gdata <- function(
+ vector_species = "Ae. albopictus",
+ climate = "temperate",
+ ##### Infection
+ # HUMAN
+ # Infection probability from vector to host
+ bMH = NULL,
+ # delta = NULL,
+ muH = NULL, # transition rate from exposed (E) to infected (I) for humans (1/days) (doi: 10.1038/s41598-019-53127-z : 1/2)
+ rhoH = NULL, # Recovery rate of humans (1/days) (doi: 10.1038/s41598-019-53127-z : 1/4)
+ # eggs
+ muE = NULL,
+ TE = NULL,
+ TDDE = NULL,
+ # larva
+ # muL = 0.08,
+ q1L = NULL,
+ q2L = NULL,
+ q3L = NULL,
+ muEM = NULL,
+ q1P = NULL,
+ q2P = NULL,
+ q3P = NULL,
+ # emerging adults
+ muA = NULL,
+ gammaAem = NULL,
+ sigma = NULL,
+ # host-seeking adults
+ gammaAh = NULL,
+ muR = NULL,
+ # engorged adults
+ TAG = NULL,
+ TDDAG = NULL,
+ # oviposition-site-seeking adults
+ gammaAo = NULL,
+ beta1 = NULL,
+ beta2 = NULL,
+ # mosquitoes infection
+ bHM = NULL,
+ startFav = NULL,
+ endFav = NULL,
+ verbose = T) {
+
+ if(! vector_species %in% c(NULL, "Ae. albopictus"))
+ stop(
+ "Only Ae. albopictus in temperate climate has been implemented yet, nevertheless all parameters can be modified."
+ )
+
+ if(! climate %in% c(NULL, "temperate"))
+ stop(
+ "Only Ae. albopictus in temperate climate has been implemented yet, nevertheless all parameters can be modified."
+ )
+
+ if(vector_species == "Ae. albopictus" & climate == "temperate"){
+ ##### Infection
+ if(is.null(bMH)) bMH = 0.5
+ # if(is.null(delta)) delta = 4.5
+ if(is.null(muH)) muH = 1/5
+ if(is.null(rhoH)) rhoH = 1/5
+ if(is.null(muE)) muE = 0.05
+ if(is.null(TE)) TE = 10
+ if(is.null(TDDE)) TDDE = 110
+ if(is.null(q1L)) q1L = -0.0007
+ if(is.null(q2L)) q2L = 0.0392
+ if(is.null(q3L)) q3L = -0.3911
+ if(is.null(muEM)) muEM = 0.1
+ if(is.null(q1P)) q1P = 0.0008
+ if(is.null(q2P)) q2P = -0.0051
+ if(is.null(q3P)) q3P = 0.0319
+ if(is.null(muA)) muA = 0.025
+ if(is.null(gammaAem)) gammaAem = 0.4
+ if(is.null(sigma)) sigma = 0.5
+ if(is.null(gammaAh))gammaAh = 0.3
+ if(is.null(muR)) muR = 0.08
+ if(is.null(TAG)) TAG = 10
+ if(is.null(TDDAG)) TDDAG = 77
+ if(is.null(gammaAo)) gammaAo = 0.28
+ if(is.null(beta1)) beta1 = 95
+ if(is.null(beta2)) beta2 = 75
+ if(is.null(bHM)) bHM = 0.31
+ if(is.null(startFav)) startFav = as.Date("10/03/2020", format = "%d/%m/%Y")
+ if(is.null(endFav)) endFav = as.Date("30/09/2020", format = "%d/%m/%Y")
+ }
+
+ ### Checks
+
+ sapply(c(
+ bMH,
+ muH,
+ muE,
+ muEM,
+ muA,
+ gammaAem,
+ sigma,
+ gammaAh,
+ muR,
+ gammaAo,
+ bHM
+ ), function(x)
+ if (x < 0 | x > 1)
+ stop(
+ paste(x, " must fall between 0 and 1.")
+ ))
+
+ if(!is.null(startFav) & !is.null(startFav)){
+ if(class(startFav) != "Date")
+ startFav %<>% as.Date
+
+ if(class(endFav) != "Date")
+ endFav %<>% as.Date
+
+ if(verbose)
+ message("Diapause period is from ", startFav, " to ", endFav)
+ }
+
+ sapply(c(beta1, beta2), function(x)
+ if (x < 0)
+ stop(
+ "Number of eggs layed beta1 and beta2 must be positive numbers."
+ ))
+
+ gdata = list(
+ # humans
+ bMH = bMH, # Infection probability from vector to host - 0.5
+ # delta = delta,
+ muH = muH, # transition rate from exposed (E) to infected (I) for humans (1/days) (doi: 10.1038/s41598-019-53127-z : 1/2)
+ rhoH = rhoH, # Recovery rate of humans (1/days) (doi: 10.1038/s41598-019-53127-z : 1/4)
+ # eggs
+ muE = muE,
+ TE = TE, # 10.4 in Tran 2013
+ TDDE = TDDE,
+ # larva
+ q1L = q1L,
+ q2L = q2L,
+ q3L = q3L,
+ # pupa
+ muEM = muEM,
+ q1P = q1P,
+ q2P = q2P,
+ q3P = q3P,
+ # emerging adults
+ muA = muA, # 0.02 dans le papier // 0.025 dans ocelet
+ gammaAem = gammaAem,
+ sigma = sigma, # proportion of females
+ # host-seeking adults
+ gammaAh = gammaAh, # 0.2 dans Tran 2013
+ muR = muR,
+ # engorged adults
+ TAG = TAG,
+ TDDAG = TDDAG,
+ # oviposition-site-seeking adults
+ gammaAo = gammaAo,
+ beta1 = beta1,
+ beta2 = beta2,
+ # mosquitoes infection
+ bHM = bHM,
+ startFav = startFav,
+ endFav = endFav
+ )
+
+ return(gdata)
+
+}
diff --git a/R/build_ldata.R b/R/build_ldata.R
new file mode 100644
index 0000000000000000000000000000000000000000..f3b4228b5b232e2f8fd6707c561a643022911dd1
--- /dev/null
+++ b/R/build_ldata.R
@@ -0,0 +1,133 @@
+#' Build ldata
+#'
+#' @description function to build ldata
+#'
+#' @usage build_ldata(PARCELLE, METEO)
+#'
+#' @param PARCELLE data.frame or data.table
+#' @param METEO data.frame or data.table
+#' @param gdata list
+#' @param vector_species string
+#' @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
+#' @importFrom pbapply pblapply
+#'
+#' @return matrix
+#'
+#' @export
+
+
+build_ldata <- function(PARCELLE,
+ METEO,
+ gdata = NULL,
+ vector_species = "Ae. albopictus",
+ mMov = NULL,
+ start_date = NULL,
+ end_date = NULL,
+ nYR_init = 1){
+
+
+ ####################
+ ### 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) %>% as.Date
+
+ if(is.null(end_date))
+ end_date <- max(METEO$DATE) %>% as.Date
+
+ if(start_date < min(METEO$DATE) | end_date > max(METEO$DATE))
+ stop("start_date and end_date must be included in meteorological records")
+
+ message(paste("Simulation period is from",start_date , "to",end_date))
+
+ #### Simulation duration ####
+ time_max <- seq(from = start_date,
+ to = end_date, by = "days") %>% length
+
+ ###############################################################
+ ## 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 <- check_gdata(gdata, vector_species, verbose =F)
+
+ ##############
+ ## Diapause ##
+ ##############
+
+ # Calculates days for diapause over all the simulated period and build pts_fun
+ gdata %<>% .[-which(names(gdata) %in% c("startFav","endFav"))] %>% unlist
+
+ mode(gdata) <- "numeric"
+
+ ##################
+ ## 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).
+
+ if(is.null(mMov)){
+ mMov <- diag(nrow(PARCELLE))
+ rownames(mMov) <- PARCELLE[, ID]
+ colnames(mMov) <- PARCELLE[, ID]
+ message("The current scenario do not consider any human mobility (mMov = NULL)")
+ }
+
+ ldata <- create_ldata(PARCELLE, METEO, gdata, mMov, time_max, nYR_init)
+
+ ldata
+ }
diff --git a/R/build_mMov.R b/R/build_mMov.R
new file mode 100644
index 0000000000000000000000000000000000000000..25138511a8997bc9e64ff6b87fb5ea18b1233816
--- /dev/null
+++ b/R/build_mMov.R
@@ -0,0 +1,136 @@
+#' Creates a synthetic network
+#'
+#' @description This function creates a synthetic human mobility network.
+#' You can use this script to generate different input data to be used as examples of "network structures".
+#'
+#' @param PARCELLE data.table
+#' @param group string
+#' @param clust_ratio_inout numeric. Ratio of intra-building clustering (relative to inter-building clustering)
+#' @param within_clust_lev numeric. Probability of creating links between wards of the same building
+#' @param between_clust_lev numeric. Probability of creating links between wards of different building
+#' @param verbose logical. Print plot if TRUE
+#'
+#' @importFrom igraph erdos.renyi.game
+#' @importFrom igraph as_adjacency_matrix
+#' @importFrom fields image.plot
+#' @import igraph
+#'
+#' @usage build_mMov()
+#'
+#' @return a matrix with random probabilities of contact
+#'
+#' @export
+
+build_mMov <- function(PARCELLE,
+ group,
+ within_clust_lev = 0.8,
+ between_clust_lev = 0.1,
+ clust_ratio_inout = 0.8,
+ verbose = F){
+
+ #######################################################################
+
+ ## Checks
+ if(clust_ratio_inout < 0 | clust_ratio_inout > 1) stop("\'clust_ratio_inout\' is a ratio and must be between 0 and 1")
+
+ ## total number of patchs
+ tot_n <- nrow(PARCELLE)
+
+ ## matContact (contact matrix between wards)
+ ## this is the matrix of the (proportion of) time spent by HCW in a given ward (row) in any other ward in the hospital (columns)
+ ## sum of each row must be equal to 100 (i.e. 100%)
+
+ ###############################################################################
+ # Function that takes a graph as input and generates adjacency matrix as output
+ ###############################################################################
+
+ contact_matrix_generator <- function(g) {
+ adj_M <- as_adjacency_matrix(g,sparse = F) %>% as.matrix()
+ res <- apply(adj_M,1,function(x){
+ if(sum(x) == 0) normalized_p = x else normalized_p = x/sum(x)
+ normalized_p
+ })
+ return(res)
+ }
+
+ # The idea in order to build the synthetic hospital ward network is to assume that hospital wards are clustered
+ # (for example because they share the same building).
+ # Wards within a building(=cluster) will be highly connected, while connection across buildings will be less frequent.
+ # To implement this, we create a first layer with connection within buildings, a second layer with links across the buildings,
+ # and we sum the two (with weights)
+
+ ## the following function generates a matrix from the weighted sum of two networks
+ generate_network <- function(PARCELLE,
+ group,
+ tot_n,
+ p = c(within_clust_lev,between_clust_lev),
+ dist_within_between = c(clust_ratio_inout,(1-clust_ratio_inout))){
+
+
+ grouped <- PARCELLE %>% split(., by=group)
+
+ # vector indicating the parameters of the erdos renyi graph for within building and at the hospital level networks
+ # build a erdos.renyi.game network for each building
+ networks_list <- lapply(grouped, function(group){
+ clust_network <- erdos.renyi.game(nrow(group), p = p[1], type = "gnp")
+ V(clust_network)$name <- group$ID
+ clust_network
+ })
+
+ # build a erdos.renyi.game network for the whole hospital
+ networks_full <- erdos.renyi.game(tot_n, p = p[2], type = "gnp")
+ V(networks_full)$name <- PARCELLE$ID
+
+ # create contact matrix
+ M_full <- contact_matrix_generator(networks_full)
+
+ M <- list()
+ for(i in 1:length(networks_list)){
+ M[[i]] <- contact_matrix_generator(networks_list[[i]])
+ }
+
+ # create a matrix of size n_wards*n_buidings with the blocks only
+ M_cluster <- matrix(0, nrow=tot_n,ncol=tot_n)
+ rownames(M_cluster) <- colnames(M_cluster) <- PARCELLE$ID
+
+ for(i in length(networks_list)){
+ range <- colnames(M[[i]])
+ M_cluster[range,range] <- M[[i]]
+ }
+
+ # Generate final matrix as the sum of the two matrix
+ # distribution of contacts that occur within cluster and outside cluster
+
+ M_final = M_cluster*dist_within_between[1] + M_full*dist_within_between[2]
+ # image(M_final)
+
+ ## merge with a diagonal matrix so that final network has a strong diagonal component
+ ## elements on the diagonal represent the time spent within its own word
+ ## (which realistically should be higher than the time spent in any other ward; we assume at least 50% of the time)
+ for(i in colnames(M_final)){
+ ## time spent within their own patch
+ t_patch <- runif(n = 1, min = 0.5, max = 0.8)
+ M_final[i,] <- M_final[i,]*(1.- t_patch)
+ M_final[i,i] <- M_final[i,i]+t_patch
+ M_final[i,] <- M_final[i,]/sum(M_final[i,])
+ }
+
+ M_final %<>% t
+
+ return(M_final)
+ }
+
+ matContact = generate_network(PARCELLE,
+ group,
+ tot_n,
+ p = c(within_clust_lev,between_clust_lev),
+ dist_within_between = c(clust_ratio_inout,(1-clust_ratio_inout)))
+
+ if(verbose)
+ image.plot(matContact)
+
+ rownames(matContact) <- PARCELLE$ID
+ colnames(matContact) <- PARCELLE$ID
+
+ return(matContact)
+}
diff --git a/R/build_pts_fun.R b/R/build_pts_fun.R
index 333f804113244d0fcd3e2f4b275aa6f70dbf3bfe..0a90087e73efaf22540128689fcdfd95caeac903 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 <- ""
@@ -36,7 +36,12 @@ build_pts_fun <- function(u0, v0, gdata, diapause_interv){
###########################################################
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;
@@ -47,87 +52,94 @@ build_pts_fun <- function(u0, v0, gdata, diapause_interv){
int i;
for (i = 1; i < 4; ++i)
{
- v_new[i] = ldata[1 + (int)ldata[0] * (i-1) + (int)t];
+ v_new[i] = ldata[2 + (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
// eggs layed by infected mosquitoes
int neweggI;
neweggI = u[', which(names(u0) == "Neggs") - 1,'] - v[', which(rownames(v0) == "prevEggs") - 1,'];
- v_new[', which(rownames(v0) == "prevEggs") - 1,'] = u[', which(names(u0) == "Neggs") - 1,']/1000;
+ v_new[', which(rownames(v0) == "prevEggs") - 1,'] = u[', which(names(u0) == "Neggs") - 1,'];
// 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 +147,180 @@ 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"], ') ;
+ double infA1h;
+ // number of h1 gorged in stochastic
+ infA1h = u[', which(names(u0) == "ninfm1") - 1,'] - v[', which(rownames(v0) == "nIm1") - 1,'];
+ v_new[', which(rownames(v0) == "nIm1") - 1,'] = u[', which(names(u0) == "ninfm1") - 1,'];
+ varA1hm = v[', which(rownames(v0) == "Aemm") - 1,'] * ', gdata["gammaAem"], ' - (v[', which(rownames(v0) == "A1hm") - 1,'] - infA1h) * (fmAr + ', gdata["gammaAh"], ') - infA1h ;
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"], ') ;
+ double infA2h;
+ // number of h1 gorged in stochastic
+ infA2h = u[', which(names(u0) == "ninfm2") - 1,'] - v[', which(rownames(v0) == "nIm2") - 1,'];
+ v_new[', which(rownames(v0) == "nIm2") - 1,'] = u[', which(names(u0) == "ninfm2") - 1,'];
+ varA2hm = (v[', which(rownames(v0) == "A1om") - 1,'] + v[', which(rownames(v0) == "A2om") - 1,']) * ', gdata["gammaAo"], ' - (v[', which(rownames(v0) == "A2hm") - 1,'] - infA2h) * (fmAr + ', gdata["gammaAh"], ') - infA2h ;
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 patch j
+ int nMh;
+ // proportion of infected mosquitoes looking for host in patch j
+ 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[2 + (int)ldata[0] * 3 + j];
+ wIm = wIm + pIMh;
+ }
+
+ // prop of infected mosquitoes in patch i
+ nMh = u[',which(colnames(u0) == "A2hmI") - 1,'] + v[', which(rownames(v0) == "A1hm") - 1,'] + v[', which(rownames(v0) == "A2hm") - 1,'];
+ if(nMh > 0){
+ wIm = (double) u[',which(colnames(u0) == "A2hmI") - 1,'] / nMh;
+ }
+
+ 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[2 + (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..a61cc9c9c0098fee25715dbbb86b7026f840d18e 100644
--- a/R/build_transitions.R
+++ b/R/build_transitions.R
@@ -34,16 +34,24 @@ build_transitions <- function(){
stoch_transitions <- c(
#### HUMANS INFECTION ####
- "A2hmI + Sh -> A2hmI * gammaAh * bMH * Sh/(Sh+Eh+Ih+Rh) -> A2gmI + Eh + ninfh", ## ninfh count the number of autochtonous transitions
- "A2hmI -> A2hmI * gammaAh * (1 - bMH * Sh/(Sh+Eh+Ih+Rh)) -> A2gmI",
+ ## ninfh count the number of autochtonous transitions
+
+ # Infection of host in the patch
+ "A2hmI + Sh -> A2hmI * gammaAh * bMH * Sh/(Sh+Eh+Ih+Rh) * pRes -> A2gmI + Eh + ninfh",
+ # Infection of host in another patch
+ "Sh -> gammaAh * bMH * pIm -> Eh + ninfh",
+
+ # Infected mosquito is gorged without infecting anyone
+ "A2hmI -> A2hmI * gammaAh * (1 - bMH * (Sh/(Sh+Eh+Ih+Rh) * pRes + (1 - pRes))) -> A2gmI",
"Eh -> Eh * muH -> Ih",
"Ih -> Ih * rhoH -> Rh",
#### 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",
+ "@ -> (A1hm - (nIm1 - ninfm1)) * gammaAh * pIh * bHM -> A1gmI + ninfm1",
+ "@ -> (A2hm - (nIm2 - ninfm2)) * gammaAh * pIh * bHM -> A2gmI + ninfm2",
"A1gmI -> temperature > TAG ? A1gmI * (temperature - TAG) / TDDAG : 0 -> A1omI",
"A1omI -> A1omI * gammaAo -> A2hmI + 95 * Neggs", ## New eggs
diff --git a/R/check_gdata.R b/R/check_gdata.R
new file mode 100644
index 0000000000000000000000000000000000000000..c0b639fc952cb992eeea9a443c118dc41a900e12
--- /dev/null
+++ b/R/check_gdata.R
@@ -0,0 +1,31 @@
+#' Check and build gdata
+#'
+#' @description function to build gdata
+#'
+#' @usage check_gdata(gdata, vector_species)
+#'
+#' @param gdata list
+#' @param vector_species string
+#'
+#' @return gdata
+
+check_gdata <- function(gdata, vector_species, verbose = T){
+ 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 <- build_gdata(vector_species = "Ae. albopictus",
+ climate = "temperate",
+ verbose = verbose)
+ } 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(build_gdata(verbose = F) %>% names, 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" \n
+ You can use the function build_gdata to generate a well formated dataset')
+
+ gdata
+}
diff --git a/R/create_ldata.R b/R/create_ldata.R
index ecfa3fcbe56ec2ed847b064c9913a16576b0d5c5..2eca89517e3de5f691988231cdfea9f8e1226943 100644
--- a/R/create_ldata.R
+++ b/R/create_ldata.R
@@ -1,63 +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
- ## CHECKS
+ METEO %<>% copy
+ setDT(METEO)
- try(if(NA %in% match(PARCELLE$STATION, meteo$POSTE))
- stop("Some stations are missing from the meteorological dataset"))
+ ## CHECKS
+ 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
## for each administrative unit
+ message("## Generating parameters for all patches can take time, please be patient. ##")
+
+ #make it parallel
- #make it paralell
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
@@ -96,13 +109,25 @@ 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 %<>% .[c(1, seq(nrow(ldata))),]
+
+ldata[2,names(diag(mMov))] <- diag(mMov)
+
+ldata %<>% rbind(., mMov[,colnames(ldata)])
+
+rownames(ldata) <- c("n_days", "pRes",
+ paste0("Tp_", seq(nYR_init * 365 + time_max)),
+ paste0("kL_", seq(nYR_init * 365 + time_max)),
+ paste0("kP_", seq(nYR_init * 365 + time_max)),
+ PARCELLE$ID)
return(ldata)
}
diff --git a/R/diapause.R b/R/diapause.R
index 5158e234a31997b65269ee27a8a9c003581256b9..88f2fba9f3420cf9ab432cb096383e299abfbf09 100644
--- a/R/diapause.R
+++ b/R/diapause.R
@@ -19,10 +19,11 @@
diapause <- function(startFav, endFav, time_max, nYR_init = 2){
# FIX ME add check for date format
+ # FIX ME run it on real date sequence
-x <- startFav %<>% as.POSIXlt(., tryFormats = c("%d/%m/%Y")) %>% yday
+x <- startFav %>% yday
-DurFav <- endFav %>% as.POSIXlt(., tryFormats = c("%d/%m/%Y")) %>% yday %>% subtract(startFav)
+DurFav <- endFav %>% yday %>% subtract(x)
DurUnFav <- 365 - DurFav
while(max(x) < (time_max + nYR_init * 365)) {
diff --git a/R/generate_events.R b/R/generate_events.R
deleted file mode 100644
index 5b97e5ae26a6a38dc6005420f8defef86621c898..0000000000000000000000000000000000000000
--- a/R/generate_events.R
+++ /dev/null
@@ -1,57 +0,0 @@
-#' Generate event objects for introduction of exposed humans
-#'
-#' @description Function used to generate the event matrix and associated objects required by SimInf package
-#'
-#' @usage build_E(timetosample, compartments, nadm)
-#'
-#' @param timetosample Numerical vector. Successive interval of periods to sample introduction
-#' @param nintro Numerical value. Number of random introductions.
-#' @param nodeintro Numerical vector. List of patch to seed random introductions.
-#' @param compartments String vector. Compartments included in the stochastic model.
-#'
-#' @return List with E, N and events matrix required by SimInf. Scheduled introductions are stored in the 'events' matrix.
-#'
-#' @keywords events
-#'
-#' @export
-
-build_E <- function(timetosample, nintro, 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
- 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, #3
- 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, #4
- 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, #5
- 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0 #6
- ),
- .Dim = c(length(compartments),6),
- .Dimnames = list(compartments,
- c("1", "2", "3", "4", "5", "6"))
- )
-
-
-# N matrix for the shift process
-N <- matrix(rep(0, length(compartments)), nrow = length(compartments), ncol = 2,
- dimnames = list(compartments,
- c("1", "2")))
-
-timetosample %<>% split(., ceiling(seq_along(.)/2)) %>%
- lapply(., function(x) paste(x, collapse = ":")) %>%
- unlist %>% paste(., collapse = ",") %>% c("c(", ., ")") %>%
- paste(., collapse ="") %>% parse(text=.) %>% eval(.)
-
-events <- data.frame(
- event = rep("enter", nintro), ## Event "extTrans" is a movement between nodes
- 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
- 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
- )
-
-return(list(E = E,
- N = N,
- events = events))
-}
diff --git a/R/iniState.R b/R/iniState.R
index 41291589aa80b3990953fe4fbd033a6ec8f27d3e..53c74754e5b9436a383681331177a133627d044f 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(
@@ -36,7 +34,9 @@ u0 <- data.frame(
Neggs = rep(0, nadm), #new eggs from infected mosquitoes
# Number of autochtonous human infection
- ninfh = rep(0, nadm)
+ ninfh = rep(0, nadm),
+ ninfm1 = rep(0, nadm),
+ ninfm2 = rep(0, nadm)
)
###############
@@ -51,14 +51,15 @@ u0 <- data.frame(
cont_var <- c(
"z", "temperature", "kL", "kP",
"Em", "Lm", "Pm",
- "Aemm", "A1hm", "A1gm", "A1om", "A2hm", "A2gm", "A2om", "prevEggs",
- "R0")
+ "Aemm", "A1hm", "A1gm", "A1om", "A2hm", "A2gm", "A2om",
+ "prevEggs", "nIm1", "nIm2",
+ "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..d6bdac986954ce49e153fdef960eece174b585f8
--- /dev/null
+++ b/R/runArboRisk.R
@@ -0,0 +1,298 @@
+#' 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 vector_species string
+#' @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)
+#' @param n_sim integer
+#' @param introduction_pts data.frame or data.table
+#'
+#'
+#' @importFrom SimInf mparse
+#' @importFrom SimInf run
+#' @importFrom SimInf trajectory
+#' @importFrom pbapply pblapply
+#'
+#' @return data.table
+#'
+#' @export
+
+
+runArboRisk <- function(PARCELLE,
+ vector_species = "Ae. albopictus",
+ gdata = NULL,
+ ldata = NULL,
+ METEO = NULL,
+ mMov = NULL,
+ start_date = NULL,
+ end_date = NULL,
+ nYR_init = 1,
+ n_sim = 1,
+ introduction_pts = NULL){
+
+ if(is.null(ldata) & is.null(METEO)) stop("you must provide either ldata matrix or meteorological dataset (METEO)")
+
+ if(!is.null(ldata) & !is.null(mMov)) stop("you can not provide a new human contact matrix if you provide ldata. Please generate a new ldata with the current mMov and set mMov to NULL")
+
+ if((is.null(start_date) | is.null(end_date)) & is.null(METEO)) stop("if the meteorological dataset (METEO) is not provided, you must provide start en end dates for simulations (arguments: start_date, end_date)")
+
+ if(n_sim < 1)
+ stop('The number of simulations must be positive a positive integer')
+
+ ####################
+ ### 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"')
+
+ ## 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(!is.null(METEO))
+ 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 = start_date,
+ to = end_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) ####
+
+ gdata <- check_gdata(gdata, vector_species)
+
+ ##################
+ ## 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).
+
+
+ if(is.null(ldata)){
+ message("## Format local data ##")
+ ldata <- build_ldata(PARCELLE,
+ METEO,
+ gdata,
+ vector_species,
+ mMov,
+ start_date,
+ end_date,
+ nYR_init)
+ } else {
+ if(ldata[1,1] != (time_max + 365 * nYR_init))
+ stop(paste0("The number of days considered in ldata (",ldata[1,1],") is not consistent with the simulated period and the number of initialization years (",time_max + 365 * nYR_init," days)."))
+
+ if(ncol(ldata) != nrow(PARCELLE))
+ stop(paste0("The number of patch considered in ldata (",ncol(ldata),") is not consistent with the number of patchs in PARCELLE (",nrow(PARCELLE),")."))
+ }
+
+ ##############
+ ## Diapause ##
+ ##############
+
+ # Calculates days for diapause over all the simulated period and build pts_fun
+ if( !is.null(gdata$startFav) & !is.null(gdata$endFav)){
+ diapause_interv <- diapause(startFav = gdata$startFav,
+ endFav = gdata$endFav,
+ time_max,
+ nYR_init)
+ }
+
+ gdata %<>% .[-which(names(gdata) %in% c("startFav","endFav"))]
+ gdata %<>% unlist
+ mode(gdata) <- "numeric"
+
+
+ message("## Build stochastic transitions ##")
+ #### Stochastic transitions ####
+ stoch_transitions <- build_transitions()
+
+ ###########################################################
+ ### 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)
+
+ ###################
+ ### Simulations ###
+ ###################
+
+ # u0[1, 2] <- 1000
+ # u0[, 7] <- 1000
+
+ # start_time <- Sys.time()
+
+ message("## Build the model ##")
+ # message(Sys.time())
+
+ if(is.null(introduction_pts)){
+
+ 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
+ )
+
+ } else {
+
+ ###########################################
+ ### Introduction of exposed individuals ###
+ ###########################################
+
+ # E Matrix for the select process
+ E <- structure(.Data = c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, #1
+ 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, #2
+ 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, #3
+ 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, #4
+ 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, #5
+ 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0 #6
+ ),
+ .Dim = c(length(u0 %>% colnames),6),
+ .Dimnames = list(u0 %>% colnames,
+ c("1", "2", "3", "4", "5", "6"))
+ )
+
+ # N matrix for the shift process
+ N <- matrix(rep(0, length(u0 %>% colnames)), nrow = length(u0 %>% colnames), ncol = 2,
+ dimnames = list(u0 %>% colnames,
+ c("1", "2")))
+
+ # turn dates of introduction_pts into time points
+ introduction_pts$time %<>% match(., seq.Date(from = (start_date - (365*nYR_init)),
+ to = end_date,
+ by ="day"))
+
+ # turn IDs into numerical IDs
+ introduction_pts$node %<>% match(., PARCELLE$ID)
+
+ # select the stage
+ introduction_pts$select %<>% match(., u0 %>% colnames)
+
+ introduction_pts$event = "enter" ## Event "extTrans" is a movement between nodes
+ introduction_pts$dest = 0 ## Which node is the destination node
+ introduction_pts$proportion = 0 ## This is not used when n > 0
+ introduction_pts$shift = 0 ## Not used in this example
+
+ if(introduction_pts %>% is.na %>% sum %>% is_greater_than(0))
+ stop("introduction_pts format is not respected. \n
+ ids in 'node' columns must be in PARCELLE ID. \n
+ Dates in 'time' must be included in the simulated period. \n
+ Epidemiological stages in 'select' must be among the following stages: 'Sh', 'Eh', 'Ih', 'Rh', 'A1gmI', 'A1omI', 'A2hmI', 'A2gmI', 'A2omI'.")
+
+ 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 = introduction_pts,
+ E = E,
+ N = N
+ )
+ }
+
+
+ message("## Run simulations ##")
+ output <- pblapply(1:n_sim, function(x){
+ traj <- model %>% run %>% trajectory
+ setDT(traj)
+ traj[, DATE := seq.Date(from = (start_date - (365*nYR_init)),
+ to = end_date,
+ by ="day"), by = node]
+
+ traj %<>% .[DATE >= start_date]
+ traj[,time := time - min(time)]
+
+ traj[, ID := PARCELLE$ID[node]]
+
+ setcolorder(traj, c("ID", "node",
+ "DATE", "time",
+ "Sh", "Eh", "Ih", "Rh",
+ "Em", "Lm", "Pm", "Aemm", "A1hm", "A1gm", "A1om", "A2hm", "A2gm", "A2om", "A1gmI", "A1omI", "A2hmI", "A2gmI", "A2omI",
+ "R0", "ninfh", "ninfm1", "ninfm2",
+ "z", "temperature", "kL", "kP"))
+
+ traj[, `:=`(Neggs = NULL,
+ prevEggs = NULL,
+ nIm1 = NULL,
+ nIm2 = NULL,
+ pIh = NULL,
+ pIm = NULL
+ )]
+
+ if(is.null(introduction_pts)){
+ traj[, `:=`(ninfh = NULL,
+ ninfm1 = NULL,
+ ninfm2 = NULL
+ )]
+ }
+
+ })
+
+ return(output)
+}
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 a8fa94bdddd52aa0e592add03843cc3a4c2bffd3..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/build_E.Rd b/man/build_E.Rd
deleted file mode 100644
index a393ce2181d7ce051fb3aed4d71a8ad1ead87701..0000000000000000000000000000000000000000
--- a/man/build_E.Rd
+++ /dev/null
@@ -1,24 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/generate_events.R
-\name{build_E}
-\alias{build_E}
-\title{Generate event objects for introduction of exposed humans}
-\usage{
-build_E(timetosample, compartments, nadm)
-}
-\arguments{
-\item{timetosample}{Numerical vector. Successive interval of periods to sample introduction}
-
-\item{nintro}{Numerical value. Number of random introductions.}
-
-\item{compartments}{String vector. Compartments included in the stochastic model.}
-
-\item{nodeintro}{Numerical vector. List of patch to seed random introductions.}
-}
-\value{
-List with E, N and events matrix required by SimInf. Scheduled introductions are stored in the 'events' matrix.
-}
-\description{
-Function used to generate the event matrix and associated objects required by SimInf package
-}
-\keyword{events}
diff --git a/man/build_E_random.Rd b/man/build_E_random.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..cfe002588f6e826abd488ee12318e19dfffc9928
--- /dev/null
+++ b/man/build_E_random.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/build_E_random.R
+\name{build_E_random}
+\alias{build_E_random}
+\title{Generate the list of introduction event}
+\usage{
+build_E_random(period_start, period_end, n_ind = 1, n_events = 1, stage = "Eh", loc)
+}
+\arguments{
+\item{period_start}{date in '\%Y-\%m-\%d' format. Define the beginning of the period for introduction.}
+
+\item{period_end}{date in '\%Y-\%m-\%d' format. Define the end of the period for introduction.}
+
+\item{n_ind}{Numerical vector. Range of number of individuals introduced.}
+
+\item{n_events}{String vector. Number of introduction generated.}
+
+\item{stage}{String vector. Epidemiological stage of introduced individuals.}
+
+\item{loc}{String or numerical vector. List of patch to seed random introductions.}
+}
+\value{
+events matrix with scheduled introductions.
+}
+\description{
+Function used to generate the event matrix and associated objects required by SimInf package
+}
+\examples{
+build_E_random(period_start = as.Date("2020/03/10"), period_end = as.Date("2020/09/30"), n_ind = 1, n_events = 10, stage = "Eh", loc = LETTERS[1:3])
+
+}
+\keyword{events}
diff --git a/man/build_gdata.Rd b/man/build_gdata.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..4f7e15faa979d231280527b15a0231565fe8ab7b
--- /dev/null
+++ b/man/build_gdata.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/build_gdata.R
+\name{build_gdata}
+\alias{build_gdata}
+\title{Generate global parameters dataset}
+\usage{
+build_gdata()
+}
+\arguments{
+\item{vector_species}{}
+}
+\value{
+Named list of parameters
+}
+\description{
+Then the functions formats all global parameters required for the transitions as a named list readable by SimInf package.
+}
+\examples{
+build_gdata()
+
+}
diff --git a/man/build_ldata.Rd b/man/build_ldata.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..b3add023d68c46e455d2335fa2b9fd05cd5c1c82
--- /dev/null
+++ b/man/build_ldata.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/build_ldata.R
+\name{build_ldata}
+\alias{build_ldata}
+\title{Build ldata}
+\usage{
+build_ldata(PARCELLE, METEO)
+}
+\arguments{
+\item{PARCELLE}{data.frame or data.table}
+
+\item{METEO}{data.frame or data.table}
+
+\item{gdata}{list}
+
+\item{vector_species}{string}
+
+\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{
+matrix
+}
+\description{
+function to build ldata
+}
diff --git a/man/build_mMov.Rd b/man/build_mMov.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..122878197b11d29074af73323b03ba904bbd6af4
--- /dev/null
+++ b/man/build_mMov.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/build_mMov.R
+\name{build_mMov}
+\alias{build_mMov}
+\title{Creates a synthetic network}
+\usage{
+build_mMov()
+}
+\arguments{
+\item{PARCELLE}{data.table}
+
+\item{group}{string}
+
+\item{within_clust_lev}{numeric. Probability of creating links between wards of the same building}
+
+\item{between_clust_lev}{numeric. Probability of creating links between wards of different building}
+
+\item{clust_ratio_inout}{numeric. Ratio of intra-building clustering (relative to inter-building clustering)}
+
+\item{verbose}{logical. Print plot if TRUE}
+}
+\value{
+a matrix with random probabilities of contact
+}
+\description{
+This function creates a synthetic human mobility network.
+You can use this script to generate different input data to be used as examples of "network structures".
+}
diff --git a/man/check_gdata.Rd b/man/check_gdata.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..f33918c644e8375feee4281df4fc869097be910b
--- /dev/null
+++ b/man/check_gdata.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/check_gdata.R
+\name{check_gdata}
+\alias{check_gdata}
+\title{Check and build gdata}
+\usage{
+check_gdata(gdata, vector_species)
+}
+\arguments{
+\item{gdata}{list}
+
+\item{vector_species}{string}
+}
+\value{
+gdata
+}
+\description{
+function to build gdata
+}
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/hello.Rd b/man/hello.Rd
deleted file mode 100644
index 0fa7c4b8817591c2dff2b3997d2566320ac6d9fc..0000000000000000000000000000000000000000
--- a/man/hello.Rd
+++ /dev/null
@@ -1,12 +0,0 @@
-\name{hello}
-\alias{hello}
-\title{Hello, World!}
-\usage{
-hello()
-}
-\description{
-Prints 'Hello, world!'.
-}
-\examples{
-hello()
-}
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..2b2973639dd76fd52103ba9a5ceac1f20685fd36
--- /dev/null
+++ b/man/runArboRisk.Rd
@@ -0,0 +1,35 @@
+% 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{vector_species}{string}
+
+\item{gdata}{list}
+
+\item{METEO}{data.frame or data.table}
+
+\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)}
+
+\item{n_sim}{integer}
+
+\item{introduction_pts}{data.frame or data.table}
+}
+\value{
+data.table
+}
+\description{
+function to run simulations
+}
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000000000000000000000000000000000000..79b6b99567223433766f8cb7c02b31071f879afd
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,12 @@
+# This file is part of the standard setup for testthat.
+# It is recommended that you do not modify it.
+#
+# Where should you do additional test configuration?
+# Learn more about the roles of various files in:
+# * https://r-pkgs.org/tests.html
+# * https://testthat.r-lib.org/reference/test_package.html#special-files
+
+library(testthat)
+library(ArboRisk)
+
+test_check("ArboRisk")
diff --git a/tests/testthat/test-runArboRisk.R b/tests/testthat/test-runArboRisk.R
new file mode 100644
index 0000000000000000000000000000000000000000..051dec441344a460bb35bfb8a5cd8884753ed827
--- /dev/null
+++ b/tests/testthat/test-runArboRisk.R
@@ -0,0 +1,12 @@
+test_that("simulation works", {
+
+
+ traj <- runArboRisk(PARCELLE = PARCELLE,
+ METEO = METEO,
+ n_sim = 2)
+
+ expect_type(traj, "list")
+ expect_s3_class(traj[[1]], "data.frame")
+ expect_length(traj, 2)
+
+})
diff --git a/tests/testthat/test-test_intro.R b/tests/testthat/test-test_intro.R
new file mode 100644
index 0000000000000000000000000000000000000000..8da7c70283b78d2b7a0bc62db26fe732ca5d4d48
--- /dev/null
+++ b/tests/testthat/test-test_intro.R
@@ -0,0 +1,11 @@
+test_that("introduction works", {
+ introduction_pts <- build_E_random(
+ period_start <- "2020/03/10" %>% as.Date,
+ period_end <- "2020/09/30" %>% as.Date,
+ n_ind = 1:10,
+ n_events = 5,
+ stage = "Eh",
+ loc = letters[1:5])
+ expect_equal(dim(introduction_pts), c(5,4))
+ expect_equal(names(introduction_pts), c("time", "node", "n", "select"))
+})
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 b600b3e2b06a92703cd9159a93f23102d72afe5c..3e6ff3abe008a70568d5f1d0ee901b758e3a9088 100644
--- a/work_in_progress/main.R
+++ b/work_in_progress/main.R
@@ -1,4 +1,4 @@
-### Date of change: 30/01/2023
+### Date of change: 20/02/2023
### Authors: Pachka, Ewy, Elena
## This script is under development at the UMR Astre - Cirad.
@@ -14,161 +14,95 @@ gc() # free up memory and report the memory usage.
###################################
### Load packages and functions ###
###################################
+devtools::document()
devtools::install(quick = T)
+devtools::test()
+
+library(ArboRisk)
library(data.table)
library(magrittr)
-library(ggplot2)
-
-library(ArboRisk)
-library(SimInf)
-# Parallel processing
-library(parallel)
-library(doParallel)
-library(foreach)
+library(hydroTSM)
####################
### General data ###
####################
# input data:
# PARCELLE : table where each row is a node and columns are attributes
-# meteo: table whit four columns nodeID, date, daily rainfall, daily mean temperature
+# 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 <- 400 # number of days simulated
-nYR_init <- 2 # number of years used to initialize the population dynamics
+PARCELLE %<>% .[startsWith(PARCELLE$ID, "34"),]
+# PARCELLE %<>% .[1:100,]
+METEO
###########################################
-## Create u0, v0 and compartment objects ##
+### Introduction of exposed individuals ###
###########################################
-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_pts <- build_E_random(
+ period_start <- "2020/03/10" %>% as.Date,
+ period_end <- "2020/09/30" %>% as.Date,
+ n_ind = 1:10,
+ n_events = 5,
+ stage = "Eh",
+ loc = PARCELLE$ID)
-###########################################
-### Introduction of exposed individuals ###
-###########################################
-mov <- build_E(timetosample = diapause_interv,
- nintro = 10,
- u0 %>% colnames,
- nodeintro = nrow(PARCELLE))
+######################################
+### Introduction of human mobility ###
+######################################
-E <- mov$E
-N <- mov$N
-events <- mov$events
+mMov <- build_mMov(PARCELLE,
+ group = "NOM_COM",
+ within_clust_lev = 1,
+ between_clust_lev = 0.1,
+ clust_ratio_inout = 0.8,
+ verbose = T)
+
+ldata <- build_ldata(PARCELLE,
+ METEO,
+ vector_species = "Ae. albopictus")
###################
### Simulations ###
###################
+# without introductions
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 <- runArboRisk(PARCELLE = PARCELLE,
+ vector_species = "Ae. albopictus",
+ METEO = METEO,
+ n_sim = 30)
-traj <- trajectory(result)
+end_time <- Sys.time()
+end_time - start_time
+
+# without introductions
+start_time <- Sys.time()
+
+traj <- runArboRisk(PARCELLE = PARCELLE,
+ vector_species = "Ae. albopictus",
+ ldata = ldata,
+ n_sim = 30,
+ start_date = min(METEO$DATE) %>% as.Date,
+ end_date = max(METEO$DATE) %>% as.Date)
+
+end_time <- Sys.time()
+end_time - start_time
+
+# with introductions
+start_time <- Sys.time()
+
+traj <- runArboRisk(PARCELLE = PARCELLE,
+ vector_species = "Ae. albopictus",
+ ldata = ldata,
+ n_sim = 30,
+ start_date = min(METEO$DATE) %>% as.Date,
+ end_date = max(METEO$DATE) %>% as.Date,
+ introduction_pts = introduction_pts)
end_time <- Sys.time()
end_time - start_time
@@ -177,9 +111,14 @@ end_time - start_time
##### Explore output
#####
-traj %>% setDT
traj
+traj[, date2season := time2season(time2date, out.fmt = "seasons")]
+t2d <- rev(ymd('2022/11/30') - days(0:max(traj$time)))
+traj[, time2date := t2d[(time+1)]]
+traj[, date2season := time2season(time2date, out.fmt = "seasons")]
+traj %<>% .[time > (max(time)-time_max)]
+
# load polygons of PARCELLES
filename <- "test/data/IRIS-GE_2-0_SHP_LAMB93_FXX-2022/IRIS_GE.SHP"
PARCELLES <- vect(filename)
diff --git a/work_in_progress/processInput.R b/work_in_progress/processInput.R
index f53f17d8f09bafce7fc095053a7725a89e5db2da..be5f63f5a9de12cd90cbde2b55e0743c6c18cf79 100644
--- a/work_in_progress/processInput.R
+++ b/work_in_progress/processInput.R
@@ -1,25 +1,22 @@
library(terra)
-# library(tidyterra)
library(geodata)
-# library(raster)
-# library(sp)
-# library(sf)
#### 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 ###
-STATIONS <- read.csv2("test/data/stations.txt", sep = ";", dec = ".")
+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,
@@ -28,7 +25,7 @@ STATIONS[, `:=`(TYPE = NULL,
STATIONS[, ALTITUDE := substr(ALTITUDE, 1, 4) %>% as.numeric]
#### List of administrative unit ####
-filename <- "test/data/IRIS-GE_2-0_SHP_LAMB93_FXX-2022/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
@@ -50,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
@@ -61,10 +58,11 @@ PARCELLE$ALTITUDE_UNIT %>% is.na %>% sum %>% equals(0)
# Ignore the error
# https://github.com/sneumann/xcms/issues/288
PARCELLE[is.na(PARCELLE$ALTITUDE_UNIT)]$ALTITUDE_UNIT <- 0
+# FIX ME: thinnk about it
## 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,9 +70,9 @@ PARCELLE$POP <- terra::extract(pop, PARCELLE, fun=sum, na.rm=TRUE) %>% .$populat
# Absence of missing values?
PARCELLE$POP %>% is.na %>% sum %>% equals(0)
# Turn density into number of persons
-PARCELLE$POP %<>% multiply_by(PARCELLE$SURF_HA/100)
+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
@@ -87,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)