diff --git a/.gitignore b/.gitignore
index c833a2c666716a8e6dc01a2456b9f64ca6a2392b..20d892b5c4464af84383e6027af617dfd9785878 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,4 @@
.RData
.Ruserdata
inst/doc
+.RDataTmp
diff --git a/NAMESPACE b/NAMESPACE
index 0e970517d4c247f54b96dcb0623e7741dc921f33..7545b745280cabc17f597fbd10f919db5f9b14c0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,11 +4,8 @@ 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)
@@ -16,11 +13,21 @@ import(igraph)
importFrom(SimInf,mparse)
importFrom(SimInf,run)
importFrom(SimInf,trajectory)
+importFrom(TDLM,check_format_names)
+importFrom(TDLM,extract_opportunities)
+importFrom(TDLM,extract_spatial_information)
+importFrom(TDLM,run_law)
+importFrom(data.table,as.IDate)
+importFrom(dplyr,anti_join)
importFrom(fields,image.plot)
importFrom(igraph,as_adjacency_matrix)
importFrom(igraph,erdos.renyi.game)
importFrom(imputeTS,na_ma)
+importFrom(magrittr,"%<>%")
importFrom(magrittr,"%>%")
+importFrom(magrittr,is_greater_than)
+importFrom(magrittr,is_less_than)
+importFrom(magrittr,subtract)
importFrom(parallel,clusterEvalQ)
importFrom(parallel,clusterExport)
importFrom(parallel,detectCores)
@@ -28,3 +35,4 @@ importFrom(parallel,makeCluster)
importFrom(parallel,stopCluster)
importFrom(pbapply,pbapply)
importFrom(pbapply,pblapply)
+importFrom(zoo,rollapply)
diff --git a/R/build_E_random.R b/R/build_E_random.R
index 5d34d8f2df7dc0ee7526a377c2d9bf7652a91a39..d66a60e8f8259bb8f292d3d7ea6f8f2352a6cd6d 100644
--- a/R/build_E_random.R
+++ b/R/build_E_random.R
@@ -1,6 +1,8 @@
#' Generate the list of introduction event
#'
-#' @description Function used to generate the event matrix and associated objects required by SimInf package
+#' @description Function used to generate the random introduction.
+#' The introduction events generated can occurs in different locations and at different dates in a given period.
+#' The number of individuals introduced and their epidemiological stages can be deterministic (unique value provided) or random (randomly selected in a vector of values).
#'
#' @usage build_E_random(period_start, period_end, n_ind = 1, n_events = 1, stage = "Eh", loc)
#'
@@ -12,7 +14,7 @@
#' @param loc String or numerical vector. List of patch to seed random introductions.
#'
#'
-#' @return events matrix with scheduled introductions.
+#' @return events matrix with scheduled introductions
#'
#' @keywords events
#'
@@ -28,18 +30,42 @@ build_E_random <- function( period_start,
stage = "Eh",
loc){
+ ###
+ ### Checks
+ ###
+ if(!inherits(period_start, "Date") | !inherits(period_end, "Date"))
+ stop("start and end dates must be dates ('Date' class)")
+
+
+ if(!inherits(n_ind, c("integer","numeric")))
+ stop("n_ind must be numeric")
+
+ if(!inherits(n_events, "numeric") | length(n_events) != 1 )
+ stop("n_events must be a unique numeric value")
+
+ list_compartments() %>% .$u0
+
+ if(! inherits(stage, "character") | sum(is.na(match(stage, list_compartments() %>% .$u0))) > 0)
+ stop(paste(c("stage must be a character vector. Authorized entries:", list_compartments() %>% .$u0), collapse = " "))
+
+ ###
+
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)
+ ###
+ ###
+
+ # sample time, location and number of individuals
+ sampling <- expand.grid(period, loc, n_ind, stage)
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
+ dest = 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
+ select = sampling$Var4 ## Use the 4th column in the model select matrix
)
return(events)
diff --git a/R/build_gdata.R b/R/build_gdata.R
index 655873c535b14c9e6575cb52f540303ed1987cbb..3e33815a7b12bd942ec33d8a032d75365786b9de 100644
--- a/R/build_gdata.R
+++ b/R/build_gdata.R
@@ -1,13 +1,65 @@
-#' Generate global parameters dataset
+#' Generate global parameters list
#'
#' @description
-#' Then the functions formats all global parameters required for the transitions as a named list readable by SimInf package.
+#' The functions formats all global parameters required for the model.
#'
#' @usage build_gdata()
#'
#' @importFrom magrittr %>%
#'
-#' @param vector_species
+#' @param vector_species string. Default is "Ae. albopictus" in "temperate" climate. Unique dataset implemented to date.
+#' @param climat string. Default is "Ae. albopictus" in "temperate" climate. Unique dataset implemented to date.
+#'
+#' @param bMH numeric. Probability of infection from vector to host when a human is bitten by an infected mosquito.
+#' @param bHM numeric. Probability of infection from host to vector when an infected human is bitten by an susceptible mosquito.
+#' @param muH numeric. Daily transition rate from exposed (E) to infected (I) for humans (1/days)
+#' @param rhoH numeric. Daily recovery rate of humans (1/days)
+#'
+#' @param muE numeric. Daily mortality rate of eggs (1/days)
+#' @param TE numeric. Minimal temperature needed for egg development (°C)
+#' @param TDDE numeric. Total number of degree-day necessary for egg development (°C)
+#'
+#' @param mu1L numeric. Parameter for the larvae mortality function: \eqn{ mu1L \times e^{(temperature_{t}-10) \times mu2L} + mu3L}
+#' @param mu2L numeric. Parameter for the larvae mortality function: \eqn{ mu1L \times e^{(temperature_{t}-10) \times mu2L} + mu3L}
+#' @param mu3L numeric. Parameter for the larvae mortality function: \eqn{ mu1L \times e^{(temperature_{t}-10) \times mu2L} + mu3L}
+#'
+#' @param q1L numeric. Parameter for the function of transition from larva to pupa \eqn{q1L \times temperature_{t}^{2} + q2L \times temperature_{t} + q3L}
+#' @param q2L numeric. Parameter for the function of transition from larva to pupa \eqn{q1L \times temperature_{t}^{2} + q2L \times temperature_{t} + q3L}
+#' @param q3L numeric. Parameter for the function of transition from larva to pupa \eqn{q1L \times temperature_{t}^{2} + q2L \times temperature_{t} + q3L}
+#'
+#' @param mu1P numeric. Parameter for the larvae mortality function: \eqn{mu1P \times e^{(temperature_{t}-10) \times mu2P} + mu3P}
+#' @param mu2P numeric. Parameter for the larvae mortality function: \eqn{mu1P \times e^{(temperature_{t}-10) \times mu2P} + mu3P}
+#' @param mu3P numeric. Parameter for the larvae mortality function: \eqn{mu1P \times e^{(temperature_{t}-10) \times mu2P} + mu3P}
+#'
+#' @param muEM numeric. Daily mortality rate of emerging adults during the emergence (1/days)
+#'
+#' @param q1P numeric. Parameter for the function of transition from pupae to emerging adult \eqn{q1P \times temperature_{t}^{2} + q2P \times temperature_{t} + q3P}
+#' @param q2P numeric. Parameter for the function of transition from pupae to emerging adult \eqn{q1P \times temperature_{t}^{2} + q2P \times temperature_{t} + q3P}
+#' @param q3P numeric. Parameter for the function of transition from pupae to emerging adult \eqn{q1P \times temperature_{t}^{2} + q2P \times temperature_{t} + q3P}
+#'
+#' @param mu1A numeric. Parameter for the adult mortality function: \eqn{mu1A \times e^{(temperature_{t}-10) \times mu2A} + mu3A}
+#' @param mu2A numeric. Parameter for the adult mortality function: \eqn{mu1A \times e^{(temperature_{t}-10) \times mu2A} + mu3A}
+#' @param mu3A numeric. Parameter for the adult mortality function: \eqn{mu1A \times e^{(temperature_{t}-10) \times mu2A} + mu3A}
+#'
+#'
+#' @param gammaAem numeric. Daily development rate of emerging adults (1/days)
+#' @param sigma numeric. Sex-ratio at the emergence
+#'
+#' @param gammaAh numeric. Daily transition rate from host-seeking to engorged adults (1/days)
+#' @param muR numeric. Daily mortality rate related to seeking behavior (1/days)
+#'
+#' @param TAG numeric. Minimal temperature needed for egg maturation (°C)
+#' @param TDDAG numeric. Total number of degree-days necessary for egg maturation (°C)
+#'
+#' @param gammaAo numeric. Daily transition rate from oviposition site-seeking to host-seeking adults (1/days)
+#'
+#' @param beta1 numeric. Number of eggs laid by ovipositing nulliparous females (per female)
+#' @param beta2 numeric. Number of eggs laid by ovipositing parous females (per female)
+#'
+#' @param startFav date. First day of the favorable period for the mosquito (depend on the species and the environment).
+#' @param endFav date. Last day of the favorable period for the mosquito (depend on the species and the environment).
+#'
+#' @param verbose logical. Provide information on parameters during generation
#'
#' @return Named list of parameters
#'
@@ -19,10 +71,11 @@ build_gdata <- function(
vector_species = "Ae. albopictus",
climate = "temperate",
##### Infection
- # HUMAN
- # Infection probability from vector to host
+ # mosquitoes infection from host to vector
+ bHM = NULL,
+ # 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
@@ -30,16 +83,23 @@ build_gdata <- function(
TE = NULL,
TDDE = NULL,
# larva
- # muL = 0.08,
+ mu1L = NULL,
+ mu2L = NULL,
+ mu3L = NULL,
q1L = NULL,
q2L = NULL,
q3L = NULL,
muEM = NULL,
+ mu1P = NULL,
+ mu2P = NULL,
+ mu3P = NULL,
q1P = NULL,
q2P = NULL,
q3P = NULL,
# emerging adults
- muA = NULL,
+ mu1A = NULL,
+ mu2A = NULL,
+ mu3A = NULL,
gammaAem = NULL,
sigma = NULL,
# host-seeking adults
@@ -52,8 +112,6 @@ build_gdata <- function(
gammaAo = NULL,
beta1 = NULL,
beta2 = NULL,
- # mosquitoes infection
- bHM = NULL,
startFav = NULL,
endFav = NULL,
verbose = T) {
@@ -77,14 +135,22 @@ build_gdata <- function(
if(is.null(muE)) muE = 0.05
if(is.null(TE)) TE = 10
if(is.null(TDDE)) TDDE = 110
+ if(is.null(mu1L)) mu1L = 0.0007
+ if(is.null(mu2L)) mu2L = 0.1838
+ if(is.null(mu3L)) mu3L = 0.02
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(mu1P)) mu1P = 0.0003
+ if(is.null(mu2P)) mu2P = 0.2228
+ if(is.null(mu3P)) mu3P = 0.02
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(mu1A)) mu1A = 0.0003
+ if(is.null(mu2A)) mu2A = 0.1745
+ if(is.null(mu3A)) mu3A = 0.025
if(is.null(gammaAem)) gammaAem = 0.4
if(is.null(sigma)) sigma = 0.5
if(is.null(gammaAh))gammaAh = 0.3
@@ -106,7 +172,6 @@ build_gdata <- function(
muH,
muE,
muEM,
- muA,
gammaAem,
sigma,
gammaAh,
@@ -147,16 +212,24 @@ build_gdata <- function(
TE = TE, # 10.4 in Tran 2013
TDDE = TDDE,
# larva
+ mu1L = mu1L,
+ mu2L = mu2L,
+ mu3L = mu3L,
q1L = q1L,
q2L = q2L,
q3L = q3L,
# pupa
muEM = muEM,
+ mu1P = mu1P,
+ mu2P = mu2P,
+ mu3P = mu3P,
q1P = q1P,
q2P = q2P,
q3P = q3P,
# emerging adults
- muA = muA, # 0.02 dans le papier // 0.025 dans ocelet
+ mu1A = mu1A,
+ mu2A = mu2A,
+ mu3A = mu3A,
gammaAem = gammaAem,
sigma = sigma, # proportion of females
# host-seeking adults
diff --git a/R/build_ldata.R b/R/build_ldata.R
index f3b4228b5b232e2f8fd6707c561a643022911dd1..e5085edf8b4eb894536769c2cb18b19b52cf0014 100644
--- a/R/build_ldata.R
+++ b/R/build_ldata.R
@@ -1,33 +1,41 @@
-#' Build ldata
+#' Generate local parameters matrix
#'
-#' @description function to build ldata
+#' @description
+#' The function computes and formats for each pacth the list of local parameters required by the model. Arguments related to the simulation are required to compile daily meteorological data for the simulated period.
#'
#' @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)
+#' @param PARCELLE data.frame or data.table describing the patches. Required columns: "ID": unique identifier, "POP": population size, "SURF_HA": surface in hectare, "KLfix": fix carrying capacity for larvae, "KLvar": variable carrying capacity for larvae, "KPfix": fix carrying capacity for pupae, "KPvar": variable carrying capacity for pupae, "STATION": meteorological station identifier, "DIFF_ALT": difference between meteorological station altitude and average altitude of the patch.
+#' @param METEO data.frame or data.table reporting the daily meteorological data for each meteorological station. Required columns: 'ÌD', 'DATE', 'RR': daily precipitation (mm), 'TP': daily mean temperature (°C)
+#' @param gdata list of parameters. Can be generated using `build_gdata` function
+#' @param vector_species string. Default is "Ae. albopictus" in "temperate" climate. Unique dataset implemented to date. Only needed if gdata is not provided.
+#' @param climat string. Default is "Ae. albopictus" in "temperate" climate. Unique dataset implemented to date. Only needed if gdata is not provided.
+#' @param mMov list of two matrices. Probabilities of location of individuals.
+#' proba_ij is the probabilities of individuals from i (columns) to be bitten in j (rows).
+#' proba_ji is the probabilities for mosquito in i (rows) to bite an individual from j (columns).
+#' Row sums should be equal to 1.
+#' @param start_date date in '\%Y-\%m-\%d' format. Starting date of the simulation. If not provided, the function uses the oldest date in the meteorological dataset.
+#' @param end_date date in '\%Y-\%m-\%d' format. Ending date of the simulation. If not provided, the function uses the most recent date in the meteorological dataset.
+#' @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
+#' @importFrom dplyr anti_join
+#' @importFrom magrittr %<>%
+#' @importFrom magrittr is_greater_than
#'
#' @return matrix
#'
#' @export
-
build_ldata <- function(PARCELLE,
METEO,
gdata = NULL,
vector_species = "Ae. albopictus",
+ climate = "temperate",
mMov = NULL,
start_date = NULL,
end_date = NULL,
@@ -41,19 +49,7 @@ build_ldata <- function(PARCELLE,
# 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"')
+ check_PARCELLE(PARCELLE)
## Inform user on number of PARCELLES, include possibility to use shape
## check PARCELLE columns
@@ -68,14 +64,6 @@ build_ldata <- function(PARCELLE,
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 ###
############################
@@ -92,8 +80,18 @@ build_ldata <- function(PARCELLE,
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
+ TS_sim <- sim_time_serie(start_date, end_date, nYR_init)
+
+ ### Check if all meteorological data are available
+ test_meteo <- expand.grid(ID = METEO$ID %>% unique, DATE = seq(from = start_date,
+ to = end_date, by = "days"))
+ test_meteo$ID %<>% as.numeric
+ test_meteo$DATE %<>% as.Date
+
+ if(anti_join(test_meteo, METEO, by = c("ID", "DATE")) %>% nrow %>% is_greater_than(0))
+ stop("daily meteorological records are missing for the period")
+
+ rm(test_meteo)
###############################################################
## Define the compartments, gdata and stochastic transitions ##
@@ -102,7 +100,7 @@ build_ldata <- function(PARCELLE,
#### Global data (general parameters) ####
- gdata <- check_gdata(gdata, vector_species, verbose =F)
+ gdata <- check_gdata(gdata, vector_species, climate, verbose = F)
##############
## Diapause ##
@@ -124,10 +122,13 @@ build_ldata <- function(PARCELLE,
mMov <- diag(nrow(PARCELLE))
rownames(mMov) <- PARCELLE[, ID]
colnames(mMov) <- PARCELLE[, ID]
+
+ mMov <- list(proba_ij = mMov,
+ proba_ji = mMov)
message("The current scenario do not consider any human mobility (mMov = NULL)")
}
- ldata <- create_ldata(PARCELLE, METEO, gdata, mMov, time_max, nYR_init)
+ ldata <- create_ldata(PARCELLE, METEO, gdata, mMov, TS_sim)
ldata
}
diff --git a/R/build_mMov.R b/R/build_mMov.R
index 25138511a8997bc9e64ff6b87fb5ea18b1233816..95c02b06f742196554678200917302a146bd4755 100644
--- a/R/build_mMov.R
+++ b/R/build_mMov.R
@@ -18,10 +18,8 @@
#' @usage build_mMov()
#'
#' @return a matrix with random probabilities of contact
-#'
-#' @export
-build_mMov <- function(PARCELLE,
+build_mMov_erdosrenyigame <- function(PARCELLE,
group,
within_clust_lev = 0.8,
between_clust_lev = 0.1,
diff --git a/R/build_mMov_TDLM_based.R b/R/build_mMov_TDLM_based.R
new file mode 100644
index 0000000000000000000000000000000000000000..cea9a8e480feece6acc0e5801565c48afb46767d
--- /dev/null
+++ b/R/build_mMov_TDLM_based.R
@@ -0,0 +1,107 @@
+#' Build matrix of contact probabilities
+#'
+#' @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 SpatVec SpatVector or sf. Spatial vector of polygons representing patches. Polygons must include 'POP' and 'ID' attributes.
+#' @param law string. Law used to calculate probabilities. Default is "Unif". See documentation of TDLM::run_law function for law options details.
+#' @param param numeric. Parameter used to adjust the importance of distance or opportunity associated with the chosen law.
+#' A single value or a vector of several parameter values can be used.
+#' Not necessary for the original radiation law or the uniform law. (see TDLM package for more details)
+#' @param p2move numeric. Probability to be in the residential patch.
+#' @param pCom numeric. Proportion of cummuters per patch
+#'
+#' @importFrom TDLM extract_spatial_information
+#' @importFrom TDLM check_format_names
+#' @importFrom TDLM run_law
+#' @importFrom TDLM extract_opportunities
+#'
+#' @usage build_mMov(SpatVec)
+#'
+#' @return list of normalized matrices for probabilities of contact of agent from patch i with agent of patch j in j (proba_ij) and probabilities of contact of agent from patch j with agent of patch i in i (proba_ji)
+#'
+#' @export
+
+build_mMov <- function(SpatVec,
+ law = "Unif",
+ param = NULL,
+ p2move = 0.7,
+ pCom = 0.7){
+
+ message("Checks")
+ if(!inherits(SpatVec, c("SpatVector", "sf")))
+ stop("SpatVec must be eaither a SpatVector or a sf object")
+
+ if(!("ID" %in% names(SpatVec)) | !("POP" %in% names(SpatVec)))
+ stop("SpatVec must include 'POP' and 'ID' as polygons attributes.")
+
+ ## turn to sf
+ if(inherits(SpatVec, c("SpatVector"))){
+ message("Turn to sf")
+ SpatVec %<>% sf::st_as_sf(.)
+ }
+
+ message("Build mi, mj, Oi")
+ mass <- data.frame(Population = round(SpatVec$POP),
+ Outcommuters = round(pCom * SpatVec$POP))
+
+ row.names(mass) <- SpatVec$ID
+
+ mi <- as.numeric(mass[,1])
+ names(mi) <- rownames(mass)
+
+ mj <- mi
+
+ Oi <- as.numeric(mass[,2])
+ names(Oi) <- rownames(mass)
+
+ message("Compute distance between patches")
+ spi <- extract_spatial_information(SpatVec, id = "ID")
+ distance <- spi$distance
+
+ check_format_names(vectors = list(mi = mi, mj = mj, Oi = Oi),
+ matrices = list(distance = distance),
+ check = "format_and_names")
+
+ message("Compute opportunity vector")
+ opportunity <- mass[, 1]
+ names(opportunity) <- rownames(mass)
+
+ sij <- extract_opportunities(
+ opportunity = opportunity,
+ distance = distance,
+ check_names = TRUE
+ )
+
+ message("Run model")
+ res <- run_law(
+ law = law,
+ mass_origin = mi,
+ mass_destination = mj,
+ distance = distance,
+ opportunity = sij,
+ param = param,
+ check_names = TRUE
+ )
+
+ message("Normalize probabilities")
+ proba_ij <- res$proba
+
+ for(i in proba_ij %>% nrow %>% seq){
+ proba_ij[i,] <- (p2move * proba_ij[i,])/sum(proba_ij[i,])
+ proba_ij[i,i] <- (1 - p2move)
+ }
+
+ proba_ji <- proba_ij
+ colnames(proba_ji) <- colnames(proba_ij) <- rownames(proba_ji) <- rownames(proba_ij) <- row.names(mass)
+
+ for(i in proba_ji %>% nrow %>% seq){
+ proba_ji[,i] <- proba_ji[,i]/sum(proba_ji[,i])
+ }
+
+
+
+ proba <- list(proba_ij = proba_ij,
+ proba_ji = proba_ji)
+ return(proba)
+}
diff --git a/R/build_pts_fun.R b/R/build_pts_fun.R
index f058d60b283fcce42e3eba52f35d5eb8a82a6a15..568370479d6fa535e539a93ef4d9f08b2fa98c5c 100644
--- a/R/build_pts_fun.R
+++ b/R/build_pts_fun.R
@@ -1,20 +1,20 @@
#' Write pts_fun
#'
-#' @description Internal function to write the pts_fun code driving SimInf package to implement deterministic population dynamics of Aedes Albopictus.
+#' @description Internal function to write the pts_fun code driving SimInf package to implement deterministic population dynamics of Aedes Albopictus and synthetic estimates of risk.
#'
-#' @usage build_pts_fun(u0, v0, gdata, diapause_interv)
+#' @usage build_pts_fun(u0, v0, gdata, ldata, diapause_interv)
#'
#' @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 ldata matrix of local data
#' @param diapause_interv interval of days defining favorable period for mosquitoes
#'
#' @return String containing C code used as pts_fun by SimInf package
#'
#' @keywords internal
-#' @noRd
#'
-#' @export
+#' @noRd
build_pts_fun <- function(u0, v0, gdata, ldata, diapause_interv){
@@ -38,29 +38,41 @@ build_pts_fun <- function(u0, v0, gdata, ldata, diapause_interv){
pts_fun <- paste0('
+ //// Declare variables
+
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];
+ //// Population Dynamics
+
// update diapause
- v_new[0] = 0;
+ v_new[', which(rownames(v0) == "z") - 1,'] = 0;
if(',diapause_string,')
- v_new[0] = 1;
+ v_new[', which(rownames(v0) == "z") - 1,'] = 1;
- // update temperature and carrying capacitites
- int i;
- for (i = 1; i < 4; ++i)
- {
- v_new[i] = ldata[2 + (int)ldata[0] * (i-1) + (int)t];
- }
+ // update temperature and carrying capacities
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,'];
+ double fk;
+
+ temperature = ldata[', which(rownames(ldata) == "Tp_1") - 1,' + (int)t];
+ v_new[', which(rownames(v0) == "temperature") - 1,'] = temperature;
+
+ fk = ldata[', which(rownames(ldata) == "kvar") - 1,'] * ldata[', which(rownames(ldata) == "RR7_1") - 1,' + (int)t];
+ if(fk < ldata[', which(rownames(ldata) == "kvar") - 1,']){
+ kP = ldata[', which(rownames(ldata) == "kfix") - 1,'] + fk;
+ kL = ldata[', which(rownames(ldata) == "kfix") - 1,'] + fk;
+ }else{
+ kP = ldata[', which(rownames(ldata) == "kfix") - 1,'] + ldata[', which(rownames(ldata) == "kvar") - 1,'];
+ kL = ldata[', which(rownames(ldata) == "kfix") - 1,'] + ldata[', which(rownames(ldata) == "kvar") - 1,'];
+ }
+
+ v_new[', which(rownames(v0) == "kP") - 1,'] = kP;
+ v_new[', which(rownames(v0) == "kL") - 1,'] = kL;
// Development and mortality functions
@@ -101,7 +113,7 @@ const double * v_0 = &v[-node * nCpmt_v];
// larvae mortality
double fmL;
- fmL = (0.0007 * exp((temperature-10)*0.1838) + 0.02) * (1 + v[', which(rownames(v0) == "Lm") - 1,']/kL);
+ fmL = (', gdata["mu1L"], ' * exp((temperature-10)*', gdata["mu2L"], ') + ', gdata["mu3L"], ') * (1 + v[', which(rownames(v0) == "Lm") - 1,']/kL);
long varLm;
if((fdevL + fmL) > 1){
@@ -120,7 +132,7 @@ const double * v_0 = &v[-node * nCpmt_v];
// pupae mortality
double fmP;
- fmP = (0.0003 * exp((temperature-10)*0.2228) + 0.02) * (1 + v[', which(rownames(v0) == "Pm") - 1,']/kP);
+ fmP = (', gdata["mu1P"], ' * exp((temperature-10)*', gdata["mu2P"], ') + ', gdata["mu3P"], ') * (1 + v[', which(rownames(v0) == "Pm") - 1,']/kP);
long varPm;
if((fdevP + fmP) > 1){
@@ -133,9 +145,9 @@ const double * v_0 = &v[-node * nCpmt_v];
// aldult mortality
double fmA;
- fmA = 0.0003 * exp((temperature - 10)*0.1745) + 0.025;
- if(fmA < ', gdata["muA"], ')
- fmA = ', gdata["muA"], ';
+ fmA = ', gdata["mu1A"], ' * exp((temperature - 10)*', gdata["mu2A"], ') + ', gdata["mu3A"], ';
+ if(fmA < ', gdata["mu3A"], ')
+ fmA = ', gdata["mu3A"], ';
// emerging aldult mortality
double fmAEm;
@@ -204,7 +216,7 @@ const double * v_0 = &v[-node * nCpmt_v];
// parous oviposition-site-seeking adult
v_new[', which(rownames(v0) == "A2om") - 1,'] = v[', which(rownames(v0) == "A2om") - 1,'] + round(varA2om);
- //Calcul of R0
+ //// Calcul of R0
int atot;
atot = u[',which(colnames(u0) == "A1gmI") - 1,
@@ -255,6 +267,10 @@ const double * v_0 = &v[-node * nCpmt_v];
v_new[', which(rownames(v0) == "R0") - 1,'] = (comp_vect*capac_vect)/', gdata["rhoH"],';
}
+ //// Include mobility for infections
+
+ //// pIm infection of host by mosquitoes from another patch
+
// proportion of infected mosquitoes weighted by the probability of being in contact (probability to be in the patch)
// initialize proportion
double wIm;
@@ -262,62 +278,75 @@ const double * v_0 = &v[-node * nCpmt_v];
// 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){
+ // in each node
+ for(j = 0; j <= node; ++j){
+
+ // total number of mosquitoes seeking host
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),'];
+
+ // proportion of infected mosquitoes seeking host
if(nMh > 0){
pIMh = (double) u_0[',which(colnames(u0) == "A2hmI") - 1,' + j * ',ncol(u0),'] / nMh;
} else {
- pIMh = (double) 0;
+ pIMh = 0;
}
- pIMh = pIMh * ldata[2 + (int)ldata[0] * 3 + j];
+ // weight pIMh by the probability to be in j for individual from i
+ pIMh = pIMh * ldata[',rownames(ldata) %>% startsWith(., "Oij_") %>% which %>% min %>% subtract(1),' + 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;
- }
+ // exclude mosquitoes from 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 = wIm - ldata[', which(rownames(ldata) == "pOii") - 1,'] * (u[',which(colnames(u0) == "A2hmI") - 1,'] / nMh);
+ } else {
+ wIm = wIm ;
+ }
v_new[', which(rownames(v0) == "pIm") - 1,'] = (double) wIm;
+ //// pIh infection of mosquitoes by humans from another patch
+
// 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){
+ // in each node
+ for(k = 0; k <= node; ++k){
+
+ // total number of host
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),'];
+
+ // prportion of infected host
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];
+ pIH = pIH * ldata[',rownames(ldata) %>% startsWith(., "Dij_") %>% which %>% min %>% subtract(1),' + k];
wIh = wIh + pIH;
}
-
v_new[', which(rownames(v0) == "pIh") - 1,'] = (double) wIh;
-
return 1;
')
options(scipen=0)
diff --git a/R/build_transitions.R b/R/build_transitions.R
index a70d0aa2fdf7cd5bf97a0276468e62a4ecc82d40..6a1d1a73fc276479ba2ae74f39ee064b38e8e55b 100644
--- a/R/build_transitions.R
+++ b/R/build_transitions.R
@@ -1,8 +1,10 @@
-#' Stochatic transitions
+#' Write stochatic transitions
#'
-#' @description Internal function to build the string vector describing stochastic transitions and required by SimInf package
+#' @description Internal function to build the string vector describing stochastic transitions and counters. This vector is required by SimInf package.
#'
-#' @usage build_transitions()
+#' @usage build_transitions(gdata)
+#'
+#' @param gdata list of parameters
#'
#' @return String vector describing transitions
#'
@@ -12,7 +14,7 @@
#' @export
-build_transitions <- function(){
+build_transitions <- function(gdata){
#### Stochastic transitions ####
## HUMAN ##
@@ -37,12 +39,12 @@ stoch_transitions <- c(
## 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",
+ "A2hmI + Sh -> A2hmI * gammaAh * bMH * Sh/(Sh+Eh+Ih+Rh) * pOii -> A2gmI + Eh + ninfh",
# Infection of host in another patch
"Sh -> 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",
+ "A2hmI -> A2hmI * gammaAh * (1 - bMH * (Sh/(Sh+Eh+Ih+Rh) * pOii + (1 - pOii))) -> A2gmI",
"Eh -> Eh * muH -> Ih",
"Ih -> Ih * rhoH -> Rh",
@@ -54,16 +56,16 @@ stoch_transitions <- c(
"@ -> (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
+ paste0("A1omI -> A1omI * gammaAo -> A2hmI + ",gdata[["beta1"]]," * Neggs"), ## New eggs
"A2gmI -> temperature > TAG ? A2gmI * (temperature - TAG) / TDDAG : 0 -> A2omI",
- "A2omI -> A2omI * gammaAo -> A2hmI + 75 * Neggs", ## New eggs
+ paste0("A2omI -> A2omI * gammaAo -> A2hmI + ",gdata[["beta2"]]," * Neggs"), ## New eggs
# Mortality
- "A1gmI -> 0.0003 * exp((temperature-10)*0.1745) < 0 ? A1gmI * muA : A1gmI * (0.0003 * exp((temperature-10) * 0.1745) + muA) -> @",
- "A1omI -> 0.0003 * exp((temperature-10)*0.1745) < 0 ? A1omI * (muA + muR) : A1omI * (0.0003 * exp((temperature-10) * 0.1745) + muA + muR) -> @",
- "A2hmI -> 0.0003 * exp((temperature-10)*0.1745) < 0 ? A2hmI * (muA + muR) : A2hmI * (0.0003 * exp((temperature-10) * 0.1745) + muA + muR) -> @",
- "A2gmI -> 0.0003 * exp((temperature-10)*0.1745) < 0 ? A2gmI * muA : A2gmI * (0.0003 * exp((temperature-10) * 0.1745) + muA) -> @",
- "A2omI -> 0.0003 * exp((temperature-10)*0.1745) < 0 ? A2omI * (muA + muR) : A2omI * (0.0003 * exp((temperature-10) * 0.1745) + muA + muR) -> @"
+ "A1gmI -> mu1A * exp((temperature-10)*mu2A) < 0 ? A1gmI * mu3A : A1gmI * (mu1A * exp((temperature-10) * mu2A) + mu3A) -> @",
+ "A1omI -> mu1A * exp((temperature-10)*mu2A) < 0 ? A1omI * (mu3A + muR) : A1omI * (mu1A * exp((temperature-10) * mu2A) + mu3A + muR) -> @",
+ "A2hmI -> mu1A * exp((temperature-10)*mu2A) < 0 ? A2hmI * (mu3A + muR) : A2hmI * (mu1A * exp((temperature-10) * mu2A) + mu3A + muR) -> @",
+ "A2gmI -> mu1A * exp((temperature-10)*mu2A) < 0 ? A2gmI * mu3A : A2gmI * (mu1A * exp((temperature-10) * mu2A) + mu3A) -> @",
+ "A2omI -> mu1A * exp((temperature-10)*mu2A) < 0 ? A2omI * (mu3A + muR) : A2omI * (mu1A * exp((temperature-10) * mu2A) + mu3A + muR) -> @"
)
return(stoch_transitions)
diff --git a/R/check_PARCELLE.R b/R/check_PARCELLE.R
new file mode 100644
index 0000000000000000000000000000000000000000..68dfa7722ac2ce2e8cb31c3a61eea4f66fbb0f88
--- /dev/null
+++ b/R/check_PARCELLE.R
@@ -0,0 +1,27 @@
+#' Check PARCELLE object
+#'
+#' @description function to check format of PARCELLE object
+#'
+#' @param PARCELLE percelle object (data.frame or data.table)
+#'
+#'
+#' @noRd
+
+
+check_PARCELLE <- function(PARCELLE){
+
+ 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"')
+
+}
diff --git a/R/check_gdata.R b/R/check_gdata.R
index c0b639fc952cb992eeea9a443c118dc41a900e12..aa5d3543310d84bf59a1f78496a5d6c30527390d 100644
--- a/R/check_gdata.R
+++ b/R/check_gdata.R
@@ -2,29 +2,31 @@
#'
#' @description function to build gdata
#'
-#' @usage check_gdata(gdata, vector_species)
-#'
#' @param gdata list
-#' @param vector_species string
+#' @param vector_species string. Default is "Ae. albopictus" in "temperate" climate. Unique dataset implemented to date.
+#' @param climat string. Default is "Ae. albopictus" in "temperate" climate. Unique dataset implemented to date.
#'
#' @return gdata
+#'
+#' @noRd
+
-check_gdata <- function(gdata, vector_species, verbose = T){
+check_gdata <- function(gdata, vector_species, climate, 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",
+ if(vector_species == "Ae. albopictus" & climate == "temperate"){
+ gdata <- build_gdata(vector_species = vector_species,
+ climate = climate,
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"))
+ stop(paste("Parameters for", vector_species, "in", climate, " climate 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
+ "bMH", "muH", "rhoH", "muE", "TE", "TDDE", "mu1L", "mu2L", "mu3L", "q1L", "q2L", "q3L", "muEM", "mu1P", "mu2P", "mu3P", "q1P", "q2P", "q3P", "mu1A", "mu2A", "mu3A", "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/cleanmeteo.R b/R/cleanmeteo.R
index 93cbefaf2a5d3cb1ef2727c2d8fb1324cf63a227..4020ef9597dc0be6b7dfa9ba4d0711723d548072 100644
--- a/R/cleanmeteo.R
+++ b/R/cleanmeteo.R
@@ -12,6 +12,9 @@
#' @import data.table
#'
#' @export
+#'
+#' @noRd
+
# Clean meteorological data
diff --git a/R/create_ldata.R b/R/create_ldata.R
index 2eca89517e3de5f691988231cdfea9f8e1226943..26fc1ac3ad6e48b04cdc5b37b91e9d3d23665c08 100644
--- a/R/create_ldata.R
+++ b/R/create_ldata.R
@@ -2,14 +2,16 @@
#'
#' @description Function to define local data
#'
-#' @usage create_ldata(PARCELLE, METEO, gdata, time_max, nYR_init, nodeID)
+#' @usage create_ldata(PARCELLE, METEO, gdata, mMov, TS_sim)
#'
#' @param PARCELLE 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
+#' @param mMov list of two matrices. Probabilities of location of individuals.
+#' proba_ij is the probabilities of individuals from i (columns) to be bitten in j (rows).
+#' proba_ji is the probabilities for mosquito in i (rows) to bite an individual from j (columns).
+#' Row sums should be equal to 1.
+#' @param TS_sim list. four vectors of simulated dates (dates) and simulated days (numeric) with and without the initialization period
#'
#' @importFrom pbapply pbapply
#' @importFrom parallel makeCluster
@@ -17,15 +19,18 @@
#' @importFrom parallel clusterExport
#' @importFrom parallel clusterEvalQ
#' @importFrom parallel stopCluster
+#' @importFrom magrittr is_less_than
+#' @importFrom zoo rollapply
#'
#' @return matrix
#'
#' @keywords ldata METEO
#'
-#' @export
+#' @noRd
-create_ldata <- function(PARCELLE, METEO, gdata, mMov, time_max, nYR_init = 1){
+
+create_ldata <- function(PARCELLE, METEO, gdata, mMov, TS_sim){
## 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.
@@ -39,95 +44,83 @@ create_ldata <- function(PARCELLE, METEO, gdata, mMov, time_max, nYR_init = 1){
## CHECKS
- try(if(NA %in% match(PARCELLE$STATION, METEO$ID))
- stop("Some stations are missing from the METEOrological dataset"))
+ if(NA %in% match(PARCELLE$STATION, METEO$ID))
+ stop("Some stations are missing from the meteorological dataset")
- 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"))
+ 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[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"))
+ if(FALSE %in% (as.IDate(TS_sim$time_serie_output_d) %in% METEO[ID %in% PARCELLE$STATION, DATE]))
+ stop("meteorological data are required for each day of the simulated period")
## CALCULATION
## for each administrative unit
message("## Generating parameters for all patches can take time, please be patient. ##")
- #make it parallel
+ # Normalize precipitation
+ METEO[,RR := RR/max(RR)]
- clust <- detectCores() %>% makeCluster
- #export variables
- 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"))
+ t_sim <- length(TS_sim$time_serie_date)
-ldata <- pbapply(PARCELLE, 1, function(x){
+ # Initialization period
+ METEO <- expand.grid(DATE = TS_sim$time_serie_date,
+ ID = METEO$ID %>% unique) %>% merge(., METEO, all.x = T)
+ setDT(METEO)
- # Daily temperature corrected by the altitude
- 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[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
- kPfix <- x["KLfix"] %>% as.numeric
- kPvar <- x["KLvar"] %>% as.numeric
-
- # Daily carrying capacities corrected by the rainfall amount
- kL <- sapply(raincumul7 %>% length %>% seq, function(day){
- if(kLvar > (kLvar * (raincumul7[day]/100)))
- kLfix + kLvar * (raincumul7[day]/100)
- else
- kLfix + kLvar
- })
+ ldata <- pbapply(PARCELLE, 1, function(x){
- kP <- sapply(raincumul7 %>% length %>% seq, function(day){
- if(kPvar > (kPvar * (raincumul7[day]/100)))
- kPfix + kPvar * (raincumul7[day]/100)
- else
- kPfix+ kPvar
- })
+ while(TRUE %in% METEO[ID == x["STATION"] %>% as.numeric, is.na(RR)]){
+ NA_start <- METEO[ID == x["STATION"] %>% as.numeric & is.na(RR), DATE %>% as.IDate] %>% min %>% format(., "%m-%d")
+ NA_length <- METEO[ID == x["STATION"] %>% as.numeric & is.na(RR), DATE] %>% length
+
+ NA_start <- which(METEO[ID == x["STATION"] %>% as.numeric & !is.na(RR), DATE %>% as.IDate %>% format(., "%m-%d")] == NA_start) %>% min
- ## data for population dynamics initialization
- if((nYR_init * 365) > time_max){
- temp_init <- rep(temperature[seq(365)], nYR_init)
- kL_init <- rep(kL[seq(365)], nYR_init)
- kP_init <- rep(kP[seq(365)], nYR_init)
- } else {
- temp_init <- temperature[seq(nYR_init * 365)]
- kL_init <- kL[seq(nYR_init * 365)]
- kP_init <- kP[seq(nYR_init * 365)]
- }
+ METEO[ID == x["STATION"] %>% as.numeric & is.na(RR), c("RR", "TP")] <-
+ METEO[ID == x["STATION"] %>% as.numeric & !is.na(RR), c("RR", "TP")][seq(NA_start:(NA_start+(NA_length-1))),]
+ }
- c(time_max + 365 * nYR_init,
- # Temperature
- temp_init, temperature,
- # Carrying capacities related to rainfall
- kL_init, kL,
- kP_init, kP
- ) %>% as.numeric
+ # Daily temperature corrected by the altitude '(4.2/1000)' is in Ocelet code
+ METEO[ID == x["STATION"] %>% as.numeric, TP := TP - as.numeric(x["DIFF_ALT"])*(4.2/1000)]
-}, cl = clust)
+ temperature <- METEO[ID == x["STATION"] %>% as.numeric, TP]
-stopCluster(clust)
+ # Cumulative rainfall over 7 days
+ raincumul7 <- rollapply(c(rep(0,6), METEO[ID == x["STATION"] %>% as.numeric, RR]), 7, sum)
+
+ # Simplify with unique carrying capacity value
+ kfix <- x["KLfix"] %>% as.numeric
+ kvar <- x["KLvar"] %>% as.numeric
+
+ c(t_sim,
+ # carrying capacities
+ kfix, kvar,
+ # Temperature
+ temperature,
+ # Carrying capacities related to rainfall
+ raincumul7
+ ) %>% as.numeric
+
+ })
-colnames(ldata) <- PARCELLE$ID
+ colnames(ldata) <- PARCELLE$ID
-ldata %<>% .[c(1, seq(nrow(ldata))),]
+ ldata %<>% .[c(1,1,seq(nrow(ldata))),]
-ldata[2,names(diag(mMov))] <- diag(mMov)
+ ldata[2,] <- diag(mMov$proba_ij)
+ ldata[3,] <- diag(mMov$proba_ji)
-ldata %<>% rbind(., mMov[,colnames(ldata)])
+ ldata %<>% rbind(., t(mMov$proba_ij[,colnames(ldata)]))
+ ldata %<>% rbind(., t(mMov$proba_ji[,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)
+ rownames(ldata) <- c("n_days", "pOii","pDii",
+ "kfix", "kvar",
+ paste0("Tp_", seq(t_sim)),
+ paste0("RR7_", seq(t_sim)),
+ paste0("Oij_", PARCELLE$ID),
+ paste0("Dij_", PARCELLE$ID)
+ )
-return(ldata)
+ return(ldata)
}
diff --git a/R/diapause.R b/R/diapause.R
index 88f2fba9f3420cf9ab432cb096383e299abfbf09..6e0fd24ff38288a608704bf3b67aaa661d94b0b6 100644
--- a/R/diapause.R
+++ b/R/diapause.R
@@ -5,31 +5,24 @@
#' @usage diapause(startFav, endFav, time_max, nYR_init)
#'
#' @param startFav Date. First day of the favourable period.
-#' @param endFav Date. Last day of the favourable period.
-#' @param time_max Numerical value. Simulation length (days).
-#' @param nYR_init Numerical value. Number of years for initialization of the dynamics (default is 2).
+#' @param endFav Date. Last day of the favorable period.
+#' @param TS_sim list. four vectors of simulated dates (dates) and simulated days (numeric) with and without the initialization period
#'
-#' @return Numeric vector with days of start and end of the favourable period
+#' @return Numeric vector with days of start and end of the favorable period
#'
#' @keywords diapause
#'
-#' @export
-
-
-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 %>% yday
+#' @importFrom magrittr subtract
+#' @importFrom magrittr %>%
+#' @importFrom data.table as.IDate
+#'
+#' @noRd
-DurFav <- endFav %>% yday %>% subtract(x)
-DurUnFav <- 365 - DurFav
+diapause <- function(startFav, endFav, TS_sim){
-while(max(x) < (time_max + nYR_init * 365)) {
- x <- c(x, max(x) + DurFav)
- x <- c(x, max(x) + DurUnFav)
-}
+ startFav <- which(format(TS_sim$time_serie_date, "%m-%d") == format(as.IDate(startFav), "%m-%d"))
+ endFav <- which(format(TS_sim$time_serie_date, "%m-%d") == format(as.IDate(endFav), "%m-%d"))
+ x <- c(startFav, endFav) %>% .[order(.)]
return(x)
}
diff --git a/R/list_compartments.R b/R/list_compartments.R
new file mode 100644
index 0000000000000000000000000000000000000000..09f270f840bf5a522c02f2aee3389e173e44fa0e
--- /dev/null
+++ b/R/list_compartments.R
@@ -0,0 +1,46 @@
+#' List compartments
+#'
+#' @description function to list compartments
+#'
+#' @usage list_compartments()
+#'
+#' @return u0 and v0 compartment as list
+#'
+#' @noRd
+
+
+list_compartments <- function(){
+
+u0_compartments <- c(
+ # Human
+ "Sh",
+ "Eh",
+ "Ih",
+ "Rh",
+
+ # Infected mosquitoes (non infected mosquitoes will be simulated into v0)
+ "A1gmI",
+ "A1omI",
+ "A2hmI",
+ "A2gmI",
+ "A2omI",
+ "Neggs", #new eggs from infected mosquitoes
+
+ # Number of autochtonous human infection
+ "ninfh",
+ "ninfm1",
+ "ninfm2"
+)
+
+
+v0_compartments <- c(
+ "z", "temperature", "kL", "kP",
+ "Em", "Lm", "Pm",
+ "Aemm", "A1hm", "A1gm", "A1om", "A2hm", "A2gm", "A2om",
+ "prevEggs", "nIm1", "nIm2",
+ "R0", "pIh", "pIm")
+
+return(list(u0_compartments = u0_compartments,
+ v0_compartments = v0_compartments))
+
+}
diff --git a/R/runArboRisk.R b/R/runArboRisk.R
index 4196fea543042029f163194215e8e553f2dd737b..95ae9b4b91cf21cf4366d2b2219a0c7b6cab585d 100644
--- a/R/runArboRisk.R
+++ b/R/runArboRisk.R
@@ -11,15 +11,24 @@
#' @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 nYR_init numeric. Number of years used to initialize the population dynamics (default is 1)
#' @param n_sim integer
#' @param introduction_pts data.frame or data.table
+#' @param u0 matrix Initial population state (patches as columns and in rows: Sh, Eh, Ih, Rh, A1gmI, A1omI, A2hmI, A2gmI, A2omI, Neggs, ninfh, ninfm1, ninfm2)
+#' @param v0 matrix Initial population and time-dependant parameters state (patches as rows and in columns:
+#' z: diapause (0 = dipause, 1 = favorable period); temperature; kL and kP: carrying capacities for larvae and pupae (rainfall dependent);
+#' Em: number of eggs ; Lm: number of larvae ; Pm: number of pupae ; Aemm: number of emerging adults ;
+#' A1hm: number of nulliparous adults seeking for host ; A1gm: number of gorged nulliparous adults ; A1om: number of nulliparous adults seeking for oviposition sites ;
+#' A2hm: number of parous adults seeking for host ; A2gm: number of gorged parous adults ; A1om: number of parous adults seeking for oviposition sites ;
+#' prevEggs, nIm1, nIm2, R0, pIh and pIm are continous variables calculated over time)
+#' @param initMosq numeric. Number of eggs in each patch at t0 (default is 100000).
#'
#'
#' @importFrom SimInf mparse
#' @importFrom SimInf run
#' @importFrom SimInf trajectory
#' @importFrom pbapply pblapply
+#' @importFrom magrittr is_greater_than
#'
#' @return data.table
#'
@@ -27,16 +36,25 @@
runArboRisk <- function(PARCELLE,
+
vector_species = "Ae. albopictus",
+ climate = "temperate",
gdata = NULL,
+
ldata = NULL,
+
METEO = NULL,
mMov = NULL,
+
start_date = NULL,
end_date = NULL,
nYR_init = 1,
n_sim = 1,
- introduction_pts = NULL){
+ introduction_pts = NULL,
+
+ u0 = NULL,
+ v0 = NULL,
+ initMosq = 100000){
if(is.null(ldata) & is.null(METEO)) stop("you must provide either ldata matrix or meteorological dataset (METEO)")
@@ -54,26 +72,7 @@ runArboRisk <- function(PARCELLE,
# 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
+ check_PARCELLE(PARCELLE)
############################
### Model initialization ###
@@ -89,9 +88,8 @@ runArboRisk <- function(PARCELLE,
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
+ ## Simulation period ##
+ TS_sim <- sim_time_serie(start_date, end_date, nYR_init)
###########################################
## Create u0, v0 and compartment objects ##
@@ -99,9 +97,11 @@ runArboRisk <- function(PARCELLE,
message("## Initialization ##")
- init <- iniState(PARCELLE, initMosq = 100000)
+ init <- iniState(PARCELLE, initMosq = initMosq)
+ if(is.null(u0))
u0 <- init$u0
+ if(is.null(v0))
v0 <- init$v0
###############################################################
@@ -111,7 +111,7 @@ runArboRisk <- function(PARCELLE,
#### Global data (general parameters) ####
- gdata <- check_gdata(gdata, vector_species)
+ gdata <- check_gdata(gdata, vector_species, climate)
##################
## Define ldata ##
@@ -127,13 +127,14 @@ runArboRisk <- function(PARCELLE,
METEO,
gdata,
vector_species,
+ climate,
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(ldata[1,1] != length(TS_sim$time_serie_date))
+ 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 (",length(TS_sim$time_serie_date)," 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),")."))
@@ -147,8 +148,7 @@ runArboRisk <- function(PARCELLE,
if( !is.null(gdata$startFav) & !is.null(gdata$endFav)){
diapause_interv <- diapause(startFav = gdata$startFav,
endFav = gdata$endFav,
- time_max,
- nYR_init)
+ TS_sim = TS_sim)
}
gdata %<>% .[-which(names(gdata) %in% c("startFav","endFav"))]
@@ -158,7 +158,7 @@ runArboRisk <- function(PARCELLE,
message("## Build stochastic transitions ##")
#### Stochastic transitions ####
- stoch_transitions <- build_transitions()
+ stoch_transitions <- build_transitions(gdata)
###########################################################
### Deterministic transitions for mosquitoes life cycle ###
@@ -190,7 +190,7 @@ runArboRisk <- function(PARCELLE,
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.
+ tspan = TS_sim$time_serie_num, # number of days simulated. Start at 0 to get the indexing correct in the post-time-step function.
pts_fun = pts_fun
)
@@ -219,12 +219,14 @@ runArboRisk <- function(PARCELLE,
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"))
+
+ if(FALSE %in% (introduction_pts$time %in% seq.Date(from = start_date, to = end_date, by ="day")))
+ warning(paste0("At least one introduction occurs outside the simulation period _n", capture.output(introduction_pts[!introduction_pts$time %in% seq.Date(from = start_date, to = end_date, by ="day") ]), collapse = "\n"))
+
+ introduction_pts$time %<>% match(., TS_sim$time_serie_date)
# turn IDs into numerical IDs
- introduction_pts$node %<>% match(., PARCELLE$ID)
+ introduction_pts$node <- introduction_pts$dest %<>% match(., PARCELLE$ID)
# select the stage
introduction_pts$select %<>% match(., u0 %>% colnames)
@@ -247,7 +249,7 @@ runArboRisk <- function(PARCELLE,
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.
+ tspan = TS_sim$time_serie_num, # 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,
@@ -262,12 +264,10 @@ runArboRisk <- function(PARCELLE,
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 := TS_sim$time_serie_date, by = node]
- traj %<>% .[DATE >= start_date]
- traj[,time := time - min(time)]
+ traj %<>% .[DATE %in% TS_sim$time_serie_output_d]
+ traj[,time := TS_sim$time_serie_output_t, by = node]
traj[, ID := PARCELLE$ID[node]]
diff --git a/R/sim_time_serie.R b/R/sim_time_serie.R
new file mode 100644
index 0000000000000000000000000000000000000000..3a80979cfa98a6f4acfe96f7024bf9acd2ffc620
--- /dev/null
+++ b/R/sim_time_serie.R
@@ -0,0 +1,42 @@
+#' Time serie of the simulated periods
+#'
+#' @description Function to calculates the dates and time points for simulations
+#'
+#' @usage sim_time_serie(start_date, end_date, nYR_init)
+#'
+#' @param start_date Date. First day of output desired.
+#' @param end_date Date. Last day of output desired.
+#' @param nYR_init Numerical value. Number of years for initialization of the dynamics.
+#'
+#' @return list with dates and numeric vectors of simulated days
+#'
+#' @keywords time serie
+#'
+#' @importFrom magrittr subtract
+#' @importFrom magrittr %>%
+#' @importFrom data.table as.IDate
+#'
+#' @noRd
+
+sim_time_serie <- function(start_date, end_date, nYR_init){
+
+ time_serie_output_d <- start_date:end_date
+
+ time_serie_output_t <- seq(1:length(time_serie_output_d))
+
+ year_start <- start_date %>% format(., format="%Y") %>% as.numeric %>% subtract(nYR_init) %>% format(., format="%Y")
+ start_date <- year_start %>% paste0(., "-", start_date %>% format(., format="%m-%d")) %>% as.IDate
+
+ time_serie_date <- start_date : (end_date %>% as.IDate)
+
+ time_serie_num <- c(0, seq(length(time_serie_date)-1))
+
+ return(list(time_serie_date = time_serie_date %>% as.IDate,
+ time_serie_num = time_serie_num,
+ time_serie_output_d = time_serie_output_d %>% as.IDate,
+ time_serie_output_t = time_serie_output_t)
+
+ )
+}
+
+
diff --git a/data-raw/DATASET.R b/data-raw/DATASET.R
index ee906c4f29edc7bcb7fdc4b867c43d506b2f3475..42f92ed4fa0aba1d452f1bce502d5e5f6c70f4a9 100644
--- a/data-raw/DATASET.R
+++ b/data-raw/DATASET.R
@@ -1,8 +1,6 @@
#' Parcelle and meteo data
#'
-#' Useful metadata about airports.
-#'
-#' @name Parcelle and meteo data
+#' @name PARCELLE and METEO data
#' @docType data
#' @format two data frame with columns:
#' \describe{
diff --git a/data/SpatVec.cpg b/data/SpatVec.cpg
new file mode 100644
index 0000000000000000000000000000000000000000..3ad133c048f2189041151425a73485649e6c32c0
--- /dev/null
+++ b/data/SpatVec.cpg
@@ -0,0 +1 @@
+UTF-8
\ No newline at end of file
diff --git a/data/SpatVec.dbf b/data/SpatVec.dbf
new file mode 100644
index 0000000000000000000000000000000000000000..c16292f0335e83b91a143ba4f6803254e113d6fc
Binary files /dev/null and b/data/SpatVec.dbf differ
diff --git a/data/SpatVec.prj b/data/SpatVec.prj
new file mode 100644
index 0000000000000000000000000000000000000000..ae0206b68de2ed81139b89a08ddd36a6b0ed7e35
--- /dev/null
+++ b/data/SpatVec.prj
@@ -0,0 +1 @@
+PROJCS["RGF_1993_Lambert_93",GEOGCS["GCS_RGF_1993",DATUM["D_RGF_1993",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",700000.0],PARAMETER["False_Northing",6600000.0],PARAMETER["Central_Meridian",3.0],PARAMETER["Standard_Parallel_1",49.0],PARAMETER["Standard_Parallel_2",44.0],PARAMETER["Latitude_Of_Origin",46.5],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/data/SpatVec.shp b/data/SpatVec.shp
new file mode 100644
index 0000000000000000000000000000000000000000..29c9dcf0af942000cc41f9830edc414328939713
Binary files /dev/null and b/data/SpatVec.shp differ
diff --git a/data/SpatVec.shx b/data/SpatVec.shx
new file mode 100644
index 0000000000000000000000000000000000000000..91f466fd64db47f061303adff6e7f20504d75333
Binary files /dev/null and b/data/SpatVec.shx differ
diff --git a/man/build_E_random.Rd b/man/build_E_random.Rd
index cfe002588f6e826abd488ee12318e19dfffc9928..dd2380c6e053ce8b02f1aaeab5afeacdb4c57240 100644
--- a/man/build_E_random.Rd
+++ b/man/build_E_random.Rd
@@ -20,10 +20,12 @@ build_E_random(period_start, period_end, n_ind = 1, n_events = 1, stage = "Eh",
\item{loc}{String or numerical vector. List of patch to seed random introductions.}
}
\value{
-events matrix with scheduled introductions.
+events matrix with scheduled introductions
}
\description{
-Function used to generate the event matrix and associated objects required by SimInf package
+Function used to generate the random introduction.
+The introduction events generated can occurs in different locations and at different dates in a given period.
+The number of individuals introduced and their epidemiological stages can be deterministic (unique value provided) or random (randomly selected in a vector of values).
}
\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])
diff --git a/man/build_gdata.Rd b/man/build_gdata.Rd
index 4f7e15faa979d231280527b15a0231565fe8ab7b..104b562a48fe0e40a1b52dbd022c7a4be90c39d0 100644
--- a/man/build_gdata.Rd
+++ b/man/build_gdata.Rd
@@ -2,18 +2,90 @@
% Please edit documentation in R/build_gdata.R
\name{build_gdata}
\alias{build_gdata}
-\title{Generate global parameters dataset}
+\title{Generate global parameters list}
\usage{
build_gdata()
}
\arguments{
-\item{vector_species}{}
+\item{vector_species}{string. Default is "Ae. albopictus" in "temperate" climate. Unique dataset implemented to date.}
+
+\item{bHM}{numeric. Probability of infection from host to vector when an infected human is bitten by an susceptible mosquito.}
+
+\item{bMH}{numeric. Probability of infection from vector to host when a human is bitten by an infected mosquito.}
+
+\item{muH}{numeric. Daily transition rate from exposed (E) to infected (I) for humans (1/days)}
+
+\item{rhoH}{numeric. Daily recovery rate of humans (1/days)}
+
+\item{muE}{numeric. Daily mortality rate of eggs (1/days)}
+
+\item{TE}{numeric. Minimal temperature needed for egg development (°C)}
+
+\item{TDDE}{numeric. Total number of degree-day necessary for egg development (°C)}
+
+\item{mu1L}{numeric. Parameter for the larvae mortality function: \eqn{ mu1L \times e^{(temperature_{t}-10) \times mu2L} + mu3L}}
+
+\item{mu2L}{numeric. Parameter for the larvae mortality function: \eqn{ mu1L \times e^{(temperature_{t}-10) \times mu2L} + mu3L}}
+
+\item{mu3L}{numeric. Parameter for the larvae mortality function: \eqn{ mu1L \times e^{(temperature_{t}-10) \times mu2L} + mu3L}}
+
+\item{q1L}{numeric. Parameter for the function of transition from larva to pupa \eqn{q1L \times temperature_{t}^{2} + q2L \times temperature_{t} + q3L}}
+
+\item{q2L}{numeric. Parameter for the function of transition from larva to pupa \eqn{q1L \times temperature_{t}^{2} + q2L \times temperature_{t} + q3L}}
+
+\item{q3L}{numeric. Parameter for the function of transition from larva to pupa \eqn{q1L \times temperature_{t}^{2} + q2L \times temperature_{t} + q3L}}
+
+\item{muEM}{numeric. Daily mortality rate of emerging adults during the emergence (1/days)}
+
+\item{mu1P}{numeric. Parameter for the larvae mortality function: \eqn{mu1P \times e^{(temperature_{t}-10) \times mu2P} + mu3P}}
+
+\item{mu2P}{numeric. Parameter for the larvae mortality function: \eqn{mu1P \times e^{(temperature_{t}-10) \times mu2P} + mu3P}}
+
+\item{mu3P}{numeric. Parameter for the larvae mortality function: \eqn{mu1P \times e^{(temperature_{t}-10) \times mu2P} + mu3P}}
+
+\item{q1P}{numeric. Parameter for the function of transition from pupae to emerging adult \eqn{q1P \times temperature_{t}^{2} + q2P \times temperature_{t} + q3P}}
+
+\item{q2P}{numeric. Parameter for the function of transition from pupae to emerging adult \eqn{q1P \times temperature_{t}^{2} + q2P \times temperature_{t} + q3P}}
+
+\item{q3P}{numeric. Parameter for the function of transition from pupae to emerging adult \eqn{q1P \times temperature_{t}^{2} + q2P \times temperature_{t} + q3P}}
+
+\item{mu1A}{numeric. Parameter for the adult mortality function: \eqn{mu1A \times e^{(temperature_{t}-10) \times mu2A} + mu3A}}
+
+\item{mu2A}{numeric. Parameter for the adult mortality function: \eqn{mu1A \times e^{(temperature_{t}-10) \times mu2A} + mu3A}}
+
+\item{mu3A}{numeric. Parameter for the adult mortality function: \eqn{mu1A \times e^{(temperature_{t}-10) \times mu2A} + mu3A}}
+
+\item{gammaAem}{numeric. Daily development rate of emerging adults (1/days)}
+
+\item{sigma}{numeric. Sex-ratio at the emergence}
+
+\item{gammaAh}{numeric. Daily transition rate from host-seeking to engorged adults (1/days)}
+
+\item{muR}{numeric. Daily mortality rate related to seeking behavior (1/days)}
+
+\item{TAG}{numeric. Minimal temperature needed for egg maturation (°C)}
+
+\item{TDDAG}{numeric. Total number of degree-days necessary for egg maturation (°C)}
+
+\item{gammaAo}{numeric. Daily transition rate from oviposition site-seeking to host-seeking adults (1/days)}
+
+\item{beta1}{numeric. Number of eggs laid by ovipositing nulliparous females (per female)}
+
+\item{beta2}{numeric. Number of eggs laid by ovipositing parous females (per female)}
+
+\item{startFav}{date. First day of the favorable period for the mosquito (depend on the species and the environment).}
+
+\item{endFav}{date. Last day of the favorable period for the mosquito (depend on the species and the environment).}
+
+\item{verbose}{logical. Provide information on parameters during generation}
+
+\item{climat}{string. Default is "Ae. albopictus" in "temperate" climate. Unique dataset implemented to date.}
}
\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.
+The functions formats all global parameters required for the model.
}
\examples{
build_gdata()
diff --git a/man/build_ldata.Rd b/man/build_ldata.Rd
index b3add023d68c46e455d2335fa2b9fd05cd5c1c82..1118bffa4090e2efaa1274a6be72c1c2ffa294f4 100644
--- a/man/build_ldata.Rd
+++ b/man/build_ldata.Rd
@@ -2,30 +2,35 @@
% Please edit documentation in R/build_ldata.R
\name{build_ldata}
\alias{build_ldata}
-\title{Build ldata}
+\title{Generate local parameters matrix}
\usage{
build_ldata(PARCELLE, METEO)
}
\arguments{
-\item{PARCELLE}{data.frame or data.table}
+\item{PARCELLE}{data.frame or data.table describing the patches. Required columns: "ID": unique identifier, "POP": population size, "SURF_HA": surface in hectare, "KLfix": fix carrying capacity for larvae, "KLvar": variable carrying capacity for larvae, "KPfix": fix carrying capacity for pupae, "KPvar": variable carrying capacity for pupae, "STATION": meteorological station identifier, "DIFF_ALT": difference between meteorological station altitude and average altitude of the patch.}
-\item{METEO}{data.frame or data.table}
+\item{METEO}{data.frame or data.table reporting the daily meteorological data for each meteorological station. Required columns: 'ÌD', 'DATE', 'RR': daily precipitation (mm), 'TP': daily mean temperature (°C)}
-\item{gdata}{list}
+\item{gdata}{list of parameters. Can be generated using `build_gdata` function}
-\item{vector_species}{string}
+\item{vector_species}{string. Default is "Ae. albopictus" in "temperate" climate. Unique dataset implemented to date. Only needed if gdata is not provided.}
-\item{mMov}{double matrix}
+\item{mMov}{list of two matrices. Probabilities of location of individuals.
+proba_ij is the probabilities of individuals from i (columns) to be bitten in j (rows).
+proba_ji is the probabilities for mosquito in i (rows) to bite an individual from j (columns).
+Row sums should be equal to 1.}
-\item{start_date}{date in '\%Y-\%m-\%d' format}
+\item{start_date}{date in '\%Y-\%m-\%d' format. Starting date of the simulation. If not provided, the function uses the oldest date in the meteorological dataset.}
-\item{end_date}{date in '\%Y-\%m-\%d' format}
+\item{end_date}{date in '\%Y-\%m-\%d' format. Ending date of the simulation. If not provided, the function uses the most recent date in the meteorological dataset.}
-\item{nYR_init}{numerical value. Number of years used to initialize the population dynamics (default is 1)}
+\item{nYR_init}{numerical value. Number of years used to initialize the population dynamics (default is 1).}
+
+\item{climat}{string. Default is "Ae. albopictus" in "temperate" climate. Unique dataset implemented to date. Only needed if gdata is not provided.}
}
\value{
matrix
}
\description{
-function to build ldata
+The function computes and formats for each pacth the list of local parameters required by the model. Arguments related to the simulation are required to compile daily meteorological data for the simulated period.
}
diff --git a/man/build_mMov.Rd b/man/build_mMov.Rd
index 122878197b11d29074af73323b03ba904bbd6af4..d8f785e74678422c2d728530ef6812d7c20de87f 100644
--- a/man/build_mMov.Rd
+++ b/man/build_mMov.Rd
@@ -1,26 +1,26 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/build_mMov.R
+% Please edit documentation in R/build_mMov_TDLM_based.R
\name{build_mMov}
\alias{build_mMov}
-\title{Creates a synthetic network}
+\title{Build matrix of contact probabilities}
\usage{
-build_mMov()
+build_mMov(SpatVec)
}
\arguments{
-\item{PARCELLE}{data.table}
+\item{SpatVec}{SpatVector or sf. Spatial vector of polygons representing patches. Polygons must include 'POP' and 'ID' attributes.}
-\item{group}{string}
+\item{law}{string. Law used to calculate probabilities. Default is "Unif". See documentation of TDLM::run_law function for law options details.}
-\item{within_clust_lev}{numeric. Probability of creating links between wards of the same building}
+\item{param}{numeric. Parameter used to adjust the importance of distance or opportunity associated with the chosen law.
+A single value or a vector of several parameter values can be used.
+Not necessary for the original radiation law or the uniform law. (see TDLM package for more details)}
-\item{between_clust_lev}{numeric. Probability of creating links between wards of different building}
+\item{p2move}{numeric. Probability to be in the residential patch.}
-\item{clust_ratio_inout}{numeric. Ratio of intra-building clustering (relative to inter-building clustering)}
-
-\item{verbose}{logical. Print plot if TRUE}
+\item{pCom}{numeric. Proportion of cummuters per patch}
}
\value{
-a matrix with random probabilities of contact
+list of normalized matrices for probabilities of contact of agent from patch i with agent of patch j in j (proba_ij) and probabilities of contact of agent from patch j with agent of patch i in i (proba_ji)
}
\description{
This function creates a synthetic human mobility network.
diff --git a/man/build_mMov_erdosrenyigame.Rd b/man/build_mMov_erdosrenyigame.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..38e3d1ab3a7c6863e7b46314d492e9aa24d2bfde
--- /dev/null
+++ b/man/build_mMov_erdosrenyigame.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/build_mMov.R
+\name{build_mMov_erdosrenyigame}
+\alias{build_mMov_erdosrenyigame}
+\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
deleted file mode 100644
index f33918c644e8375feee4281df4fc869097be910b..0000000000000000000000000000000000000000
--- a/man/check_gdata.Rd
+++ /dev/null
@@ -1,19 +0,0 @@
-% 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/cleanmeteo.Rd b/man/cleanmeteo.Rd
deleted file mode 100644
index bad62c93e148e2f4b28d3c822d84d2368ea1e312..0000000000000000000000000000000000000000
--- a/man/cleanmeteo.Rd
+++ /dev/null
@@ -1,17 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/cleanmeteo.R
-\name{cleanmeteo}
-\alias{cleanmeteo}
-\title{Clean meteorological data}
-\usage{
-cleanmeteo(meteo)
-}
-\arguments{
-\item{meteo}{}
-}
-\value{
-meteo
-}
-\description{
-Clean meteorological data
-}
diff --git a/man/create_ldata.Rd b/man/create_ldata.Rd
deleted file mode 100644
index b99cd1149a02d6634d4fe811922320ae83e9a480..0000000000000000000000000000000000000000
--- a/man/create_ldata.Rd
+++ /dev/null
@@ -1,29 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/create_ldata.R
-\name{create_ldata}
-\alias{create_ldata}
-\title{Write local data matrix}
-\usage{
-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{gdata}{list}
-
-\item{time_max}{numerical value}
-
-\item{nYR_init}{numerical value}
-
-\item{nodeID}{string}
-}
-\value{
-matrix
-}
-\description{
-Function to define local data
-}
-\keyword{METEO}
-\keyword{ldata}
diff --git a/man/diapause.Rd b/man/diapause.Rd
deleted file mode 100644
index d75899d82f72725fa2723b73cc2e1e67cc54fbcd..0000000000000000000000000000000000000000
--- a/man/diapause.Rd
+++ /dev/null
@@ -1,24 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/diapause.R
-\name{diapause}
-\alias{diapause}
-\title{Diapause periods}
-\usage{
-diapause(startFav, endFav, time_max, nYR_init)
-}
-\arguments{
-\item{startFav}{Date. First day of the favourable period.}
-
-\item{endFav}{Date. Last day of the favourable period.}
-
-\item{time_max}{Numerical value. Simulation length (days).}
-
-\item{nYR_init}{Numerical value. Number of years for initialization of the dynamics (default is 2).}
-}
-\value{
-Numeric vector with days of start and end of the favourable period
-}
-\description{
-Function to calculates days for diapause over all the simulated period
-}
-\keyword{diapause}
diff --git a/man/runArboRisk.Rd b/man/runArboRisk.Rd
index 2b2973639dd76fd52103ba9a5ceac1f20685fd36..14b9780cdad59abe8767918e97464e92cb48bfb3 100644
--- a/man/runArboRisk.Rd
+++ b/man/runArboRisk.Rd
@@ -21,11 +21,22 @@ runArboRisk(PARCELLE, METEO, gdata, mMov = NULL, start_date, end_date, nYR_init
\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{nYR_init}{numeric. Number of years used to initialize the population dynamics (default is 1)}
\item{n_sim}{integer}
\item{introduction_pts}{data.frame or data.table}
+
+\item{u0}{matrix Initial population state (patches as columns and in rows: Sh, Eh, Ih, Rh, A1gmI, A1omI, A2hmI, A2gmI, A2omI, Neggs, ninfh, ninfm1, ninfm2)}
+
+\item{v0}{matrix Initial population and time-dependant parameters state (patches as rows and in columns:
+z: diapause (0 = dipause, 1 = favorable period); temperature; kL and kP: carrying capacities for larvae and pupae (rainfall dependent);
+Em: number of eggs ; Lm: number of larvae ; Pm: number of pupae ; Aemm: number of emerging adults ;
+A1hm: number of nulliparous adults seeking for host ; A1gm: number of gorged nulliparous adults ; A1om: number of nulliparous adults seeking for oviposition sites ;
+A2hm: number of parous adults seeking for host ; A2gm: number of gorged parous adults ; A1om: number of parous adults seeking for oviposition sites ;
+prevEggs, nIm1, nIm2, R0, pIh and pIm are continous variables calculated over time)}
+
+\item{initMosq}{numeric. Number of eggs in each patch at t0 (default is 100000).}
}
\value{
data.table
diff --git a/tests/testthat/test-runArboRisk.R b/tests/testthat/test-runArboRisk.R
index 051dec441344a460bb35bfb8a5cd8884753ed827..8e96ef6c802cd655ec2ba1916c4368781f078e77 100644
--- a/tests/testthat/test-runArboRisk.R
+++ b/tests/testthat/test-runArboRisk.R
@@ -1,5 +1,6 @@
test_that("simulation works", {
+ PARCELLE %<>% .[startsWith(PARCELLE$ID, "34"),]
traj <- runArboRisk(PARCELLE = PARCELLE,
METEO = METEO,
@@ -10,3 +11,66 @@ test_that("simulation works", {
expect_length(traj, 2)
})
+
+test_that("simulation works with introduction", {
+
+ PARCELLE %<>% .[startsWith(PARCELLE$ID, "34"),]
+
+ introduction_pts <- build_E_random(
+ period_start <- "2022/08/10" %>% as.Date,
+ period_end <- "2022/08/30" %>% as.Date,
+ n_ind = 1:10,
+ n_events = 5,
+ stage = "Eh",
+ loc = PARCELLE$ID[1])
+
+ traj <- runArboRisk(PARCELLE = PARCELLE,
+ METEO = METEO,
+ n_sim = 2,
+ introduction_pts = introduction_pts)
+
+ expect_type(traj, "list")
+ expect_s3_class(traj[[1]], "data.frame")
+ expect_length(traj, 2)
+
+})
+
+
+test_that("simulation works with introduction and human mobility", {
+
+ f <- system.file("data/SpatVec.shp", package="ArboRisk")
+ SpatVec <- terra::vect(f)
+
+PARCELLE %<>% .[startsWith(PARCELLE$ID, "34"),]
+SpatVec <- SpatVec[SpatVec$ID %in% PARCELLE$ID, ]
+
+ introduction_pts <- build_E_random(
+ period_start <- "2022/08/10" %>% as.Date,
+ period_end <- "2022/08/30" %>% as.Date,
+ n_ind = 1:10,
+ n_events = 5,
+ stage = "Eh",
+ loc = PARCELLE$ID[1])
+
+mMov <- build_mMov(SpatVec,
+ law = "GravExp",
+ param = 0.00001,
+ p2move = 0.9,
+ pCom = 0.7)
+
+ldata <- build_ldata(PARCELLE,
+ METEO,
+ mMov = mMov)
+
+traj <- runArboRisk(PARCELLE = PARCELLE,
+ ldata = ldata,
+ n_sim = 2,
+ start_date = min(METEO$DATE) %>% as.Date,
+ end_date = max(METEO$DATE) %>% as.Date,
+ introduction_pts = introduction_pts)
+
+ 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
index 8da7c70283b78d2b7a0bc62db26fe732ca5d4d48..2a90eb418a6b9468c5b6bc8056e815825ee97f38 100644
--- a/tests/testthat/test-test_intro.R
+++ b/tests/testthat/test-test_intro.R
@@ -7,5 +7,5 @@ test_that("introduction works", {
stage = "Eh",
loc = letters[1:5])
expect_equal(dim(introduction_pts), c(5,4))
- expect_equal(names(introduction_pts), c("time", "node", "n", "select"))
+ expect_equal(names(introduction_pts), c("time", "dest", "n", "select"))
})
diff --git a/work_in_progress/main.R b/work_in_progress/main.R
index a22829e4972e236454e1f685141a1e5b5d6daad9..732f74dabf48063811cb239906ca06bb9fc720b2 100644
--- a/work_in_progress/main.R
+++ b/work_in_progress/main.R
@@ -28,9 +28,10 @@ library(magrittr)
####################
# 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
+# METEO: table with four columns nodeID, date, daily rainfall, daily mean temperature
data(PARCELLE)
+data(METEO)
PARCELLE %<>% .[startsWith(PARCELLE$ID, "11"),]
# PARCELLE %<>% .[1:100,]
METEO
@@ -40,28 +41,30 @@ METEO
###########################################
introduction_pts <- build_E_random(
- period_start <- "2020/03/10" %>% as.Date,
- period_end <- "2020/09/30" %>% as.Date,
+ period_start <- "2022/08/10" %>% as.Date,
+ period_end <- "2022/08/30" %>% as.Date,
n_ind = 1:10,
n_events = 5,
stage = "Eh",
- loc = PARCELLE$ID)
-
+ loc = PARCELLE$ID[1])
######################################
### Introduction of human mobility ###
######################################
-mMov <- build_mMov(PARCELLE,
- group = "NOM_COM",
- within_clust_lev = 1,
- between_clust_lev = 0.1,
- clust_ratio_inout = 0.8,
- verbose = T)
+SpatVec <- "data/SpatVec.shp" %>% terra::vect(.)
+SpatVec <- SpatVec[SpatVec$ID %in% PARCELLE$ID, ]
+
+mMov <- build_mMov(SpatVec,
+ law = "GravExp",
+ param = 0.00001,
+ p2move = 0.2,
+ pCom = 0.7)
ldata <- build_ldata(PARCELLE,
METEO,
- vector_species = "Ae. albopictus")
+ vector_species = "Ae. albopictus",
+ mMov = mMov)
###################
### Simulations ###
@@ -97,7 +100,7 @@ start_time <- Sys.time()
traj <- runArboRisk(PARCELLE = PARCELLE,
vector_species = "Ae. albopictus",
ldata = ldata,
- n_sim = 30,
+ n_sim = 2,
start_date = min(METEO$DATE) %>% as.Date,
end_date = max(METEO$DATE) %>% as.Date,
introduction_pts = introduction_pts)
@@ -105,63 +108,5 @@ traj <- runArboRisk(PARCELLE = PARCELLE,
end_time <- Sys.time()
end_time - start_time
-#####
-##### Explore output
-#####
-
-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)
-PARCELLES$R0 <- traj[, mean(R0), by = node][, V1]
-plot(PARCELLES, "R0", col=rev(heat.colors(12)))
-
-
-plot(result, y = c("Em","Lm", "Pm"))
-plot(result, y = c("A1hm", "A1gm", "A1om"))
-plot(result, y = c("A2hm", "A2gm", "A2om"))
-
-plot(result, y = c("Em","Lm", "Pm"))
-plot(result, y = c("Em","Lm", "Pm","Aemm"))
-plot(result, y = c("Lm", "Pm","Aemm"))
-plot(result, y = c("Em","Lm", "Pm","Aemm","A1hm","A1gm","A1om","A2hm","A2gm","A2om"))
-plot(result, y = c("A1hm","A1gm","A1om","A2hm","A2gm","A2om"))
-plot(result, y = c("A1gmI","A1omI", "A2hmI","A2gmI","A2omI"), index = 57)
-plot(result, y = c("Em","Lm", "Pm"))
-plot(result, y = c("Sh","Eh", "Ih", "Rh"), index = 57)
-plot(result, y = c("Eh", "Ih", "Rh"))
-plot(result, y = c("temperature"))
-
-plot(result, y = c("z"))
-
-traj <- trajectory(result)
-traj %>% setDT
-traj
-
-
-final_data <- melt(traj, id.vars = c("node", "time"),
- measure.vars = c("Em","Lm", "Pm","Aemm","A1hm","A1gm","A1om","A2hm","A2gm","A2om"))
-
-final_data <- melt(traj,
- id.vars = c("node", "time"),
- measure.vars = c("A1gmI", "A1omI","A2hmI", "A2gmI", "A2omI"))
-
-final_data <- melt(traj, id.vars = c("node", "time"),
- measure.vars = c("Eh", "Ih"))
-
-ggplot() +
- geom_line(data = final_data[value > 0],
- aes(x = time, y = value, color = variable, group = variable)) +
- ylim(0, 2e+07) + # max(final_data$value)) + #
- # xlim(100,max(tspan)) +
- annotate("rect", xmin = 67, xmax = 271, ymin = 0, ymax = 2e+07, alpha = .2, fill = "green") +
- annotate("rect", xmin = 432, xmax = 636, ymin = 0, ymax = 2e+07, alpha = .2, fill = "green") +
- annotate("rect", xmin = 797, xmax = 1001, ymin = 0, ymax = 2e+07, alpha = .2, fill = "green")
-
+# in each simulation see if the infection spread other IRIS
+lapply(traj, function(x) x[ninfh>1 & ID != introduction_pts$dest %>% unique, ID] %>% unique %>% length) %>% unlist
diff --git a/work_in_progress/test_TDLM.R b/work_in_progress/test_TDLM.R
new file mode 100644
index 0000000000000000000000000000000000000000..3b1c366556d3a8d44719239c83a21eb228487e8b
--- /dev/null
+++ b/work_in_progress/test_TDLM.R
@@ -0,0 +1,218 @@
+### Date of change: 27/02/2023
+### Authors: Pachka, Ewy, Elena
+
+## This script is under development at the UMR Astre - Cirad.
+
+## The aim of the script is to develop a model to simulate
+## viral dynamics of dengue serotypes in a disease-free context
+## taking into account spatial variability, vector population dynamics
+## and human commuting.
+
+rm(list = ls(all.names = TRUE)) # clear all objects includes hidden objects.
+gc() # free up memory and report the memory usage.
+
+###################################
+### Load packages and functions ###
+###################################
+# devtools::install(quick = T)
+
+library(data.table)
+library(magrittr)
+library(ggplot2)
+
+library(terra)
+library(geodata)
+
+library(TDLM)
+library(sf)
+####################
+### General data ###
+####################
+# input data:
+# PARCELLE : table where each row is a node and columns are attributes
+
+### TEST ARBOCARTO
+filename <- "work_in_progress/data/IRIS-GE/IRIS_GE.SHP"
+#### List of administrative unit ####
+PARCELLES_spa <- vect(filename)
+terra::crs(PARCELLES_spa, describe = T)
+
+### For each administrative unit define:
+## Surface ##
+PARCELLES_spa$SURF_HA <- expanse(PARCELLES_spa, unit="km", transform=TRUE)
+
+## Average altitude ##
+# download raster
+elevation <- elevation_30s(country="FRA", path = "work_in_progress/data/")
+# projet raster
+elevation %<>% project(., "epsg:2154")
+
+## Population ##
+# download raster
+pop <- geodata::population(2020, 0.5, path = "work_in_progress/data/")
+# projet raster
+pop %<>% project(., elevation)
+PARCELLES_spa %<>% project(., elevation)
+# compute average
+PARCELLES_spa$POP <- terra::extract(pop, PARCELLES_spa, fun=sum, na.rm=TRUE) %>% .$population_density
+# Absence of missing values?
+PARCELLES_spa$POP %>% is.na %>% sum %>% equals(0)
+
+# Turn density into number of persons
+PARCELLES_spa$POP %<>% multiply_by(PARCELLES_spa$SURF_HA)
+
+## turn to sf
+PARCELLES_spa %<>% sf::st_as_sf(.)
+# Choose dpt 34
+PARCELLES_spa <- subset(PARCELLES_spa, startsWith(CODE_IRIS, "34"))
+## Dimension of the object
+dim(PARCELLES_spa)
+
+PARCELLES_spa[1:10,]
+
+plot(PARCELLES_spa, max.plot = ncol(PARCELLES_spa))
+
+mass <- data.frame(Population = round(PARCELLES_spa$POP),
+ Outcommuters = round(0.7*PARCELLES_spa$POP))
+
+row.names(mass) <- PARCELLES_spa$CODE_IRIS
+
+dim(mass)
+
+mi <- as.numeric(mass[,1])
+names(mi) <- rownames(mass)
+
+mj <- mi
+
+Oi <- as.numeric(mass[,2])
+names(Oi) <- rownames(mass)
+
+spi <- extract_spatial_information(PARCELLES_spa, id = "CODE_IRIS")
+
+
+distance2 <- spi$distance
+
+
+distance2[1:10, 1:10]
+
+check_format_names(vectors = list(mi = mi, mj = mj, Oi = Oi),
+ matrices = list(distance = distance2),
+ check = "format_and_names")
+
+res <- run_law_model(law = "NGravExp",
+ mass_origin = mi,
+ mass_destination = mj,
+ distance = distance2,
+ opportunity = NULL,
+ param = 0.01,
+ write_proba = TRUE,
+
+ model = "DCM",
+ nb_trips = sum(Oi),
+ out_trips = Oi,
+ average = FALSE,
+ nbrep = 3)
+
+res$proba[1:10,1:10]
+
+res <- run_law(
+ law = "GravExp",
+ mass_origin = mi,
+ distance = distance2,
+ param = 0.01
+)
+
+res$proba
+
+res <- run_law(
+ law = "GravExp",
+ mass_origin = mi,
+ distance = distance2,
+ param = c(0.01,0.02)
+)
+
+res$parameter_1 %>% str
+res$parameter_2
+
+
+res <- run_law(
+ law = "NGravExp",
+ mass_origin = mi,
+ distance = distance2,
+ param = 0.01
+)
+res$proba %>% dim
+res$proba[1:10, 1:10]
+
+res$parameter_1
+res$parameter_2
+
+
+res <- run_law(
+ law = "GravPow",
+ mass_origin = mi,
+ distance = distance2,
+ param = c(0.01,0.02)
+)
+
+res$parameter_1
+res$parameter_2
+
+res <- run_law(
+ law = "NGravPow",
+ mass_origin = mi,
+ distance = distance2,
+ param = c(0.01,0.02)
+)
+
+res$parameter_1
+res$parameter_2
+
+#####
+opportunity <- mass[, 1]
+names(opportunity) <- rownames(mass)
+
+sij <- extract_opportunities(
+ opportunity = opportunity,
+ distance = distance2,
+ check_names = TRUE
+)
+
+
+res <- run_law(
+ law = "Schneider",
+ mass_origin = mi,
+ opportunity = sij,
+ param = c(0.01,0.02)
+)
+
+res$parameter_1
+res$parameter_2
+
+
+
+res <- run_law(
+ law = "Rad",
+ mass_origin = mi,
+ opportunity = sij
+)
+
+res$proba
+
+
+res <- run_law(
+ law = "RadExt",
+ mass_origin = mi,
+ opportunity = sij,
+ param = c(0.01,0.02)
+)
+
+res$parameter_1
+res$parameter_2
+
+proba <- res$proba
+
+for(i in proba %>% nrow %>% seq){
+ proba[i,] <- (0.3 * proba[i,])/sum(proba[i,])
+ proba[i,i] <- 0.7
+ }