diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 818415a3..5844e49b 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -60,14 +60,16 @@ jobs: libudunits2-dev \ libgdal-dev \ libgeos-dev \ - libproj-dev + libproj-dev \ + libharfbuzz-dev \ + libfribidi-dev - name: Install system dependencies if: runner.os == 'macOS' run: | # install spatial dependencies rm '/usr/local/bin/gfortran' - brew install pkg-config gdal proj geos + brew install pkg-config gdal proj geos harfbuzz fribidi - name: Install renv shell: Rscript {0} diff --git a/CHANGELOG.md b/CHANGELOG.md index 127b572c..30fbef2f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,11 @@ this repository. - `quantity_allocation_disagreement` now takes in the use_distance parameter which is FALSE by default. This allows the model to commute the minimum total distance between observed and simulated infestations (@ChrisJones, #130). + +- `validate`, `calibrate`, `pops_multirun`, `auto_manage` and `pops` now take in + network_min_distance, network_max_distance, and network_filename parameters these are + used when the anthropogenic_kernel_type = "network". This allows directed spread along + a network such as a railroad (@ChrisJones and @wenzeslaus, #131) ### Changed diff --git a/DESCRIPTION b/DESCRIPTION index b6ef7495..a8006db0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -58,5 +58,5 @@ Suggests: pkgdown LinkingTo: Rcpp -RoxygenNote: 7.1.1 +RoxygenNote: 7.1.2 VignetteBuilder: knitr diff --git a/R/RcppExports.R b/R/RcppExports.R index 3be02bc6..e8d646fa 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,8 +1,8 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -pops_model_cpp <- function(random_seed, lethal_temperature, lethal_temperature_month, infected, total_exposed, exposed, susceptible, total_populations, total_hosts, mortality_tracker, mortality, quarantine_areas, treatment_maps, treatment_dates, pesticide_duration, resistant, movements, movements_dates, temperature, weather_coefficient, res, rows_cols, reproductive_rate, spatial_indices, season_month_start_end, frequency_config, bool_config, mortality_rate = 0.0, mortality_time_lag = 2L, start_date = "2018-01-01", end_date = "2018-12-31", treatment_method = "ratio", natural_kernel_type = "cauchy", anthropogenic_kernel_type = "cauchy", percent_natural_dispersal = 0.0, natural_distance_scale = 21, anthropogenic_distance_scale = 0.0, natural_dir = "NONE", natural_kappa = 0, anthropogenic_dir = "NONE", anthropogenic_kappa = 0, output_frequency_n = 1L, quarantine_frequency_n = 1L, spreadrate_frequency_n = 1L, mortality_frequency_n = 1L, model_type_ = "SI", latency_period = 0L, establishment_probability = 0, dispersal_percentage = 0.99, overpopulation_config = NULL) { - .Call(`_PoPS_pops_model_cpp`, random_seed, lethal_temperature, lethal_temperature_month, infected, total_exposed, exposed, susceptible, total_populations, total_hosts, mortality_tracker, mortality, quarantine_areas, treatment_maps, treatment_dates, pesticide_duration, resistant, movements, movements_dates, temperature, weather_coefficient, res, rows_cols, reproductive_rate, spatial_indices, season_month_start_end, frequency_config, bool_config, mortality_rate, mortality_time_lag, start_date, end_date, treatment_method, natural_kernel_type, anthropogenic_kernel_type, percent_natural_dispersal, natural_distance_scale, anthropogenic_distance_scale, natural_dir, natural_kappa, anthropogenic_dir, anthropogenic_kappa, output_frequency_n, quarantine_frequency_n, spreadrate_frequency_n, mortality_frequency_n, model_type_, latency_period, establishment_probability, dispersal_percentage, overpopulation_config) +pops_model_cpp <- function(random_seed, lethal_temperature, lethal_temperature_month, infected, total_exposed, exposed, susceptible, total_populations, total_hosts, mortality_tracker, mortality, quarantine_areas, treatment_maps, treatment_dates, pesticide_duration, resistant, movements, movements_dates, temperature, weather_coefficient, bbox, res, rows_cols, reproductive_rate, spatial_indices, season_month_start_end, frequency_config, bool_config, mortality_rate = 0.0, mortality_time_lag = 2L, start_date = "2018-01-01", end_date = "2018-12-31", treatment_method = "ratio", natural_kernel_type = "cauchy", anthropogenic_kernel_type = "cauchy", percent_natural_dispersal = 0.0, natural_distance_scale = 21, anthropogenic_distance_scale = 0.0, natural_dir = "NONE", natural_kappa = 0, anthropogenic_dir = "NONE", anthropogenic_kappa = 0, frequencies_n_config = NULL, model_type_ = "SI", latency_period = 0L, establishment_probability = 0, dispersal_percentage = 0.99, overpopulation_config = NULL, network_config = NULL, network_data_config = NULL) { + .Call(`_PoPS_pops_model_cpp`, random_seed, lethal_temperature, lethal_temperature_month, infected, total_exposed, exposed, susceptible, total_populations, total_hosts, mortality_tracker, mortality, quarantine_areas, treatment_maps, treatment_dates, pesticide_duration, resistant, movements, movements_dates, temperature, weather_coefficient, bbox, res, rows_cols, reproductive_rate, spatial_indices, season_month_start_end, frequency_config, bool_config, mortality_rate, mortality_time_lag, start_date, end_date, treatment_method, natural_kernel_type, anthropogenic_kernel_type, percent_natural_dispersal, natural_distance_scale, anthropogenic_distance_scale, natural_dir, natural_kappa, anthropogenic_dir, anthropogenic_kappa, frequencies_n_config, model_type_, latency_period, establishment_probability, dispersal_percentage, overpopulation_config, network_config, network_data_config) } # Register entry points for exported C++ functions diff --git a/R/auto_manage.R b/R/auto_manage.R index ca569d24..b9441c29 100644 --- a/R/auto_manage.R +++ b/R/auto_manage.R @@ -119,7 +119,10 @@ auto_manage <- function(infected_files, exposed_file = "", mask = NULL, write_outputs = "None", - output_folder_path = "") { + output_folder_path = "", + network_min_distance = 0, + network_max_distance = 0, + network_filename = "") { config <- c() config$random_seed <- random_seed @@ -186,6 +189,9 @@ auto_manage <- function(infected_files, config$mask <- mask config$write_outputs <- write_outputs config$output_folder_path <- output_folder_path + config$network_min_distance <- network_min_distance + config$network_max_distance <- network_max_distance + config$network_filename <- network_filename config <- configuration(config) @@ -321,7 +327,11 @@ auto_manage <- function(infected_files, use_overpopulation_movements = config$use_overpopulation_movements, overpopulation_percentage = overpopulation_percentage, leaving_percentage = leaving_percentage, - leaving_scale_coefficient = leaving_scale_coefficient + leaving_scale_coefficient = leaving_scale_coefficient, + bbox = config$bounding_box, + network_min_distance = config$network_min_distance, + network_max_distance = config$network_max_distance, + network_filename = config$network_filename ) infected_runs <- raster::stack(lapply(1:length(data$infected), function(x) host)) diff --git a/R/calibrate.R b/R/calibrate.R index e5cf3cee..6d38708d 100644 --- a/R/calibrate.R +++ b/R/calibrate.R @@ -175,6 +175,9 @@ calibrate <- function(infected_years_file, verbose = TRUE, write_outputs = "None", output_folder_path = "", + network_min_distance = 0, + network_max_distance = 0, + network_filename = "", use_distance = FALSE, use_rmse = FALSE) { @@ -252,6 +255,9 @@ calibrate <- function(infected_years_file, config$output_folder_path <- output_folder_path config$mortality_frequency <- mortality_frequency config$mortality_frequency_n <- mortality_frequency_n + config$network_min_distance <- network_min_distance + config$network_max_distance <- network_max_distance + config$network_filename <- network_filename config$use_distance <- use_distance config$use_rmse <- use_rmse @@ -340,7 +346,11 @@ calibrate <- function(infected_years_file, use_overpopulation_movements = config$use_overpopulation_movements, overpopulation_percentage = config$overpopulation_percentage, leaving_percentage = config$leaving_percentage, - leaving_scale_coefficient = config$leaving_scale_coefficient + leaving_scale_coefficient = config$leaving_scale_coefficient, + bbox = config$bounding_box, + network_min_distance = config$network_min_distance, + network_max_distance = config$network_max_distance, + network_filename = config$network_filename ) return(data) } diff --git a/R/configuration.R b/R/configuration.R index 93050c3e..6dd9510f 100644 --- a/R/configuration.R +++ b/R/configuration.R @@ -172,6 +172,7 @@ configuration <- function(config) { if (terra::nlyr(infected) > 1) { infected <- output_from_raster_mean_and_sd(infected) } + infected <- terra::classify(infected, matrix(c(NA, 0), ncol = 2, byrow = TRUE), right = NA) } else { config$failure <- infected_check$failed_check return(config) @@ -487,7 +488,9 @@ configuration <- function(config) { # "Log normal", # "Log Normal", "logistic", - "Logistic" + "Logistic", + "network", + "Network" ) if (config$natural_kernel_type %notin% kernel_list) { @@ -668,6 +671,12 @@ configuration <- function(config) { config$xmin <- terra::xmin(config$host) config$ymax <- terra::ymax(config$host) config$ymin <- terra::ymin(config$host) + bounding_box <- c() + bounding_box$north <- config$ymax + bounding_box$south <- config$ymin + bounding_box$west <- config$xmin + bounding_box$east <- config$xmax + config$bounding_box <- bounding_box return(config) } diff --git a/R/helpers.R b/R/helpers.R index 7875ee67..05bdb5cd 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -107,8 +107,8 @@ natural_kernel_error <- "Natural kernel type not one of 'cauchy', 'exponential', 'uniform', 'deterministic neighbor', 'power law', 'hyperbolic secant', 'gamma', 'weibull', 'logistic'" anthropogenic_kernel_error <- - "Anthropogenic kernel type not one of 'cauchy', 'exponential', 'uniform', - 'deterministic neighbor', 'power law', 'hyperbolic secant', 'gamma', 'weibull', 'logistic'" + "Anthropogenic kernel type not one of 'cauchy', 'exponential', 'uniform', 'deterministic + neighbor', 'power law', 'hyperbolic secant', 'gamma', 'weibull', 'logistic', 'network'" covariance_mat_error <- "parameter covariance matrix is not 6 x 6" paramter_means_error <- diff --git a/R/pops.r b/R/pops.r index 2f1191b3..e902828d 100644 --- a/R/pops.r +++ b/R/pops.r @@ -164,6 +164,9 @@ #' true negatives from comparisons (e.g. mask out lakes and oceans from statics #' if modeling terrestrial species). This can also be used to mask out areas #' that can't be managed in the auto_manage function. +#' @param network_min_distance minimum time a propagule rides on the network +#' @param network_max_distance maximum time a propagule rides on the network +#' @param network_filename entire file path for the network file #' #' @useDynLib PoPS, .registration = TRUE #' @importFrom terra app rast xres yres classify extract ext as.points ncol nrow @@ -233,7 +236,10 @@ pops <- function(infected_file, leaving_percentage = 0, leaving_scale_coefficient = 1, exposed_file = "", - mask = NULL) { + mask = NULL, + network_min_distance = 0, + network_max_distance = 0, + network_filename = "") { config <- c() config$random_seed <- random_seed @@ -303,6 +309,10 @@ pops <- function(infected_file, config$mortality_frequency <- mortality_frequency config$mortality_frequency_n <- mortality_frequency_n + config$network_min_distance <- network_min_distance + config$network_max_distance <- network_max_distance + config$network_filename <- network_filename + config <- configuration(config) if (!is.null(config$failure)) { @@ -345,15 +355,11 @@ pops <- function(infected_file, end_date = config$end_date, treatment_method = config$treatment_method, natural_kernel_type = config$natural_kernel_type, - anthropogenic_kernel_type = - config$anthropogenic_kernel_type, - use_anthropogenic_kernel = - config$use_anthropogenic_kernel, - percent_natural_dispersal = - config$percent_natural_dispersal[1], + anthropogenic_kernel_type = config$anthropogenic_kernel_type, + use_anthropogenic_kernel = config$use_anthropogenic_kernel, + percent_natural_dispersal = config$percent_natural_dispersal[1], natural_distance_scale = config$natural_distance_scale[1], - anthropogenic_distance_scale = - config$anthropogenic_distance_scale[1], + anthropogenic_distance_scale = config$anthropogenic_distance_scale[1], natural_dir = config$natural_dir, natural_kappa = config$natural_kappa[1], anthropogenic_dir = config$anthropogenic_dir, @@ -370,22 +376,20 @@ pops <- function(infected_file, use_spreadrates = config$use_spreadrates, model_type_ = config$model_type, latency_period = config$latency_period, - generate_stochasticity = - config$generate_stochasticity, - establishment_stochasticity = - config$establishment_stochasticity, + generate_stochasticity = config$generate_stochasticity, + establishment_stochasticity = config$establishment_stochasticity, movement_stochasticity = config$movement_stochasticity, deterministic = config$deterministic, - establishment_probability = - config$establishment_probability, + establishment_probability = config$establishment_probability, dispersal_percentage = config$dispersal_percentage, - use_overpopulation_movements = - config$use_overpopulation_movements, - overpopulation_percentage = - config$overpopulation_percentage, + use_overpopulation_movements = config$use_overpopulation_movements, + overpopulation_percentage = config$overpopulation_percentage, leaving_percentage = config$leaving_percentage, - leaving_scale_coefficient = - config$leaving_scale_coefficient + leaving_scale_coefficient = config$leaving_scale_coefficient, + bbox = config$bounding_box, + network_min_distance = config$network_min_distance, + network_max_distance = config$network_max_distance, + network_filename = config$network_filename ) return(data) diff --git a/R/pops_model.R b/R/pops_model.R index 5e9a0b31..ef103961 100644 --- a/R/pops_model.R +++ b/R/pops_model.R @@ -65,6 +65,8 @@ #' 'day' or 'time step') in which to calculate and export spread rate #' statistics. #' @param spatial_indices list of all spatial locations with suitable hosts +#' @param bbox bounding box for network kernel +#' #' @return list of vector matrices of infected and susceptible hosts per #' simulated year and associated statistics (e.g. spread rate) #' @export @@ -137,7 +139,11 @@ pops_model <- use_overpopulation_movements = FALSE, overpopulation_percentage = 0.0, leaving_percentage = 0.0, - leaving_scale_coefficient = 1.0) { + leaving_scale_coefficient = 1.0, + bbox = NULL, + network_min_distance = 0, + network_max_distance = 0, + network_filename = "") { # List of overpopulation parameters of type double overpopulation_config <- c() @@ -145,6 +151,26 @@ pops_model <- overpopulation_config$leaving_percentage <- leaving_percentage overpopulation_config$leaving_scale_coefficient <- leaving_scale_coefficient + + # List of frequency n parameters + frequencies_n_config <- c() + frequencies_n_config$output_frequency_n <- output_frequency_n + frequencies_n_config$quarantine_frequency_n <- quarantine_frequency_n + frequencies_n_config$spreadrate_frequency_n <- spreadrate_frequency_n + frequencies_n_config$mortality_frequency_n <- mortality_frequency_n + + # Network configuration + network_config <- NULL; + network_data_config <- NULL; + if (!(is.na(network_filename) || is.null(network_filename) || network_filename == '')) { + network_config <- c() + network_config$network_min_distance <- network_min_distance + network_config$network_max_distance <- network_max_distance + + network_data_config <- c() + network_data_config$network_filename <- network_filename + } + # List of frequencies type string frequency_config <- c() frequency_config$time_step <- time_step @@ -168,6 +194,7 @@ pops_model <- bool_config$deterministic <- deterministic bool_config$use_overpopulation_movements <- use_overpopulation_movements + data <- pops_model_cpp(random_seed = random_seed, lethal_temperature = lethal_temperature, @@ -189,6 +216,7 @@ pops_model <- movements_dates = movements_dates, temperature = temperature, weather_coefficient = weather_coefficient, + bbox = bbox, res = res, rows_cols = rows_cols, reproductive_rate = reproductive_rate, @@ -210,14 +238,14 @@ pops_model <- natural_kappa = natural_kappa, anthropogenic_dir = anthropogenic_dir, anthropogenic_kappa = anthropogenic_kappa, - output_frequency_n = output_frequency_n, - quarantine_frequency_n = quarantine_frequency_n, - spreadrate_frequency_n = spreadrate_frequency_n, - mortality_frequency_n = mortality_frequency_n, + frequencies_n_config = frequencies_n_config, model_type_ = model_type_, latency_period = latency_period, dispersal_percentage = dispersal_percentage, - overpopulation_config = overpopulation_config + establishment_probability = establishment_probability, + overpopulation_config = overpopulation_config, + network_config = network_config, + network_data_config = network_data_config ) } diff --git a/R/pops_multirun.R b/R/pops_multirun.R index aecca215..d8aa1014 100644 --- a/R/pops_multirun.R +++ b/R/pops_multirun.R @@ -91,7 +91,10 @@ pops_multirun <- function(infected_file, exposed_file = "", mask = NULL, write_outputs = "None", - output_folder_path = "") { + output_folder_path = "", + network_min_distance = 0, + network_max_distance = 0, + network_filename = "") { config <- c() config$random_seed <- random_seed config$infected_file <- infected_file @@ -148,7 +151,7 @@ pops_multirun <- function(infected_file, config$number_of_iterations <- number_of_iterations config$number_of_cores <- number_of_cores # add function name for use in configuration function to skip - # function specific specifc configurations namely for validation and + # function specific specific configurations namely for validation and # calibration. config$function_name <- "multirun" config$failure <- NULL @@ -158,6 +161,9 @@ pops_multirun <- function(infected_file, config$output_folder_path <- output_folder_path config$mortality_frequency <- mortality_frequency config$mortality_frequency_n <- mortality_frequency_n + config$network_min_distance <- network_min_distance + config$network_max_distance <- network_max_distance + config$network_filename <- network_filename config <- configuration(config) @@ -182,8 +188,7 @@ pops_multirun <- function(infected_file, random_seed = config$random_seed, use_lethal_temperature = config$use_lethal_temperature, lethal_temperature = config$lethal_temperature, - lethal_temperature_month = - config$lethal_temperature_month, + lethal_temperature_month = config$lethal_temperature_month, infected = config$infected, total_exposed = config$total_exposed, exposed = config$exposed, @@ -216,16 +221,11 @@ pops_multirun <- function(infected_file, end_date = config$end_date, treatment_method = config$treatment_method, natural_kernel_type = config$natural_kernel_type, - anthropogenic_kernel_type = - config$anthropogenic_kernel_type, - use_anthropogenic_kernel = - config$use_anthropogenic_kernel, - percent_natural_dispersal = - config$percent_natural_dispersal[i], - natural_distance_scale = - config$natural_distance_scale[i], - anthropogenic_distance_scale = - config$anthropogenic_distance_scale[i], + anthropogenic_kernel_type = config$anthropogenic_kernel_type, + use_anthropogenic_kernel = config$use_anthropogenic_kernel, + percent_natural_dispersal = config$percent_natural_dispersal[i], + natural_distance_scale = config$natural_distance_scale[i], + anthropogenic_distance_scale = config$anthropogenic_distance_scale[i], natural_dir = config$natural_dir, natural_kappa = config$natural_kappa[i], anthropogenic_dir = config$anthropogenic_dir, @@ -242,19 +242,20 @@ pops_multirun <- function(infected_file, use_spreadrates = config$use_spreadrates, model_type_ = config$model_type, latency_period = config$latency_period, - generate_stochasticity = - config$generate_stochasticity, - establishment_stochasticity = - config$establishment_stochasticity, + generate_stochasticity = config$generate_stochasticity, + establishment_stochasticity = config$establishment_stochasticity, movement_stochasticity = config$movement_stochasticity, deterministic = config$deterministic, - establishment_probability = - config$establishment_probability, + establishment_probability = config$establishment_probability, dispersal_percentage = config$dispersal_percentage, use_overpopulation_movements = config$use_overpopulation_movements, overpopulation_percentage = config$overpopulation_percentage, leaving_percentage = config$leaving_percentage, - leaving_scale_coefficient = config$leaving_scale_coefficient + leaving_scale_coefficient = config$leaving_scale_coefficient, + bbox = config$bounding_box, + network_min_distance = config$network_min_distance, + network_max_distance = config$network_max_distance, + network_filename = config$network_filename ) run <- c() diff --git a/R/validate.R b/R/validate.R index 9c3cbcc5..8149bceb 100644 --- a/R/validate.R +++ b/R/validate.R @@ -106,6 +106,9 @@ validate <- function(infected_years_file, write_outputs = "None", output_folder_path = "", point_file = "", + network_min_distance = 0, + network_max_distance = 0, + network_filename = "", use_distance = FALSE, use_configuration = FALSE) { config <- c() @@ -165,7 +168,7 @@ validate <- function(infected_years_file, config$number_of_iterations <- number_of_iterations config$number_of_cores <- number_of_cores # add function name for use in configuration function to skip - # function specific specifc configurations namely for validation and + # function specific specific configurations namely for validation and # calibration. config$function_name <- "validate" config$failure <- NULL @@ -175,6 +178,9 @@ validate <- function(infected_years_file, config$mortality_frequency <- mortality_frequency config$mortality_frequency_n <- mortality_frequency_n config$point_file <- point_file + config$network_min_distance <- network_min_distance + config$network_max_distance <- network_max_distance + config$network_filename <- network_filename config$use_configuration <- use_configuration config$use_distance <- use_distance @@ -201,8 +207,7 @@ validate <- function(infected_years_file, random_seed = config$random_seed, use_lethal_temperature = config$use_lethal_temperature, lethal_temperature = config$lethal_temperature, - lethal_temperature_month = - config$lethal_temperature_month, + lethal_temperature_month = config$lethal_temperature_month, infected = config$infected, total_exposed = config$total_exposed, exposed = config$exposed, @@ -235,16 +240,11 @@ validate <- function(infected_years_file, end_date = config$end_date, treatment_method = config$treatment_method, natural_kernel_type = config$natural_kernel_type, - anthropogenic_kernel_type = - config$anthropogenic_kernel_type, - use_anthropogenic_kernel = - config$use_anthropogenic_kernel, - percent_natural_dispersal = - config$percent_natural_dispersal[i], - natural_distance_scale = - config$natural_distance_scale[i], - anthropogenic_distance_scale = - config$anthropogenic_distance_scale[i], + anthropogenic_kernel_type = config$anthropogenic_kernel_type, + use_anthropogenic_kernel = config$use_anthropogenic_kernel, + percent_natural_dispersal = config$percent_natural_dispersal[i], + natural_distance_scale = config$natural_distance_scale[i], + anthropogenic_distance_scale = config$anthropogenic_distance_scale[i], natural_dir = config$natural_dir, natural_kappa = config$natural_kappa[i], anthropogenic_dir = config$anthropogenic_dir, @@ -261,19 +261,20 @@ validate <- function(infected_years_file, use_spreadrates = config$use_spreadrates, model_type_ = config$model_type, latency_period = config$latency_period, - generate_stochasticity = - config$generate_stochasticity, - establishment_stochasticity = - config$establishment_stochasticity, + generate_stochasticity = config$generate_stochasticity, + establishment_stochasticity = config$establishment_stochasticity, movement_stochasticity = config$movement_stochasticity, deterministic = config$deterministic, - establishment_probability = - config$establishment_probability, + establishment_probability = config$establishment_probability, dispersal_percentage = config$dispersal_percentage, use_overpopulation_movements = config$use_overpopulation_movements, overpopulation_percentage = config$overpopulation_percentage, leaving_percentage = config$leaving_percentage, - leaving_scale_coefficient = config$leaving_scale_coefficient + leaving_scale_coefficient = config$leaving_scale_coefficient, + bbox = config$bounding_box, + network_min_distance = config$network_min_distance, + network_max_distance = config$network_max_distance, + network_filename = config$network_filename ) diff --git a/inst/cpp/pops-core b/inst/cpp/pops-core index f127a7b1..8c24e4e2 160000 --- a/inst/cpp/pops-core +++ b/inst/cpp/pops-core @@ -1 +1 @@ -Subproject commit f127a7b166e4f7cfb88c066a5665f2e1f66adc2f +Subproject commit 8c24e4e228e73e062049c3750b0a32c3e1b20e82 diff --git a/inst/extdata/simple20x20/nodes.csv b/inst/extdata/simple20x20/nodes.csv new file mode 100644 index 00000000..1f213384 --- /dev/null +++ b/inst/extdata/simple20x20/nodes.csv @@ -0,0 +1,6 @@ +1,50,1950 +2,950,1950 +3,450,1550 +4,50,1050 +5,950,1050 +6,1450,550 diff --git a/inst/extdata/simple20x20/segments.csv b/inst/extdata/simple20x20/segments.csv new file mode 100644 index 00000000..453c1a59 --- /dev/null +++ b/inst/extdata/simple20x20/segments.csv @@ -0,0 +1,5 @@ +1,2,50;1950;100;1950;150;1950;200;1950;250;1950;300;1950;350;1950;400;1950;450;1950;500;1950;550;1950;600;1950;650;1950;700;1950;750;1950;800;1950;850;1950;900;1950;950;1950 +1,4,50;1950;50;1900;50;1850;50;1800;50;1750;50;1700;50;1650;50;1600;50;1550;50;1500;50;1450;50;1400;50;1350;50;1300;50;1250;50;1200;50;1150;50;1100;50;1050 +1,3,50;1950;100;1900;150;1850;200;1800;250;1750;300;1700;350;1650;400;1600;450;1550 +3,5,450;1550;500;1500;550;1450;600;1400;650;1350;700;1300;750;1250;800;1200;850;1150;900;1100;950;1050 +5,6,950;1050;1000;1000;1050;950;1100;900;1150;850;1200;800;1250;750;1300;700;1350;650;1300;600;1450;550 diff --git a/inst/include/PoPS_RcppExports.h b/inst/include/PoPS_RcppExports.h index 5881274d..356314b1 100644 --- a/inst/include/PoPS_RcppExports.h +++ b/inst/include/PoPS_RcppExports.h @@ -24,17 +24,17 @@ namespace PoPS { } } - inline List pops_model_cpp(int random_seed, double lethal_temperature, int lethal_temperature_month, IntegerMatrix infected, IntegerMatrix total_exposed, std::vector exposed, IntegerMatrix susceptible, IntegerMatrix total_populations, IntegerMatrix total_hosts, std::vector mortality_tracker, IntegerMatrix mortality, IntegerMatrix quarantine_areas, std::vector treatment_maps, std::vector treatment_dates, std::vector pesticide_duration, IntegerMatrix resistant, std::vector> movements, std::vector movements_dates, std::vector temperature, std::vector weather_coefficient, List res, List rows_cols, double reproductive_rate, std::vector> spatial_indices, List season_month_start_end, List frequency_config, List bool_config, double mortality_rate = 0.0, int mortality_time_lag = 2, std::string start_date = "2018-01-01", std::string end_date = "2018-12-31", std::string treatment_method = "ratio", std::string natural_kernel_type = "cauchy", std::string anthropogenic_kernel_type = "cauchy", double percent_natural_dispersal = 0.0, double natural_distance_scale = 21, double anthropogenic_distance_scale = 0.0, std::string natural_dir = "NONE", double natural_kappa = 0, std::string anthropogenic_dir = "NONE", double anthropogenic_kappa = 0, int output_frequency_n = 1, int quarantine_frequency_n = 1, int spreadrate_frequency_n = 1, int mortality_frequency_n = 1, std::string model_type_ = "SI", int latency_period = 0, double establishment_probability = 0, double dispersal_percentage = 0.99, Nullable overpopulation_config = R_NilValue) { + inline List pops_model_cpp(int random_seed, double lethal_temperature, int lethal_temperature_month, IntegerMatrix infected, IntegerMatrix total_exposed, std::vector exposed, IntegerMatrix susceptible, IntegerMatrix total_populations, IntegerMatrix total_hosts, std::vector mortality_tracker, IntegerMatrix mortality, IntegerMatrix quarantine_areas, std::vector treatment_maps, std::vector treatment_dates, std::vector pesticide_duration, IntegerMatrix resistant, std::vector> movements, std::vector movements_dates, std::vector temperature, std::vector weather_coefficient, List bbox, List res, List rows_cols, double reproductive_rate, std::vector> spatial_indices, List season_month_start_end, List frequency_config, List bool_config, double mortality_rate = 0.0, int mortality_time_lag = 2, std::string start_date = "2018-01-01", std::string end_date = "2018-12-31", std::string treatment_method = "ratio", std::string natural_kernel_type = "cauchy", std::string anthropogenic_kernel_type = "cauchy", double percent_natural_dispersal = 0.0, double natural_distance_scale = 21, double anthropogenic_distance_scale = 0.0, std::string natural_dir = "NONE", double natural_kappa = 0, std::string anthropogenic_dir = "NONE", double anthropogenic_kappa = 0, Nullable frequencies_n_config = R_NilValue, std::string model_type_ = "SI", int latency_period = 0, double establishment_probability = 0, double dispersal_percentage = 0.99, Nullable overpopulation_config = R_NilValue, Nullable network_config = R_NilValue, Nullable network_data_config = R_NilValue) { typedef SEXP(*Ptr_pops_model_cpp)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP); static Ptr_pops_model_cpp p_pops_model_cpp = NULL; if (p_pops_model_cpp == NULL) { - validateSignature("List(*pops_model_cpp)(int,double,int,IntegerMatrix,IntegerMatrix,std::vector,IntegerMatrix,IntegerMatrix,IntegerMatrix,std::vector,IntegerMatrix,IntegerMatrix,std::vector,std::vector,std::vector,IntegerMatrix,std::vector>,std::vector,std::vector,std::vector,List,List,double,std::vector>,List,List,List,double,int,std::string,std::string,std::string,std::string,std::string,double,double,double,std::string,double,std::string,double,int,int,int,int,std::string,int,double,double,Nullable)"); + validateSignature("List(*pops_model_cpp)(int,double,int,IntegerMatrix,IntegerMatrix,std::vector,IntegerMatrix,IntegerMatrix,IntegerMatrix,std::vector,IntegerMatrix,IntegerMatrix,std::vector,std::vector,std::vector,IntegerMatrix,std::vector>,std::vector,std::vector,std::vector,List,List,List,double,std::vector>,List,List,List,double,int,std::string,std::string,std::string,std::string,std::string,double,double,double,std::string,double,std::string,double,Nullable,std::string,int,double,double,Nullable,Nullable,Nullable)"); p_pops_model_cpp = (Ptr_pops_model_cpp)R_GetCCallable("PoPS", "_PoPS_pops_model_cpp"); } RObject rcpp_result_gen; { RNGScope RCPP_rngScope_gen; - rcpp_result_gen = p_pops_model_cpp(Shield(Rcpp::wrap(random_seed)), Shield(Rcpp::wrap(lethal_temperature)), Shield(Rcpp::wrap(lethal_temperature_month)), Shield(Rcpp::wrap(infected)), Shield(Rcpp::wrap(total_exposed)), Shield(Rcpp::wrap(exposed)), Shield(Rcpp::wrap(susceptible)), Shield(Rcpp::wrap(total_populations)), Shield(Rcpp::wrap(total_hosts)), Shield(Rcpp::wrap(mortality_tracker)), Shield(Rcpp::wrap(mortality)), Shield(Rcpp::wrap(quarantine_areas)), Shield(Rcpp::wrap(treatment_maps)), Shield(Rcpp::wrap(treatment_dates)), Shield(Rcpp::wrap(pesticide_duration)), Shield(Rcpp::wrap(resistant)), Shield(Rcpp::wrap(movements)), Shield(Rcpp::wrap(movements_dates)), Shield(Rcpp::wrap(temperature)), Shield(Rcpp::wrap(weather_coefficient)), Shield(Rcpp::wrap(res)), Shield(Rcpp::wrap(rows_cols)), Shield(Rcpp::wrap(reproductive_rate)), Shield(Rcpp::wrap(spatial_indices)), Shield(Rcpp::wrap(season_month_start_end)), Shield(Rcpp::wrap(frequency_config)), Shield(Rcpp::wrap(bool_config)), Shield(Rcpp::wrap(mortality_rate)), Shield(Rcpp::wrap(mortality_time_lag)), Shield(Rcpp::wrap(start_date)), Shield(Rcpp::wrap(end_date)), Shield(Rcpp::wrap(treatment_method)), Shield(Rcpp::wrap(natural_kernel_type)), Shield(Rcpp::wrap(anthropogenic_kernel_type)), Shield(Rcpp::wrap(percent_natural_dispersal)), Shield(Rcpp::wrap(natural_distance_scale)), Shield(Rcpp::wrap(anthropogenic_distance_scale)), Shield(Rcpp::wrap(natural_dir)), Shield(Rcpp::wrap(natural_kappa)), Shield(Rcpp::wrap(anthropogenic_dir)), Shield(Rcpp::wrap(anthropogenic_kappa)), Shield(Rcpp::wrap(output_frequency_n)), Shield(Rcpp::wrap(quarantine_frequency_n)), Shield(Rcpp::wrap(spreadrate_frequency_n)), Shield(Rcpp::wrap(mortality_frequency_n)), Shield(Rcpp::wrap(model_type_)), Shield(Rcpp::wrap(latency_period)), Shield(Rcpp::wrap(establishment_probability)), Shield(Rcpp::wrap(dispersal_percentage)), Shield(Rcpp::wrap(overpopulation_config))); + rcpp_result_gen = p_pops_model_cpp(Shield(Rcpp::wrap(random_seed)), Shield(Rcpp::wrap(lethal_temperature)), Shield(Rcpp::wrap(lethal_temperature_month)), Shield(Rcpp::wrap(infected)), Shield(Rcpp::wrap(total_exposed)), Shield(Rcpp::wrap(exposed)), Shield(Rcpp::wrap(susceptible)), Shield(Rcpp::wrap(total_populations)), Shield(Rcpp::wrap(total_hosts)), Shield(Rcpp::wrap(mortality_tracker)), Shield(Rcpp::wrap(mortality)), Shield(Rcpp::wrap(quarantine_areas)), Shield(Rcpp::wrap(treatment_maps)), Shield(Rcpp::wrap(treatment_dates)), Shield(Rcpp::wrap(pesticide_duration)), Shield(Rcpp::wrap(resistant)), Shield(Rcpp::wrap(movements)), Shield(Rcpp::wrap(movements_dates)), Shield(Rcpp::wrap(temperature)), Shield(Rcpp::wrap(weather_coefficient)), Shield(Rcpp::wrap(bbox)), Shield(Rcpp::wrap(res)), Shield(Rcpp::wrap(rows_cols)), Shield(Rcpp::wrap(reproductive_rate)), Shield(Rcpp::wrap(spatial_indices)), Shield(Rcpp::wrap(season_month_start_end)), Shield(Rcpp::wrap(frequency_config)), Shield(Rcpp::wrap(bool_config)), Shield(Rcpp::wrap(mortality_rate)), Shield(Rcpp::wrap(mortality_time_lag)), Shield(Rcpp::wrap(start_date)), Shield(Rcpp::wrap(end_date)), Shield(Rcpp::wrap(treatment_method)), Shield(Rcpp::wrap(natural_kernel_type)), Shield(Rcpp::wrap(anthropogenic_kernel_type)), Shield(Rcpp::wrap(percent_natural_dispersal)), Shield(Rcpp::wrap(natural_distance_scale)), Shield(Rcpp::wrap(anthropogenic_distance_scale)), Shield(Rcpp::wrap(natural_dir)), Shield(Rcpp::wrap(natural_kappa)), Shield(Rcpp::wrap(anthropogenic_dir)), Shield(Rcpp::wrap(anthropogenic_kappa)), Shield(Rcpp::wrap(frequencies_n_config)), Shield(Rcpp::wrap(model_type_)), Shield(Rcpp::wrap(latency_period)), Shield(Rcpp::wrap(establishment_probability)), Shield(Rcpp::wrap(dispersal_percentage)), Shield(Rcpp::wrap(overpopulation_config)), Shield(Rcpp::wrap(network_config)), Shield(Rcpp::wrap(network_data_config))); } if (rcpp_result_gen.inherits("interrupted-error")) throw Rcpp::internal::InterruptedException(); diff --git a/inst/include/anthropogenic_kernel.hpp b/inst/include/anthropogenic_kernel.hpp new file mode 100644 index 00000000..1831ced5 --- /dev/null +++ b/inst/include/anthropogenic_kernel.hpp @@ -0,0 +1,94 @@ +/* + * PoPS model - main disperal kernel + * + * Copyright (C) 2019-2021 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_ANTHROPOGENIC_KERNEL_HPP +#define POPS_ANTHROPOGENIC_KERNEL_HPP + +#include "radial_kernel.hpp" +#include "deterministic_kernel.hpp" +#include "uniform_kernel.hpp" +#include "neighbor_kernel.hpp" +#include "network_kernel.hpp" +#include "kernel_types.hpp" +#include "kernel_base.hpp" +#include "config.hpp" + +#include + +namespace pops { + +/** + * @brief Create anthropogenic kernel from configuration + * + * Same structure as the natural kernel, but the parameters are for anthropogenic + * kernel when available. + * + * @param dispersers The disperser raster (reference, for deterministic kernel) + * @param network Network (initialized or not) + * @return Created kernel + */ +template +std::unique_ptr> create_anthro_kernel( + const Config& config, + const IntegerRaster& dispersers, + const Network& network) +{ + auto anthro_kernel = kernel_type_from_string(config.anthro_kernel_type); + if (anthro_kernel == DispersalKernelType::Uniform) { + using Kernel = DynamicWrapperKernel; + return std::unique_ptr(new Kernel(config.rows, config.cols)); + } + else if (anthro_kernel == DispersalKernelType::DeterministicNeighbor) { + using Kernel = + DynamicWrapperKernel; + return std::unique_ptr( + new Kernel(direction_from_string(config.anthro_direction))); + } + else if (anthro_kernel == DispersalKernelType::Network) { + using Kernel = + DynamicWrapperKernel, Generator>; + return std::unique_ptr(new Kernel( + network, config.network_min_distance, config.network_max_distance)); + } + else if (config.deterministic) { + using Kernel = DynamicWrapperKernel< + DeterministicDispersalKernel, + Generator>; + return std::unique_ptr(new Kernel( + anthro_kernel, + dispersers, + config.dispersal_percentage, + config.ew_res, + config.ns_res, + config.anthro_scale, + config.shape)); + } + else { + using Kernel = + DynamicWrapperKernel, Generator>; + return std::unique_ptr(new Kernel( + config.ew_res, + config.ns_res, + anthro_kernel, + config.anthro_scale, + direction_from_string(config.anthro_direction), + config.anthro_kappa, + config.shape)); + } +} + +} // namespace pops + +#endif // POPS_ANTHROPOGENIC_KERNEL_HPP diff --git a/inst/include/config.hpp b/inst/include/config.hpp index 31da92d0..d45e75c7 100644 --- a/inst/include/config.hpp +++ b/inst/include/config.hpp @@ -25,6 +25,7 @@ #define POPS_CONFIG_HPP #include "scheduling.hpp" +#include "utils.hpp" #include @@ -40,6 +41,7 @@ class Config int cols{0}; double ew_res{0}; double ns_res{0}; + BBox bbox; // Reduced stochasticity bool generate_stochasticity{true}; bool establishment_stochasticity{true}; @@ -65,6 +67,8 @@ class Config std::string anthro_kernel_type; double anthro_scale{0}; std::string anthro_direction; + double network_min_distance{0}; + double network_max_distance{0}; double anthro_kappa{0}; double shape{1.0}; // Treatments diff --git a/inst/include/deterministic_kernel.hpp b/inst/include/deterministic_kernel.hpp index 10d30435..e3b2ba24 100644 --- a/inst/include/deterministic_kernel.hpp +++ b/inst/include/deterministic_kernel.hpp @@ -287,6 +287,15 @@ class DeterministicDispersalKernel return std::make_tuple(row + row_movement, col + col_movement); } + /*! \copydoc RadialDispersalKernel::is_cell_eligible() + */ + bool is_cell_eligible(int row, int col) + { + UNUSED(row); + UNUSED(col); + return true; + } + /*! Returns true if the kernel class support a given kernel type * * \warning This function is experimental and may be removed or diff --git a/inst/include/kernel.hpp b/inst/include/kernel.hpp index 173b9a86..67e7c6c9 100644 --- a/inst/include/kernel.hpp +++ b/inst/include/kernel.hpp @@ -1,7 +1,7 @@ /* * PoPS model - main disperal kernel * - * Copyright (C) 2019-2020 by the authors. + * Copyright (C) 2019-2021 by the authors. * * Authors: Vaclav Petras (wenzeslaus gmail com) * @@ -46,8 +46,10 @@ #ifndef POPS_KERNEL_HPP #define POPS_KERNEL_HPP -#include "switch_kernel.hpp" +#include "natural_kernel.hpp" +#include "anthropogenic_kernel.hpp" #include "natural_anthropogenic_kernel.hpp" +#include "config.hpp" namespace pops { @@ -75,10 +77,25 @@ namespace pops { * See NaturalAnthropogenicDispersalKernel and SwitchDispersalKernel for further * documentation. */ -template +template using DispersalKernel = NaturalAnthropogenicDispersalKernel< - SwitchDispersalKernel, - SwitchDispersalKernel>; + KernelInterface, + KernelInterface>; + +template +DispersalKernel create_dynamic_kernel( + const Config& config, + const IntegerRaster& dispersers, + const Network& network) +{ + return DispersalKernel( + create_natural_kernel( + config, dispersers), + create_anthro_kernel( + config, dispersers, network), + config.use_anthropogenic_kernel, + config.percent_natural_dispersal); +} } // namespace pops diff --git a/inst/include/kernel_base.hpp b/inst/include/kernel_base.hpp new file mode 100644 index 00000000..3e60d828 --- /dev/null +++ b/inst/include/kernel_base.hpp @@ -0,0 +1,88 @@ +/* + * PoPS model - disperal kernels + * + * Copyright (C) 2021 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef KERNEL_BASE_HPP +#define KERNEL_BASE_HPP + +#include "kernel_types.hpp" + +namespace pops { + +/** + * Interface which all dynamically used kernels need to inherit from. + */ +template +class KernelInterface +{ +public: + /*! \copydoc RadialDispersalKernel::operator()() + */ + virtual std::tuple operator()(Generator& generator, int row, int col) = 0; + + virtual bool is_cell_eligible(int row, int col) = 0; + + /*! \copydoc RadialDispersalKernel::supports_kernel() + */ + virtual bool supports_kernel(const DispersalKernelType type) = 0; + + virtual ~KernelInterface() = default; +}; + +/** + * A class which turns any kernel into a kernel usable with KernelInterface. + * + * Kernel needs to implement operator() and is_cell_eligible function. + */ +template +class DynamicWrapperKernel : public KernelInterface +{ +public: + /** Creates an interal copy of the provided kernel */ + DynamicWrapperKernel(const ActualKernel& kernel) : kernel_(kernel) {} + + /** Creates a new internal kernel instance from the provided parameters */ + template + DynamicWrapperKernel(Args&&... args) : kernel_(std::forward(args)...) + {} + + /*! \copydoc RadialDispersalKernel::operator()() + */ + std::tuple operator()(Generator& generator, int row, int col) override + { + return kernel_.operator()(generator, row, col); + } + + /*! \copydoc RadialDispersalKernel::is_cell_eligible() + */ + bool is_cell_eligible(int row, int col) override + { + return DynamicWrapperKernel::kernel_.is_cell_eligible( + row, col); + } + + /*! \copydoc RadialDispersalKernel::supports_kernel() + */ + bool supports_kernel(const DispersalKernelType type) override + { + return ActualKernel::supports_kernel(type); + } + +protected: + ActualKernel kernel_; +}; + +} // namespace pops + +#endif // KERNEL_BASE_HPP diff --git a/inst/include/kernel_types.hpp b/inst/include/kernel_types.hpp index 0c0d8a72..727052fb 100644 --- a/inst/include/kernel_types.hpp +++ b/inst/include/kernel_types.hpp @@ -64,6 +64,7 @@ enum class DispersalKernelType Normal, //!< Normal dispersal kernel LogNormal, //!< Log-normal dispersal kernel Logistic, //!< Logistic dispersal kernel + Network, //!< Network spread None, //!< No dispersal kernel (no spread) }; @@ -111,6 +112,8 @@ inline DispersalKernelType kernel_type_from_string(const std::string& text) return DispersalKernelType::LogNormal; else if (text == "logistic" || text == "Logistic") return DispersalKernelType::Logistic; + else if (text == "network" || text == "Network") + return DispersalKernelType::Network; else if (text == "none" || text == "None" || text == "NONE" || text.empty()) return DispersalKernelType::None; else diff --git a/inst/include/model.hpp b/inst/include/model.hpp index 9252e197..1645a5e6 100644 --- a/inst/include/model.hpp +++ b/inst/include/model.hpp @@ -1,7 +1,7 @@ /* * Tests for the PoPS Model class. * - * Copyright (C) 2020 by the authors. + * Copyright (C) 2020-2021 by the authors. * * Authors: Vaclav Petras * @@ -28,6 +28,7 @@ #include "treatments.hpp" #include "spread_rate.hpp" #include "simulation.hpp" +#include "switch_kernel.hpp" #include "kernel.hpp" #include "scheduling.hpp" #include "quarantine.hpp" @@ -36,7 +37,13 @@ namespace pops { -template +template< + typename IntegerRaster, + typename FloatRaster, + typename RasterIndex, + typename Generator = std::default_random_engine, + typename KernelFactory = DispersalKernel( + const Config&, const IntegerRaster&, const Network&)> class Model { protected: @@ -46,46 +53,10 @@ class Model UniformDispersalKernel uniform_kernel; DeterministicNeighborDispersalKernel natural_neighbor_kernel; DeterministicNeighborDispersalKernel anthro_neighbor_kernel; - Simulation simulation_; + Simulation simulation_; + KernelFactory& kernel_factory_; unsigned last_index{0}; - /** - * @brief Create natural kernel - * - * Kernel parameters are taken from the configuration. - * - * @param dispersers The disperser raster (reference, for deterministic kernel) - * @return Created kernel - */ - SwitchDispersalKernel - create_natural_kernel(const IntegerRaster& dispersers) - { - RadialDispersalKernel radial_kernel( - config_.ew_res, - config_.ns_res, - natural_kernel, - config_.natural_scale, - direction_from_string(config_.natural_direction), - config_.natural_kappa, - config_.shape); - DeterministicDispersalKernel deterministic_kernel( - natural_kernel, - dispersers, - config_.dispersal_percentage, - config_.ew_res, - config_.ns_res, - config_.natural_scale, - config_.shape); - SwitchDispersalKernel selectable_kernel( - natural_kernel, - radial_kernel, - deterministic_kernel, - uniform_kernel, - natural_neighbor_kernel, - config_.deterministic); - return selectable_kernel; - } - /** * @brief Create overpopulation movement kernel * @@ -94,10 +65,12 @@ class Model * by the leaving scale coefficient. * * @param dispersers The disperser raster (reference, for deterministic kernel) + * @param network Network (initialized or not) * @return Created kernel */ - SwitchDispersalKernel - create_overpopulation_movement_kernel(const IntegerRaster& dispersers) + SwitchDispersalKernel + create_overpopulation_movement_kernel( + const IntegerRaster& dispersers, const Network& network) { RadialDispersalKernel radial_kernel( config_.ew_res, @@ -115,56 +88,24 @@ class Model config_.ns_res, config_.natural_scale * config_.leaving_scale_coefficient, config_.shape); - SwitchDispersalKernel selectable_kernel( + NetworkDispersalKernel network_kernel( + network, config_.network_min_distance, config_.network_max_distance); + SwitchDispersalKernel selectable_kernel( natural_kernel, radial_kernel, deterministic_kernel, uniform_kernel, + network_kernel, natural_neighbor_kernel, config_.deterministic); return selectable_kernel; } - /** - * @brief Create anthropogenic kernel - * - * Same structure as the natural kernel, but the parameters are for anthropogenic - * kernel when available. - * - * @param dispersers The disperser raster (reference, for deterministic kernel) - * @return Created kernel - */ - SwitchDispersalKernel - create_anthro_kernel(const IntegerRaster& dispersers) - { - RadialDispersalKernel radial_kernel( - config_.ew_res, - config_.ns_res, - anthro_kernel, - config_.anthro_scale, - direction_from_string(config_.anthro_direction), - config_.anthro_kappa, - config_.shape); - DeterministicDispersalKernel deterministic_kernel( - anthro_kernel, - dispersers, - config_.dispersal_percentage, - config_.ew_res, - config_.ns_res, - config_.anthro_scale, - config_.shape); - SwitchDispersalKernel selectable_kernel( - anthro_kernel, - radial_kernel, - deterministic_kernel, - uniform_kernel, - anthro_neighbor_kernel, - config_.deterministic); - return selectable_kernel; - } - public: - Model(const Config& config) + Model( + const Config& config, + KernelFactory& kernel_factory = + create_dynamic_kernel) : config_(config), natural_kernel(kernel_type_from_string(config.natural_kernel_type)), anthro_kernel(kernel_type_from_string(config.anthro_kernel_type)), @@ -179,7 +120,8 @@ class Model config.latency_period_steps, config.generate_stochasticity, config.establishment_stochasticity, - config.movement_stochasticity) + config.movement_stochasticity), + kernel_factory_(kernel_factory) {} /** @@ -207,25 +149,28 @@ class Model * infected + exposed + susceptible in the cell. * @param[in,out] total_populations All host and non-host individuals in the area * @param[out] dispersers Dispersing individuals (used internally) - * @param exposed[in,out] Exposed hosts (if SEI model is active) - * @param mortality_tracker[in,out] Mortality tracker used to generate *died*. + * @param[in,out] total_exposed Sum of all exposed hosts (if SEI model is active) + * @param[in,out] exposed Exposed hosts (if SEI model is active) + * @param[in,out] mortality_tracker Mortality tracker used to generate *died*. * Expectation is that mortality tracker is of length (1/mortality_rate + * mortality_time_lag) - * @param died[out] Infected hosts which died this step based on the mortality + * @param[out] died Infected hosts which died this step based on the mortality * schedule - * @param temperatures[in] Vector of temperatures used to evaluate lethal + * @param[in] temperatures Vector of temperatures used to evaluate lethal * temperature - * @param weather_coefficient[in] Weather coefficient (for the current step) - * @param treatments[in,out] Treatments to be applied (also tracks use of + * @param[in] weather_coefficient Weather coefficient (for the current step) + * @param[in,out] treatments Treatments to be applied (also tracks use of * treatments) - * @param resistant[in,out] Resistant hosts (host temporarily removed from + * @param[in,out] resistant Resistant hosts (host temporarily removed from * susceptible hosts) - * @param outside_dispersers[in,out] Dispersers escaping the rasters (adds to the + * @param[in,out] outside_dispersers Dispersers escaping the rasters (adds to the * vector) - * @param spread_rate[in,out] Spread rate tracker - * @param quarantine[in,out] Quarantine escape tracker - * @param quarantine_areas[in] Quarantine areas - * @param movements[in] Table of host movements + * @param[in,out] spread_rate Spread rate tracker + * @param[in,out] quarantine Quarantine escape tracker + * @param[in] quarantine_areas Quarantine areas + * @param[in] movements Table of host movements + * @param network Network (initialized or Network::null_network() if unused) + * @param[in,out] suitable_cells List of indices of cells with hosts * * @note The parameters roughly correspond to Simulation::disperse() * and Simulation::disperse_and_infect() functions, so these can be used @@ -238,6 +183,7 @@ class Model IntegerRaster& total_populations, IntegerRaster& total_hosts, IntegerRaster& dispersers, + IntegerRaster& established_dispersers, IntegerRaster& total_exposed, std::vector& exposed, std::vector& mortality_tracker, @@ -251,6 +197,7 @@ class Model QuarantineEscape& quarantine, // out const IntegerRaster& quarantine_areas, const std::vector> movements, + const Network& network, std::vector>& suitable_cells) { @@ -269,23 +216,21 @@ class Model if (config_.spread_schedule()[step]) { simulation_.generate( dispersers, + established_dispersers, infected, config_.weather, weather_coefficient, config_.reproductive_rate, suitable_cells); - DispersalKernel dispersal_kernel( - create_natural_kernel(dispersers), - create_anthro_kernel(dispersers), - config_.use_anthropogenic_kernel, - config_.percent_natural_dispersal); + auto dispersal_kernel = kernel_factory_(config_, dispersers, network); auto overpopulation_kernel = - create_overpopulation_movement_kernel(dispersers); + create_overpopulation_movement_kernel(dispersers, network); simulation_.disperse_and_infect( step, dispersers, + established_dispersers, susceptible, exposed, infected, diff --git a/inst/include/natural_anthropogenic_kernel.hpp b/inst/include/natural_anthropogenic_kernel.hpp index 0a1beb18..d7cbd321 100644 --- a/inst/include/natural_anthropogenic_kernel.hpp +++ b/inst/include/natural_anthropogenic_kernel.hpp @@ -18,6 +18,7 @@ #include "kernel_types.hpp" +#include #include #include #include @@ -38,27 +39,32 @@ namespace pops { * Bernoulli distribution is used to decide between the natural and * anthropogenic distance kernel. The anthropogenic distance dispersal can be also * competely disabled. + * + * Instances cannot be copied due to contained std::unique_ptr, but they can be + * moved, so returning new objects from functions works. This is enforced by the + * compiler. */ + template class NaturalAnthropogenicDispersalKernel { protected: bool use_anthropogenic_kernel_; - NaturalKernelType natural_kernel_; - AnthropogenicKernelType anthropogenic_kernel_; + std::unique_ptr natural_kernel_; + std::unique_ptr anthropogenic_kernel_; std::bernoulli_distribution bernoulli_distribution; public: NaturalAnthropogenicDispersalKernel( - const NaturalKernelType& natural_kernel, - const AnthropogenicKernelType& anthropogenic_kernel, + std::unique_ptr natural_kernel, + std::unique_ptr anthropogenic_kernel, bool use_anthropogenic_kernel, double percent_natural_dispersal) : use_anthropogenic_kernel_(use_anthropogenic_kernel), // Here we initialize all distributions, // although we won't use all of them. - natural_kernel_(natural_kernel), - anthropogenic_kernel_(anthropogenic_kernel), + natural_kernel_(std::move(natural_kernel)), + anthropogenic_kernel_(std::move(anthropogenic_kernel)), // use bernoulli distribution to act as the sampling with prob(gamma,1-gamma) bernoulli_distribution(percent_natural_dispersal) {} @@ -69,12 +75,12 @@ class NaturalAnthropogenicDispersalKernel std::tuple operator()(Generator& generator, int row, int col) { // switch in between the supported kernels - if (!use_anthropogenic_kernel_ || bernoulli_distribution(generator)) { - return natural_kernel_(generator, row, col); - } - else { - return anthropogenic_kernel_(generator, row, col); + if (!use_anthropogenic_kernel_ + || !anthropogenic_kernel_->is_cell_eligible(row, col) + || bernoulli_distribution(generator)) { + return natural_kernel_->operator()(generator, row, col); } + return anthropogenic_kernel_->operator()(generator, row, col); } /*! \copydoc RadialDispersalKernel::supports_kernel() diff --git a/inst/include/natural_kernel.hpp b/inst/include/natural_kernel.hpp new file mode 100644 index 00000000..f8d2d90a --- /dev/null +++ b/inst/include/natural_kernel.hpp @@ -0,0 +1,84 @@ +/* + * PoPS model - main disperal kernel + * + * Copyright (C) 2019-2021 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_NATURAL_KERNEL_HPP +#define POPS_NATURAL_KERNEL_HPP + +#include "radial_kernel.hpp" +#include "deterministic_kernel.hpp" +#include "uniform_kernel.hpp" +#include "neighbor_kernel.hpp" +#include "kernel_types.hpp" +#include "kernel_base.hpp" +#include "config.hpp" + +#include + +namespace pops { + +/** + * @brief Create natural kernel from configuration + * + * Kernel parameters are taken from the configuration. + * + * @param dispersers The disperser raster (reference, for deterministic kernel) + * @return Created kernel + */ +template +std::unique_ptr> +create_natural_kernel(const Config& config, const IntegerRaster& dispersers) +{ + auto natural_kernel = kernel_type_from_string(config.natural_kernel_type); + if (natural_kernel == DispersalKernelType::Uniform) { + using Kernel = DynamicWrapperKernel; + // This can be std::make_unique in C++14. + return std::unique_ptr(new Kernel(config.rows, config.cols)); + } + else if (natural_kernel == DispersalKernelType::DeterministicNeighbor) { + using Kernel = + DynamicWrapperKernel; + return std::unique_ptr( + new Kernel(direction_from_string(config.natural_direction))); + } + else if (config.deterministic) { + using Kernel = DynamicWrapperKernel< + DeterministicDispersalKernel, + Generator>; + return std::unique_ptr(new Kernel( + natural_kernel, + dispersers, + config.dispersal_percentage, + config.ew_res, + config.ns_res, + config.natural_scale, + config.shape)); + } + else { + using Kernel = + DynamicWrapperKernel, Generator>; + return std::unique_ptr(new Kernel( + config.ew_res, + config.ns_res, + natural_kernel, + config.natural_scale, + direction_from_string(config.natural_direction), + config.natural_kappa, + config.shape)); + } +} + +} // namespace pops + +#endif // POPS_NATURAL_KERNEL_HPP diff --git a/inst/include/neighbor_kernel.hpp b/inst/include/neighbor_kernel.hpp index de649956..a13cc9ef 100644 --- a/inst/include/neighbor_kernel.hpp +++ b/inst/include/neighbor_kernel.hpp @@ -98,6 +98,15 @@ class DeterministicNeighborDispersalKernel return std::make_tuple(row, col); } + /*! \copydoc RadialDispersalKernel::is_cell_eligible() + */ + bool is_cell_eligible(int row, int col) + { + UNUSED(row); + UNUSED(col); + return true; + } + /*! \copybrief RadialDispersalKernel::supports_kernel() */ static bool supports_kernel(const DispersalKernelType type) diff --git a/inst/include/network.hpp b/inst/include/network.hpp new file mode 100644 index 00000000..adf812fc --- /dev/null +++ b/inst/include/network.hpp @@ -0,0 +1,710 @@ +/* + * PoPS model - network dispersal kernel + * + * Copyright (C) 2020-2021 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_NETWORK_HPP +#define POPS_NETWORK_HPP + +#include "raster.hpp" +#include "utils.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace pops { + +/** + * Network structure and algorithms + * + * The network consists of nodes connected by edges. Edges are spatially represented + * as segments. Edges themselves don't carry cost and have meaning only as indicators of + * existence of a connection between nodes. Cost for each edge is a travel distance + * determined by advancing through cells in a segment. + * + * Nodes are the hop-on locations for dispersers. The dispersers can hop-off anywhere. + * + * The general workflow is contructing the object (with the constructor) and loading the + * data (with the load() function). Then the network is ready to be used for simulating + * trips over the network (with the travel() function). + * + * When the travel() function is used from a kernel, user of the network directly calls + * only the setup functions. + * + * The class exposes number of functions as public which are meant for testing or other + * special workflows. + * + * RasterIndex template parameter is an integral type used for indexing in the grid, + * specifically used as a type for rows and columns. + * + * Only non-const (write) functions for network are those loading the network, so an + * existing network can be shared by multiple kernel or model instances. + */ +template +class Network +{ +public: + using NodeId = int; ///< Type for node IDs + using Statistics = std::map; ///< Type for summary statistics + + /** + * @brief Construct empty network. + * + * @param bbox Bounding box of the raster grid (in real world coordinates) + * @param ew_res East-west resolution of the raster grid + * @param ns_res North-south resolution of the raster grid + */ + Network(BBox bbox, double ew_res, double ns_res) + : bbox_(bbox), + ew_res_(ew_res), + ns_res_(ns_res), + max_row_(0), + max_col_(0), + distance_per_cell_((ew_res + ns_res) / 2) + { + std::tie(max_row_, max_col_) = xy_to_row_col(bbox_.east, bbox_.south); + } + + /** + * @brief Create an empty network not meant for futher use. + * + * This is useful when a network object is needed to contruct a kernel, but it will + * not be used in runtime. + * + * @return New Network object + */ + static Network null_network() + { + return Network(BBox(), 0, 0); + } + + /** + * @brief Convert real world coordinates to row and column in the raster grid. + * + * @param x The X coordinate (east-west or column direction) + * @param y The Y coordinate (north-south or row direction) + * @return Row and column + */ + std::pair xy_to_row_col(double x, double y) const + { + double col = (x - bbox_.west) / ew_res_; + double row = (bbox_.north - y) / ns_res_; + return {std::floor(row), std::floor(col)}; + } + + /** + * Coverts coordinates as xy_to_row_col(double, double) but converts from strings. + */ + std::pair + xy_to_row_col(const std::string& x, const std::string& y) const + { + try { + return xy_to_row_col(std::stod(x), std::stod(y)); + } + catch (const std::invalid_argument&) { + if (string_contains(x, '"') || string_contains(x, '\'')) { + throw std::invalid_argument( + std::string("Text for coordinate X cannot contain quotes: ") + x); + } + if (string_contains(y, '"') || string_contains(y, '\'')) { + throw std::invalid_argument( + std::string("Text for coordinate Y cannot contain quotes: ") + y); + } + if (x.empty()) { + throw std::invalid_argument( + "Text for coordinate X cannot be an empty string, Y is " + y); + } + if (y.empty()) { + throw std::invalid_argument( + "Text for coordinate Y cannot be an empty string, X is " + y); + } + throw std::invalid_argument( + std::string("Text cannot be converted to X and Y coordinates: X: ") + x + + " Y: " + y); + } + catch (const std::out_of_range&) { + throw std::out_of_range( + std::string("Numerical value too large for X or Y coordinate: X: ") + x + + " Y: " + y); + } + } + + /** + * @brief Test if coordinates XY are in the bounding box. + * @param x Real world coordinate X + * @param y Real world coordinate Y + * @return True if XY is in the bounding box, false otherwise. + */ + bool xy_out_of_bbox(double x, double y) const + { + return x > bbox_.east || x < bbox_.west || y > bbox_.north || y < bbox_.south; + } + + /** + * @brief Test if a cell is in the bounding box. + * @param row Row in the raster grid + * @param col Column in the raster grid + * @return True if cells is in the bounding box, false otherwise. + */ + bool row_col_out_of_bbox(RasterIndex row, RasterIndex col) const + { + return row > max_row_ || row < 0 || col > max_col_ || col < 0; + } + + /** + * @brief Test if a cell is in the bounding box. + * @param cell Pair consisting of a row and column indices + * @return True if cell is in the bounding box, false otherwise. + */ + bool cell_out_of_bbox(const std::pair& cell) const + { + return cell.first > max_row_ || cell.first < 0 || cell.second > max_col_ + || cell.second < 0; + } + + /** + * @brief Load network from an input stream. + * + * @param stream Input stream containing text records for network + * @param allow_empty True if the loaded network can be empty + * + * Network stream is text in a custom format resembling CSV geared towards + * parsing, not readability, with rows `node_id_1,node_id_2,X1;Y1;X2;Y2;X3;Y3;...`. + * Each line represents an edge (pair of connected nodes) and the spatial + * representation of the edge (here called a segment). Node coordinates are taken + * from the first and last coordinate pair in the segment. + * + * The function can take any input stream which behaves like std::ifstream or + * std::istringstream. These are also the two expected streams to be used + * (for reading from a file and reading from a string, respectively). + * + * A simple code for reading from files without any checking may look like this: + * + * ``` + * Network network(...); + * std::string filename = "network.txt" + * std::ifstream stream{filename}; + * network.load(stream); + * ``` + * + * The input network is a list of segments where each segment consists of a pairs of + * nodes and list of coordinates. In the terminology of the this class, the pair of + * nodes is called an edge and the edge with coordinates (or rows and columns) is + * called a segment. + * + * The input is used as-is. The network is agnostic + * towards how the segments look like in terms of crossing each other or other + * nodes. + * + * All XY coordinates of segment points including nodes coordinates are converted + * to row and column using bounding box and resolution. + * + * Coordinates for nodes are taken from the beginning and the end of each segment, + * so input segment coordinates should include the node coordinates. + * + * Only edges (segments) with both nodes in the bounding box + * are considered, i.e., input network is reduced to the bounding box during + * reading, but clipping happens on the level of whole edges, not in the middle of + * a segment. + * + * Standalone nodes (without any connection) are theoretically allowed in the + * internal representation of the netwrok and handled + * as no movement from the source cell, but the input always needs to contain an + * edge. + */ + template + void load(InputStream& stream, bool allow_empty = false) + { + load_segments(stream); + if (node_matrix_.empty()) { + if (allow_empty) + return; + else + throw std::runtime_error("Network: No nodes within the extend"); + } + } + + /** + * @brief Get a list of nodes at a given cell + * + * Returns a const reference to an internally stored set. If there are no nodes + * at a given cell, a reference to an empty set is returned. + * + * @param row Row in the raster grid + * @param col Column in the raster grid + * @return List of node IDs as a set + */ + const std::set& get_nodes_at(RasterIndex row, RasterIndex col) const + { + auto it = nodes_by_row_col_.find(std::make_pair(row, col)); + if (it != nodes_by_row_col_.end()) + return it->second; + // If not found, return an empty set (const reference to local static). + static std::set empty; + return empty; + } + + /** + * @brief Get a randomly selected node at a given cell + * + * There can be more than one node at a given cell. This function selects one + * of these nodes randomly. + * + * @param row Row in the raster grid + * @param col Column in the raster grid + * @param generator Random number generator + * @return Node ID + * @throws std::invalid_argument if there is no node at a given cell + */ + template + NodeId + get_random_node_at(RasterIndex row, RasterIndex col, Generator& generator) const + { + const auto& nodes = get_nodes_at(row, col); + auto num_nodes = nodes.size(); + if (num_nodes == 1) { + return *nodes.begin(); + } + else if (num_nodes > 1) { + return pick_random_item(nodes, generator); + } + else { + throw std::invalid_argument("No nodes at a given row and column"); + } + } + + /** + * @brief Test if the network has a node at a given row and column + * @param row Row + * @param col Column + * @return True if there is at least one node, false otherwise + */ + bool has_node_at(RasterIndex row, RasterIndex col) const + { + // Replace by contains in C++20. + // return nodes_by_row_col_.count(std::make_pair(row, col)) > 0; + auto it = nodes_by_row_col_.find(std::make_pair(row, col)); + if (it == nodes_by_row_col_.end()) + return false; + return true; + } + + /** + * @brief Get row and column for a node. + * @param node Node id to get the coordinates for + * @return Row and column pair + */ + std::pair get_node_row_col(NodeId node) const + { + for (const auto& item : nodes_by_row_col_) { + if (container_contains(item.second, node)) + return item.first; + } + throw std::invalid_argument("No node with a given id"); + } + + /** + * Travel given distance in the network from given row and column. + * + * All previously visited nodes are tracked and, if possible, excluded + * from further traveling. + * + * @returns Final row and column pair + */ + template + std::tuple travel( + RasterIndex row, RasterIndex col, double distance, Generator& generator) const + { + // We assume there is a node here, i.e., that we are made decision + // to use this kernel knowing there is a node. + auto node_id = get_random_node_at(row, col, generator); + std::set visited_nodes; + while (distance >= 0) { + auto next_node_id = next_node(node_id, visited_nodes, generator); + visited_nodes.insert(node_id); + // If there is no segment from the node, return the start cell. + if (next_node_id == node_id) + return std::make_tuple(row, col); + auto segment = get_segment(node_id, next_node_id); + // nodes may need special handling + for (const auto& cell : segment) { + distance -= distance_per_cell_; + if (distance <= 0) { + return cell; + // Given the while condition, this subsequently ends the while loop + // as well. + // break; + } + } + node_id = next_node_id; + } + throw std::invalid_argument("Distance must be greater than or equal to zero"); + } + + /** + * @brief Get all nodes as vector of all ids with their row and column. + * + * If there is more than one node at a given cell (row and column), + * each of these nodes is returned as a separate item in the list. + * + * This translates the internal representation and returns a new object. + * + * @return A vector of pairs of id and row and col pair. + */ + std::vector>> + get_all_nodes() const + { + std::vector>> nodes; + nodes.reserve(nodes_by_row_col_.size()); // It will be at least this big. + for (const auto& item : nodes_by_row_col_) { + RasterIndex row = item.first.first; + RasterIndex col = item.first.second; + for (const auto node_id : item.second) { + nodes.emplace_back(node_id, std::make_pair(row, col)); + } + } + return nodes; + } + + /** + * @brief Collect statistics about the network + * + * Computes statistics such as number of nodes and segments. + * + * @return Associative array with statistics + */ + Statistics collect_statistics() const + { + std::map stats; + std::set node_ids; + for (const auto& item : nodes_by_row_col_) { + for (const auto node_id : item.second) { + node_ids.insert(node_id); + } + } + stats["num_nodes"] = node_ids.size(); + // We store segements in both directions, so each segment is stored twice. + stats["num_segments"] = node_matrix_.size() / 2; + std::set nodes_with_segments; + for (const auto& item : node_matrix_) { + nodes_with_segments.insert(item.first); + } + stats["num_nodes_with_segments"] = nodes_with_segments.size(); + // TODO: Only count the nodes here. Output them in the dump network output. + int num_standalone_nodes = 0; + for (NodeId node_id : node_ids) { + if (nodes_with_segments.find(node_id) == nodes_with_segments.end()) { + stats["standalone_node_" + std::to_string(++num_standalone_nodes)] = + node_id; + } + } + stats["num_standalone_nodes"] = num_standalone_nodes; + // TODO: This can be faster. The list is a ordered set and we iterate over + // elements above. + const auto minmax_node_id = minmax_element(node_ids.begin(), node_ids.end()); + stats["min_node_id"] = *minmax_node_id.first; + stats["max_node_id"] = *minmax_node_id.second; + return stats; + } + + /** + * @brief Output network to a stream. + * + * Writes statistics and different views of the network as a YAML file. + */ + template + void dump_yaml(OutputStream& stream) const + { + stream << "network:\n"; + stream << " statistics:\n"; + auto stats = collect_statistics(); + for (const auto& item : stats) { + stream << " " << item.first << ": " << item.second << "\n"; + } + stream << " extent:\n"; + stream << " north: " << bbox_.north << "\n"; + stream << " south: " << bbox_.south << "\n"; + stream << " east: " << bbox_.east << "\n"; + stream << " west: " << bbox_.west << "\n"; + stream << " resolution:\n"; + stream << " ns: " << ns_res_ << "\n"; + stream << " ew: " << ew_res_ << "\n"; + stream << " raster:\n"; + RasterIndex min_row; + RasterIndex min_col; + RasterIndex max_row; + RasterIndex max_col; + std::tie(min_row, min_col) = + xy_to_row_col(bbox_.west + ew_res_ / 2, bbox_.north - ns_res_ / 2); + std::tie(max_row, max_col) = + xy_to_row_col(bbox_.east - ew_res_ / 2, bbox_.south + ew_res_ / 2); + stream << " min_row: " << min_row << "\n"; + stream << " min_col: " << min_col << "\n"; + stream << " max_row: " << max_row << "\n"; + stream << " max_col: " << max_col << "\n"; + stream << " num_rows: " << max_row - min_row + 1 << "\n"; + stream << " num_cols: " << max_col - min_col + 1 << "\n"; + stream << " cost:\n"; + stream << " distance_per_cell: " << distance_per_cell_ << "\n"; + std::set> edges; + for (const auto& item : node_matrix_) { + for (const auto& node_id : item.second) { + edges.emplace(item.first, node_id); + } + } + stream << " edges:\n"; + for (const auto& item : edges) { + stream << " - [" << item.first << ", " << item.second << "]\n"; + } + stream << " nodes:\n"; + for (const auto& item : nodes_by_row_col_) { + auto row = item.first.first; + auto col = item.first.second; + for (const auto& node_id : item.second) { + stream << " - id: " << node_id << "\n"; + stream << " row: " << row << "\n"; + stream << " col: " << col << "\n"; + } + } + stream << " segments:\n"; + for (const auto& item : segments_by_nodes_) { + stream << " - start_node: " << item.first.first << "\n"; + stream << " end_node: " << item.first.second << "\n"; + stream << " cells: ["; + for (const auto& cell : item.second) { + stream << "[" << cell.first << ", " << cell.second << "], "; + } + stream << "]\n"; + } + } + +protected: + /** Node connections (edges) + * + * The matrix is sparse. One node can be connected to multiple nodes. + * Connection in both directions is stored explicitly (see load_segments() source + * code). + */ + using NodeMatrix = std::map>; + /** Cells connecting two nodes (segment between nodes) */ + using Segment = std::vector>; + /** Constant view of a segment (to iterate a segment in either direction) */ + using SegmentView = ContainerView; + /** Segments by nodes (edges) */ + using SegmentsByNodes = std::map, Segment>; + + /** + * @brief Convert string to node ID + * + * @param text String with node ID + * @return Node ID + */ + static NodeId node_id_from_text(const std::string& text) + { + try { + return std::stoi(text); + } + catch (const std::invalid_argument& err) { + if (string_contains(text, '"') || string_contains(text, '\'')) { + throw std::invalid_argument( + std::string("Text for a node ID cannot contain quotes " + "(only digits are allowed): ") + + text); + } + if (text.empty()) { + throw std::invalid_argument( + "Text for a node ID cannot be an empty string"); + } + else { + throw std::invalid_argument( + std::string("Text cannot be converted to a node ID " + "(only digits are allowed): ") + + text); + } + } + catch (const std::out_of_range& err) { + throw std::out_of_range( + std::string("Numerical value too large for a node ID: ") + text); + } + } + + /** + * @brief Read segments from a stream. + * + * @param stream Input stream with segments + * @param node_ids Container with node IDs to use + */ + template + void load_segments(InputStream& stream) + { + std::string line; + while (std::getline(stream, line)) { + std::istringstream line_stream{line}; + char delimeter{','}; + std::string node_1_text; + std::getline(line_stream, node_1_text, delimeter); + std::string node_2_text; + std::getline(line_stream, node_2_text, delimeter); + auto node_1_id = node_id_from_text(node_1_text); + auto node_2_id = node_id_from_text(node_2_text); + if (node_1_id < 1 || node_2_id < 1) { + std::runtime_error(std::string( + "Node ID must be greater than zero: " + node_1_text + " " + + node_2_text)); + } + if (node_1_id == node_2_id) { + std::runtime_error( + std::string("Edge cannot begin and end with the same node: ") + + node_1_text + " " + node_2_text); + } + std::string segment_text; + std::getline(line_stream, segment_text, delimeter); + std::istringstream segment_stream{segment_text}; + char in_cell_delimeter{';'}; + std::string x_coord_text; + std::string y_coord_text; + Segment segment; + while (std::getline(segment_stream, x_coord_text, in_cell_delimeter) + && std::getline(segment_stream, y_coord_text, in_cell_delimeter)) { + // The same cell is possibly repeated if raster resolution is lower than + // the detail ("resolution") of each segment, so we check if the last + // added coordinate pair converted to same cell. + auto new_point = xy_to_row_col(x_coord_text, y_coord_text); + if (segment.empty() || segment.back() != new_point) + segment.emplace_back(new_point); + } + // If either node of the segment is not in the extent, skip the segment. + // This means that a segment is ignored even if one of the nodes and + // significant portion of the segment is in the area of iterest. + // TODO: Part of the segment may still be out, so that needs to be checked. + // Cut to extend + if (cell_out_of_bbox(segment.front()) || cell_out_of_bbox(segment.back())) + continue; + // We are done with the segment data, so we move them to the attribute. + segments_by_nodes_.emplace( + std::make_pair(node_1_id, node_2_id), std::move(segment)); + } + + for (const auto& node_segment : segments_by_nodes_) { + nodes_by_row_col_[node_segment.second.front()].insert( + node_segment.first.first); + nodes_by_row_col_[node_segment.second.back()].insert( + node_segment.first.second); + node_matrix_[node_segment.first.first].push_back(node_segment.first.second); + node_matrix_[node_segment.first.second].push_back(node_segment.first.first); + } + } + + /** + * @brief Get a segment between the two given nodes + * + * For a given pair of nodes, the function finds the connecting segment + * and returns a view of the segment oriented according to the order of the nodes. + * + * @param start Start node id + * @param end End node id + * @returns View of the segment oriented from start to end + * @throws std::invalid_argument if there is no segment connecting the nodes + */ + SegmentView get_segment(NodeId start, NodeId end) const + { + for (const auto& item : segments_by_nodes_) { + const auto& key{item.first}; + if (key.first == start && key.second == end) + return SegmentView(item.second.cbegin(), item.second.cend()); + if (key.second == start && key.first == end) { + return SegmentView(item.second.crbegin(), item.second.crend()); + } + } + throw std::invalid_argument(std::string( + "No segment for given nodes: " + std::to_string(start) + " " + + std::to_string(end))); + } + + /** + * @brief Get nodes connected by an edge to a given node. + * + * @param node Node to get connections from + * @return List of connected nodes + */ + const std::vector& nodes_connected_to(NodeId node) const + { + return node_matrix_.at(node); + } + + /** + * @brief Pick a next node from the given node. + * + * If there is more than one edge leading from the given node, a random node is + * picked. If there are no edges leading from the given node, the node is returned. + * + * The random node is picked from candidate nodes which are nodes connected to + * a given node. The candidate nodes which are in the *ignore* list are excluded + * from the random selection. If all candiate nodes are in the *ignore* list, + * the *ignore* list is ignored and all candidate nodes are used. + * + * The function always returns a node to go to even if it means going to an ignored + * node or returning the value of the *node* parameter. + */ + template + NodeId + next_node(NodeId node, const std::set& ignore, Generator& generator) const + { + // Get all candidate nodes. + const auto& all_nodes = nodes_connected_to(node); + + // Resolve disconnected node and dead end cases. + auto num_nodes = all_nodes.size(); + if (!num_nodes) + return node; + else if (num_nodes == 1) + return all_nodes[0]; + + // Filter out the ignored nodes. + std::vector nodes; + std::back_insert_iterator> back_it(nodes); + std::remove_copy_if( + all_nodes.begin(), all_nodes.end(), back_it, [&ignore](NodeId id) { + return container_contains(ignore, id); + }); + + // Pick a random node. Fallback to all candidate nodes if all are in ignore. + num_nodes = nodes.size(); + if (!num_nodes) + return pick_random_item(all_nodes, generator); + else if (num_nodes == 1) + return nodes[0]; + return pick_random_item(nodes, generator); + } + + BBox bbox_; ///< Bounding box of the network grid in real world coordinates + double ew_res_; ///< East-west resolution of the grid + double ns_res_; ///< North-south resolution of the grid + RasterIndex max_row_; ///< Maximum row index in the grid + RasterIndex max_col_; ///< Maximum column index in the grid + double distance_per_cell_; ///< Distance to travel through one cell (cost) + /** Node IDs stored by row and column (multiple nodes per cell) */ + std::map, std::set> nodes_by_row_col_; + NodeMatrix node_matrix_; ///< List of node neighbors by node ID (edges) + SegmentsByNodes segments_by_nodes_; ///< Lists of cells connecting nodes +}; + +} // namespace pops + +#endif // POPS_NETWORK_HPP diff --git a/inst/include/network_kernel.hpp b/inst/include/network_kernel.hpp new file mode 100644 index 00000000..e52e79f2 --- /dev/null +++ b/inst/include/network_kernel.hpp @@ -0,0 +1,96 @@ +/* + * PoPS model - network dispersal kernel + * + * Copyright (C) 2020-2021 by the authors. + * + * Authors: Vaclav Petras (wenzeslaus gmail com) + * + * The code contained herein is licensed under the GNU General Public + * License. You may obtain a copy of the GNU General Public License + * Version 2 or later at the following locations: + * + * http://www.opensource.org/licenses/gpl-license.html + * http://www.gnu.org/copyleft/gpl.html + */ + +#ifndef POPS_NETWORK_KERNEL_HPP +#define POPS_NETWORK_KERNEL_HPP + +#include "kernel_types.hpp" +#include "network.hpp" + +namespace pops { + +/*! + * @brief Dispersal kernel for dispersal over a network. + * + * Network node must be present in the cell to start traveling, so is_cell_eligible() + * needs to be called first to see if the kernel can be used with the given row and + * column. + */ +template +class NetworkDispersalKernel +{ +public: + /** + * @brief Create kernel. + * + * The kernel assumes that the *network* is already initialized. It does not modify + * the network. + * + * The *min_distance* and *max_distance* parameters are used as a range for uniform + * real distribution which determines the travel distance through the network for + * one trip. + * + * @param network Existing network + * @param min_distance Minimum travel distance + * @param max_distance Maximum travel distance + */ + NetworkDispersalKernel( + const Network& network, double min_distance, double max_distance) + : network_(network), distance_distribution_(min_distance, max_distance) + {} + + /*! \copybrief RadialDispersalKernel::operator()() + * + * @throws std::invalid_argument if there is no network node at a given *row* + * and *col* (can be checked with is_cell_eligible() beforehand) + */ + template + std::tuple operator()(Generator& generator, int row, int col) + { + double distance = distance_distribution_(generator); + std::tie(row, col) = network_.travel(row, col, distance, generator); + + return std::make_tuple(row, col); + } + + /** + * @brief Test if cell is eligible to be used with the kernel. + * + * @param row Row to be used with the kernel + * @param col Column to be used with the kernel + * @return true if cell can be used, false otherwise + */ + bool is_cell_eligible(int row, int col) + { + return network_.has_node_at(row, col); + } + + /*! \copydoc RadialDispersalKernel::supports_kernel() + */ + static bool supports_kernel(const DispersalKernelType type) + { + return type == DispersalKernelType::Network; + } + +protected: + /** Reference to the network */ + const Network& network_; + /** Travel distance distribution */ + std::uniform_real_distribution distance_distribution_; +}; + +} // namespace pops + +#endif // POPS_NETWORK_KERNEL_HPP diff --git a/inst/include/radial_kernel.hpp b/inst/include/radial_kernel.hpp index 82341878..75666fe3 100644 --- a/inst/include/radial_kernel.hpp +++ b/inst/include/radial_kernel.hpp @@ -147,9 +147,9 @@ class RadialDispersalKernel /*! Generates a new position for the spread. * - * The randomness is based on the *generator*. The result may depend - * on previous calls of this operator (see e.g. - * `std::cauchy_distribution::reset()`). + * The randomness is based on the *generator*, but the result may depend + * on previous calls of this operator (see e.g. documentation of the underlying + * `std::cauchy_distribution`, specifically the `reset()` function). * Parameters *row* and *col* are row and column position of the * current disperser. The generated position will be relative to it. */ @@ -204,6 +204,15 @@ class RadialDispersalKernel return std::make_tuple(row, col); } + /*! Returns true if kernel can be used with a given cell. + */ + bool is_cell_eligible(int row, int col) + { + UNUSED(row); + UNUSED(col); + return true; + } + /*! Returns true if the kernel class support a given kernel type * * \warning This function is experimental and may be removed or diff --git a/inst/include/raster.hpp b/inst/include/raster.hpp index f1e837d0..cba00d61 100644 --- a/inst/include/raster.hpp +++ b/inst/include/raster.hpp @@ -130,7 +130,8 @@ class Raster { cols_ = other.cols_; rows_ = other.rows_; - data_ = new Number[cols_ * rows_]{value}; + data_ = new Number[cols_ * rows_]; + std::fill_n(data_, cols_ * rows_, value); } Raster(Raster&& other) : owns_(other.owns_) @@ -152,7 +153,8 @@ class Raster { this->cols_ = cols; this->rows_ = rows; - this->data_ = new Number[cols_ * rows_]{value}; + this->data_ = new Number[cols_ * rows_]; + std::fill_n(data_, cols_ * rows_, value); } /*! Use existing data storage @@ -230,7 +232,7 @@ class Raster } template - void for_each(UnaryOperation op) + void for_each(UnaryOperation op) const { std::for_each(data_, data_ + (cols_ * rows_), op); } diff --git a/inst/include/scheduling.hpp b/inst/include/scheduling.hpp index 7773d080..ec56827a 100644 --- a/inst/include/scheduling.hpp +++ b/inst/include/scheduling.hpp @@ -308,7 +308,7 @@ class Scheduler } /** * @brief Prints schedule for debugging purposes. - * @param vector of bools to print along the steps + * @param schedule vector of bools to print along the steps */ void debug_schedule(std::vector& schedule) const { @@ -443,6 +443,8 @@ inline std::vector schedule_from_string( } else if (frequency == "every_n_steps" && n > 0) return scheduler.schedule_action_nsteps(n); + else if (frequency == "every_step" || frequency == "time_step") + return scheduler.schedule_action_nsteps(1); else throw std::invalid_argument("Invalid value of output frequency"); } diff --git a/inst/include/switch_kernel.hpp b/inst/include/switch_kernel.hpp index 7b5c5707..cf74bfe4 100644 --- a/inst/include/switch_kernel.hpp +++ b/inst/include/switch_kernel.hpp @@ -1,7 +1,7 @@ /* * PoPS model - disperal kernels * - * Copyright (C) 2019 - 2020 by the authors. + * Copyright (C) 2019-2021 by the authors. * * Authors: Vaclav Petras (wenzeslaus gmail com) * @@ -20,7 +20,9 @@ #include "deterministic_kernel.hpp" #include "uniform_kernel.hpp" #include "neighbor_kernel.hpp" +#include "network_kernel.hpp" #include "kernel_types.hpp" +#include "utils.hpp" namespace pops { @@ -33,7 +35,7 @@ namespace pops { * its call in the function call operator, and extend the * supports_kernel() function. */ -template +template class SwitchDispersalKernel { protected: @@ -42,6 +44,7 @@ class SwitchDispersalKernel DeterministicDispersalKernel deterministic_kernel_; UniformDispersalKernel uniform_kernel_; DeterministicNeighborDispersalKernel deterministic_neighbor_kernel_; + NetworkDispersalKernel network_kernel_; bool deterministic_; public: @@ -50,6 +53,7 @@ class SwitchDispersalKernel const RadialDispersalKernel& radial_kernel, const DeterministicDispersalKernel& deterministic_kernel, const UniformDispersalKernel& uniform_kernel, + const NetworkDispersalKernel& network_kernel, const DeterministicNeighborDispersalKernel& deterministic_neighbor_kernel = DeterministicNeighborDispersalKernel(Direction::None), const bool deterministic = false) @@ -60,6 +64,7 @@ class SwitchDispersalKernel deterministic_kernel_(deterministic_kernel), uniform_kernel_(uniform_kernel), deterministic_neighbor_kernel_(deterministic_neighbor_kernel), + network_kernel_(network_kernel), deterministic_(deterministic) {} @@ -75,6 +80,9 @@ class SwitchDispersalKernel else if (dispersal_kernel_type_ == DispersalKernelType::DeterministicNeighbor) { return deterministic_neighbor_kernel_(generator, row, col); } + else if (dispersal_kernel_type_ == DispersalKernelType::Network) { + return network_kernel_(generator, row, col); + } else if (deterministic_) { return deterministic_kernel_(generator, row, col); } @@ -83,6 +91,27 @@ class SwitchDispersalKernel } } + bool is_cell_eligible(int row, int col) + { + // switch in between the supported kernels + if (dispersal_kernel_type_ == DispersalKernelType::Uniform) { + // TODO: Individual kernels should support this. + return true; + } + else if (dispersal_kernel_type_ == DispersalKernelType::DeterministicNeighbor) { + return true; + } + else if (dispersal_kernel_type_ == DispersalKernelType::Network) { + return network_kernel_.is_cell_eligible(row, col); + } + else if (deterministic_) { + return true; + } + else { + return true; + } + } + /*! \copydoc RadialDispersalKernel::supports_kernel() */ static bool supports_kernel(const DispersalKernelType type) diff --git a/inst/include/treatments.hpp b/inst/include/treatments.hpp index 281da2c9..fdc598be 100644 --- a/inst/include/treatments.hpp +++ b/inst/include/treatments.hpp @@ -350,9 +350,9 @@ class Treatments * \param map treatment raster * \param start_date date when treatment is applied * \param num_days for simple treatments should be 0, otherwise number of days host - * is resistant \param treatment_application if efficiency < 100% how should it be - * applied to infected/susceptible \param increase_by_step function to increase - * simulation step + * is resistant + * \param treatment_application if efficiency < 100% how should it be + * applied to infected/susceptible */ void add_treatment( const FloatRaster& map, @@ -380,8 +380,12 @@ class Treatments * * \param current simulation step * \param infected raster of infected host + * \param exposed Exposed hosts per cohort * \param susceptible raster of susceptible host * \param resistant raster of resistant host + * \param total_hosts All host individuals in the area (I+E+S in the cell) + * \param suitable_cells List of indices of cells with hosts + * * \return true if any management action was necessary */ bool manage( @@ -416,6 +420,8 @@ class Treatments * \brief Separately manage mortality infected cohorts * \param current simulation step * \param infected raster of infected host + * \param suitable_cells List of indices of cells with hosts + * * \return true if any management action was necessary */ bool manage_mortality( diff --git a/inst/include/uniform_kernel.hpp b/inst/include/uniform_kernel.hpp index 48650866..cec45dbd 100644 --- a/inst/include/uniform_kernel.hpp +++ b/inst/include/uniform_kernel.hpp @@ -62,6 +62,15 @@ class UniformDispersalKernel return std::make_tuple(row, col); } + /*! \copydoc RadialDispersalKernel::is_cell_eligible() + */ + bool is_cell_eligible(int row, int col) + { + UNUSED(row); + UNUSED(col); + return true; + } + /*! \copydoc RadialDispersalKernel::supports_kernel() */ static bool supports_kernel(const DispersalKernelType type) diff --git a/inst/include/utils.hpp b/inst/include/utils.hpp index dc9e0cb7..613873ec 100644 --- a/inst/include/utils.hpp +++ b/inst/include/utils.hpp @@ -4,7 +4,7 @@ /* * PoPS model - general utility functions (unrelated to the model) * - * Copyright (C) 2020 by the authors. + * Copyright (C) 2020-2021 by the authors. * * Authors: Vaclav Petras * @@ -35,12 +35,151 @@ #define M_PI 3.14159265358979323846 #define PI M_PI +#include #include +#include + +/** + * Return true if _container_ contains _value_. + */ +template +bool container_contains(const Container& container, const Value& value) +{ + return container.find(value) != container.end(); +} + +/** + * Return true if _string_ contains _value_ (character or substring). + */ +template +bool string_contains(const String& string, const Value& value) +{ + return string.find(value) != String::npos; +} + +// Replace by ranges::shuffle in C++20. +/** + * Reorder items in container. + */ +template +void shuffle_container(Container& container, Generator& generator) +{ + std::shuffle(container.begin(), container.end(), generator); +} + +/** + * Select a random item from a container. + * + * May be slow if it takes a long time to increase the iterator (e.g., for std::set) and + * the container is large. + */ +template +typename Container::value_type +pick_random_item(const Container& container, Generator& generator) +{ + // Replace .size() call by std::size in C++17. + std::uniform_int_distribution distribution(0, container.size() - 1); + auto index = distribution(generator); + // For small containers, this is expected to be fast for both sets and vectors. + return *std::next(container.begin(), index); +} + +/** + * \brief A const iterator which encapsulates either forward or reverse iterator. + * + * Uses the iterator was passed in the constructor, but it can contain either + * forward or reverse iterator, so it allows writing direction agnostic code + * where the direction is determined in the runtime. + */ +template +class ConstEitherWayIterator +{ +public: + ConstEitherWayIterator(typename Container::const_iterator it) + : is_forward_(true), forward_it_(it) + {} + ConstEitherWayIterator(typename Container::const_reverse_iterator it) + : is_forward_(false), reverse_it_(it) + {} + ConstEitherWayIterator& operator++() + { + if (is_forward_) + ++forward_it_; + else + ++reverse_it_; + return *this; + } + const typename Container::value_type& operator*() const + { + if (is_forward_) + return *forward_it_; + return *reverse_it_; + } + bool operator!=(const ConstEitherWayIterator& other) + { + if (is_forward_) + return forward_it_ != other.forward_it_; + return reverse_it_ != other.reverse_it_; + } + +private: + // Can be improved by std::optional or any in C++17. + bool is_forward_; + typename Container::const_iterator forward_it_; + typename Container::const_reverse_iterator reverse_it_; +}; + +/** + * \brief A const view of a container which can be iterated. + * + * Depending on which iterators are passed in the constructor, the view + * is in the forward direction of the container or in reverse. The view can also be + * a subset. The view supports only const iterators. + */ +template +class ContainerView +{ +public: + ContainerView( + typename Container::const_iterator first, + typename Container::const_iterator last) + : begin_(first), end_(last) + {} + ContainerView( + typename Container::const_reverse_iterator first, + typename Container::const_reverse_iterator last) + : begin_(first), end_(last) + {} + ConstEitherWayIterator begin() const + { + return begin_; + } + ConstEitherWayIterator end() const + { + return end_; + } + +private: + ConstEitherWayIterator begin_; + ConstEitherWayIterator end_; +}; typedef std::tuple BBoxInt; typedef std::tuple BBoxFloat; typedef std::tuple BBoxBool; +// C++20 NumericType +template +struct BBox +{ + Number north; + Number south; + Number east; + Number west; + + BBox() : north(0), south(0), east(0), west(0) {} +}; + /*! Spread direction * * Spread, typically wind, direction. @@ -123,4 +262,33 @@ std::string quarantine_enum_to_string(Direction type) } } +/** + * Create a list of suitable cells in from a host raster. + * + * Suitable cell is defined as cell with value higher than zero, i.e., there is at least + * one host. + * + * Suitable cells datastructure is vector of vectors where the nested vector always has + * size equal to two and contains row and column index. The type was chosen to work well + * with Rcpp. + * + * First template parameter is the index type for the resulting sutibale cell indices. + * Second template parameter is deduced automatically from the function parameter. + */ +template +std::vector> find_suitable_cells(const RasterType& raster) +{ + std::vector> cells; + // The assumption is that the raster is sparse (otherwise we would not be doing + // this), so we have no number for the reserve method. + for (RasterIndex row = 0; row < raster.rows(); ++row) { + for (RasterIndex col = 0; col < raster.cols(); ++col) { + if (raster(row, col) > 0) { + cells.push_back({row, col}); + } + } + } + return cells; +} + #endif // POPS_UTILS_HPP diff --git a/man/auto_manage.Rd b/man/auto_manage.Rd index 4b29bd2e..f1573a6b 100644 --- a/man/auto_manage.Rd +++ b/man/auto_manage.Rd @@ -75,7 +75,10 @@ auto_manage( exposed_file = "", mask = NULL, write_outputs = "None", - output_folder_path = "" + output_folder_path = "", + network_min_distance = 0, + network_max_distance = 0, + network_filename = "" ) } \arguments{ @@ -332,6 +335,12 @@ that can't be managed in the auto_manage function.} \item{output_folder_path}{this is the full path with either / or \\ (e.g., "C:/user_name/desktop/pops_sod_2020_2023/outputs/")} + +\item{network_min_distance}{minimum time a propagule rides on the network} + +\item{network_max_distance}{maximum time a propagule rides on the network} + +\item{network_filename}{entire file path for the network file} } \value{ list of infected and susceptible per year diff --git a/man/calibrate.Rd b/man/calibrate.Rd index c9236031..46680699 100644 --- a/man/calibrate.Rd +++ b/man/calibrate.Rd @@ -73,6 +73,9 @@ calibrate( verbose = TRUE, write_outputs = "None", output_folder_path = "", + network_min_distance = 0, + network_max_distance = 0, + network_filename = "", use_distance = FALSE, use_rmse = FALSE ) @@ -360,6 +363,12 @@ Defaults if FALSE.} \item{output_folder_path}{this is the full path with either / or \\ (e.g., "C:/user_name/desktop/pops_sod_2020_2023/outputs/")} +\item{network_min_distance}{minimum time a propagule rides on the network} + +\item{network_max_distance}{maximum time a propagule rides on the network} + +\item{network_filename}{entire file path for the network file} + \item{use_distance}{Boolean if you want to compare distance between simulations and observations. Default is FALSE.} diff --git a/man/pops.Rd b/man/pops.Rd index 78f14192..25b0da1e 100644 --- a/man/pops.Rd +++ b/man/pops.Rd @@ -60,7 +60,10 @@ pops( leaving_percentage = 0, leaving_scale_coefficient = 1, exposed_file = "", - mask = NULL + mask = NULL, + network_min_distance = 0, + network_max_distance = 0, + network_filename = "" ) } \arguments{ @@ -275,6 +278,12 @@ the natural kernel (if applicable)} true negatives from comparisons (e.g. mask out lakes and oceans from statics if modeling terrestrial species). This can also be used to mask out areas that can't be managed in the auto_manage function.} + +\item{network_min_distance}{minimum time a propagule rides on the network} + +\item{network_max_distance}{maximum time a propagule rides on the network} + +\item{network_filename}{entire file path for the network file} } \value{ list of infected and susceptible per year diff --git a/man/pops_model.Rd b/man/pops_model.Rd index a21cb32c..f4bc145b 100644 --- a/man/pops_model.Rd +++ b/man/pops_model.Rd @@ -71,7 +71,11 @@ pops_model( use_overpopulation_movements = FALSE, overpopulation_percentage = 0, leaving_percentage = 0, - leaving_scale_coefficient = 1 + leaving_scale_coefficient = 1, + bbox = NULL, + network_min_distance = 0, + network_max_distance = 0, + network_filename = "" ) } \arguments{ @@ -282,6 +286,14 @@ is considered to be overpopulated} \item{leaving_scale_coefficient}{coefficient to multiply scale parameter of the natural kernel (if applicable)} + +\item{bbox}{bounding box for network kernel} + +\item{network_min_distance}{minimum time a propagule rides on the network} + +\item{network_max_distance}{maximum time a propagule rides on the network} + +\item{network_filename}{entire file path for the network file} } \value{ list of vector matrices of infected and susceptible hosts per diff --git a/man/pops_multirun.Rd b/man/pops_multirun.Rd index 332dfb87..fe0b7d35 100644 --- a/man/pops_multirun.Rd +++ b/man/pops_multirun.Rd @@ -64,7 +64,10 @@ pops_multirun( exposed_file = "", mask = NULL, write_outputs = "None", - output_folder_path = "" + output_folder_path = "", + network_min_distance = 0, + network_max_distance = 0, + network_filename = "" ) } \arguments{ @@ -291,6 +294,12 @@ that can't be managed in the auto_manage function.} \item{output_folder_path}{this is the full path with either / or \\ (e.g., "C:/user_name/desktop/pops_sod_2020_2023/outputs/")} + +\item{network_min_distance}{minimum time a propagule rides on the network} + +\item{network_max_distance}{maximum time a propagule rides on the network} + +\item{network_filename}{entire file path for the network file} } \value{ list of infected and susceptible per year diff --git a/man/validate.Rd b/man/validate.Rd index 9d38f8fb..57994a08 100644 --- a/man/validate.Rd +++ b/man/validate.Rd @@ -67,6 +67,9 @@ validate( write_outputs = "None", output_folder_path = "", point_file = "", + network_min_distance = 0, + network_max_distance = 0, + network_filename = "", use_distance = FALSE, use_configuration = FALSE ) @@ -293,6 +296,12 @@ the natural kernel (if applicable)} \item{point_file}{file for point comparison if not provided skips calculations} +\item{network_min_distance}{minimum time a propagule rides on the network} + +\item{network_max_distance}{maximum time a propagule rides on the network} + +\item{network_filename}{entire file path for the network file} + \item{use_distance}{Boolean if you want to compare distance between simulations and observations. Default is FALSE.} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 5719e522..03acb88b 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -14,8 +14,8 @@ Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif // pops_model_cpp -List pops_model_cpp(int random_seed, double lethal_temperature, int lethal_temperature_month, IntegerMatrix infected, IntegerMatrix total_exposed, std::vector exposed, IntegerMatrix susceptible, IntegerMatrix total_populations, IntegerMatrix total_hosts, std::vector mortality_tracker, IntegerMatrix mortality, IntegerMatrix quarantine_areas, std::vector treatment_maps, std::vector treatment_dates, std::vector pesticide_duration, IntegerMatrix resistant, std::vector> movements, std::vector movements_dates, std::vector temperature, std::vector weather_coefficient, List res, List rows_cols, double reproductive_rate, std::vector> spatial_indices, List season_month_start_end, List frequency_config, List bool_config, double mortality_rate, int mortality_time_lag, std::string start_date, std::string end_date, std::string treatment_method, std::string natural_kernel_type, std::string anthropogenic_kernel_type, double percent_natural_dispersal, double natural_distance_scale, double anthropogenic_distance_scale, std::string natural_dir, double natural_kappa, std::string anthropogenic_dir, double anthropogenic_kappa, int output_frequency_n, int quarantine_frequency_n, int spreadrate_frequency_n, int mortality_frequency_n, std::string model_type_, int latency_period, double establishment_probability, double dispersal_percentage, Nullable overpopulation_config); -static SEXP _PoPS_pops_model_cpp_try(SEXP random_seedSEXP, SEXP lethal_temperatureSEXP, SEXP lethal_temperature_monthSEXP, SEXP infectedSEXP, SEXP total_exposedSEXP, SEXP exposedSEXP, SEXP susceptibleSEXP, SEXP total_populationsSEXP, SEXP total_hostsSEXP, SEXP mortality_trackerSEXP, SEXP mortalitySEXP, SEXP quarantine_areasSEXP, SEXP treatment_mapsSEXP, SEXP treatment_datesSEXP, SEXP pesticide_durationSEXP, SEXP resistantSEXP, SEXP movementsSEXP, SEXP movements_datesSEXP, SEXP temperatureSEXP, SEXP weather_coefficientSEXP, SEXP resSEXP, SEXP rows_colsSEXP, SEXP reproductive_rateSEXP, SEXP spatial_indicesSEXP, SEXP season_month_start_endSEXP, SEXP frequency_configSEXP, SEXP bool_configSEXP, SEXP mortality_rateSEXP, SEXP mortality_time_lagSEXP, SEXP start_dateSEXP, SEXP end_dateSEXP, SEXP treatment_methodSEXP, SEXP natural_kernel_typeSEXP, SEXP anthropogenic_kernel_typeSEXP, SEXP percent_natural_dispersalSEXP, SEXP natural_distance_scaleSEXP, SEXP anthropogenic_distance_scaleSEXP, SEXP natural_dirSEXP, SEXP natural_kappaSEXP, SEXP anthropogenic_dirSEXP, SEXP anthropogenic_kappaSEXP, SEXP output_frequency_nSEXP, SEXP quarantine_frequency_nSEXP, SEXP spreadrate_frequency_nSEXP, SEXP mortality_frequency_nSEXP, SEXP model_type_SEXP, SEXP latency_periodSEXP, SEXP establishment_probabilitySEXP, SEXP dispersal_percentageSEXP, SEXP overpopulation_configSEXP) { +List pops_model_cpp(int random_seed, double lethal_temperature, int lethal_temperature_month, IntegerMatrix infected, IntegerMatrix total_exposed, std::vector exposed, IntegerMatrix susceptible, IntegerMatrix total_populations, IntegerMatrix total_hosts, std::vector mortality_tracker, IntegerMatrix mortality, IntegerMatrix quarantine_areas, std::vector treatment_maps, std::vector treatment_dates, std::vector pesticide_duration, IntegerMatrix resistant, std::vector> movements, std::vector movements_dates, std::vector temperature, std::vector weather_coefficient, List bbox, List res, List rows_cols, double reproductive_rate, std::vector> spatial_indices, List season_month_start_end, List frequency_config, List bool_config, double mortality_rate, int mortality_time_lag, std::string start_date, std::string end_date, std::string treatment_method, std::string natural_kernel_type, std::string anthropogenic_kernel_type, double percent_natural_dispersal, double natural_distance_scale, double anthropogenic_distance_scale, std::string natural_dir, double natural_kappa, std::string anthropogenic_dir, double anthropogenic_kappa, Nullable frequencies_n_config, std::string model_type_, int latency_period, double establishment_probability, double dispersal_percentage, Nullable overpopulation_config, Nullable network_config, Nullable network_data_config); +static SEXP _PoPS_pops_model_cpp_try(SEXP random_seedSEXP, SEXP lethal_temperatureSEXP, SEXP lethal_temperature_monthSEXP, SEXP infectedSEXP, SEXP total_exposedSEXP, SEXP exposedSEXP, SEXP susceptibleSEXP, SEXP total_populationsSEXP, SEXP total_hostsSEXP, SEXP mortality_trackerSEXP, SEXP mortalitySEXP, SEXP quarantine_areasSEXP, SEXP treatment_mapsSEXP, SEXP treatment_datesSEXP, SEXP pesticide_durationSEXP, SEXP resistantSEXP, SEXP movementsSEXP, SEXP movements_datesSEXP, SEXP temperatureSEXP, SEXP weather_coefficientSEXP, SEXP bboxSEXP, SEXP resSEXP, SEXP rows_colsSEXP, SEXP reproductive_rateSEXP, SEXP spatial_indicesSEXP, SEXP season_month_start_endSEXP, SEXP frequency_configSEXP, SEXP bool_configSEXP, SEXP mortality_rateSEXP, SEXP mortality_time_lagSEXP, SEXP start_dateSEXP, SEXP end_dateSEXP, SEXP treatment_methodSEXP, SEXP natural_kernel_typeSEXP, SEXP anthropogenic_kernel_typeSEXP, SEXP percent_natural_dispersalSEXP, SEXP natural_distance_scaleSEXP, SEXP anthropogenic_distance_scaleSEXP, SEXP natural_dirSEXP, SEXP natural_kappaSEXP, SEXP anthropogenic_dirSEXP, SEXP anthropogenic_kappaSEXP, SEXP frequencies_n_configSEXP, SEXP model_type_SEXP, SEXP latency_periodSEXP, SEXP establishment_probabilitySEXP, SEXP dispersal_percentageSEXP, SEXP overpopulation_configSEXP, SEXP network_configSEXP, SEXP network_data_configSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::traits::input_parameter< int >::type random_seed(random_seedSEXP); @@ -38,6 +38,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< std::vector >::type movements_dates(movements_datesSEXP); Rcpp::traits::input_parameter< std::vector >::type temperature(temperatureSEXP); Rcpp::traits::input_parameter< std::vector >::type weather_coefficient(weather_coefficientSEXP); + Rcpp::traits::input_parameter< List >::type bbox(bboxSEXP); Rcpp::traits::input_parameter< List >::type res(resSEXP); Rcpp::traits::input_parameter< List >::type rows_cols(rows_colsSEXP); Rcpp::traits::input_parameter< double >::type reproductive_rate(reproductive_rateSEXP); @@ -59,24 +60,23 @@ BEGIN_RCPP Rcpp::traits::input_parameter< double >::type natural_kappa(natural_kappaSEXP); Rcpp::traits::input_parameter< std::string >::type anthropogenic_dir(anthropogenic_dirSEXP); Rcpp::traits::input_parameter< double >::type anthropogenic_kappa(anthropogenic_kappaSEXP); - Rcpp::traits::input_parameter< int >::type output_frequency_n(output_frequency_nSEXP); - Rcpp::traits::input_parameter< int >::type quarantine_frequency_n(quarantine_frequency_nSEXP); - Rcpp::traits::input_parameter< int >::type spreadrate_frequency_n(spreadrate_frequency_nSEXP); - Rcpp::traits::input_parameter< int >::type mortality_frequency_n(mortality_frequency_nSEXP); + Rcpp::traits::input_parameter< Nullable >::type frequencies_n_config(frequencies_n_configSEXP); Rcpp::traits::input_parameter< std::string >::type model_type_(model_type_SEXP); Rcpp::traits::input_parameter< int >::type latency_period(latency_periodSEXP); Rcpp::traits::input_parameter< double >::type establishment_probability(establishment_probabilitySEXP); Rcpp::traits::input_parameter< double >::type dispersal_percentage(dispersal_percentageSEXP); Rcpp::traits::input_parameter< Nullable >::type overpopulation_config(overpopulation_configSEXP); - rcpp_result_gen = Rcpp::wrap(pops_model_cpp(random_seed, lethal_temperature, lethal_temperature_month, infected, total_exposed, exposed, susceptible, total_populations, total_hosts, mortality_tracker, mortality, quarantine_areas, treatment_maps, treatment_dates, pesticide_duration, resistant, movements, movements_dates, temperature, weather_coefficient, res, rows_cols, reproductive_rate, spatial_indices, season_month_start_end, frequency_config, bool_config, mortality_rate, mortality_time_lag, start_date, end_date, treatment_method, natural_kernel_type, anthropogenic_kernel_type, percent_natural_dispersal, natural_distance_scale, anthropogenic_distance_scale, natural_dir, natural_kappa, anthropogenic_dir, anthropogenic_kappa, output_frequency_n, quarantine_frequency_n, spreadrate_frequency_n, mortality_frequency_n, model_type_, latency_period, establishment_probability, dispersal_percentage, overpopulation_config)); + Rcpp::traits::input_parameter< Nullable >::type network_config(network_configSEXP); + Rcpp::traits::input_parameter< Nullable >::type network_data_config(network_data_configSEXP); + rcpp_result_gen = Rcpp::wrap(pops_model_cpp(random_seed, lethal_temperature, lethal_temperature_month, infected, total_exposed, exposed, susceptible, total_populations, total_hosts, mortality_tracker, mortality, quarantine_areas, treatment_maps, treatment_dates, pesticide_duration, resistant, movements, movements_dates, temperature, weather_coefficient, bbox, res, rows_cols, reproductive_rate, spatial_indices, season_month_start_end, frequency_config, bool_config, mortality_rate, mortality_time_lag, start_date, end_date, treatment_method, natural_kernel_type, anthropogenic_kernel_type, percent_natural_dispersal, natural_distance_scale, anthropogenic_distance_scale, natural_dir, natural_kappa, anthropogenic_dir, anthropogenic_kappa, frequencies_n_config, model_type_, latency_period, establishment_probability, dispersal_percentage, overpopulation_config, network_config, network_data_config)); return rcpp_result_gen; END_RCPP_RETURN_ERROR } -RcppExport SEXP _PoPS_pops_model_cpp(SEXP random_seedSEXP, SEXP lethal_temperatureSEXP, SEXP lethal_temperature_monthSEXP, SEXP infectedSEXP, SEXP total_exposedSEXP, SEXP exposedSEXP, SEXP susceptibleSEXP, SEXP total_populationsSEXP, SEXP total_hostsSEXP, SEXP mortality_trackerSEXP, SEXP mortalitySEXP, SEXP quarantine_areasSEXP, SEXP treatment_mapsSEXP, SEXP treatment_datesSEXP, SEXP pesticide_durationSEXP, SEXP resistantSEXP, SEXP movementsSEXP, SEXP movements_datesSEXP, SEXP temperatureSEXP, SEXP weather_coefficientSEXP, SEXP resSEXP, SEXP rows_colsSEXP, SEXP reproductive_rateSEXP, SEXP spatial_indicesSEXP, SEXP season_month_start_endSEXP, SEXP frequency_configSEXP, SEXP bool_configSEXP, SEXP mortality_rateSEXP, SEXP mortality_time_lagSEXP, SEXP start_dateSEXP, SEXP end_dateSEXP, SEXP treatment_methodSEXP, SEXP natural_kernel_typeSEXP, SEXP anthropogenic_kernel_typeSEXP, SEXP percent_natural_dispersalSEXP, SEXP natural_distance_scaleSEXP, SEXP anthropogenic_distance_scaleSEXP, SEXP natural_dirSEXP, SEXP natural_kappaSEXP, SEXP anthropogenic_dirSEXP, SEXP anthropogenic_kappaSEXP, SEXP output_frequency_nSEXP, SEXP quarantine_frequency_nSEXP, SEXP spreadrate_frequency_nSEXP, SEXP mortality_frequency_nSEXP, SEXP model_type_SEXP, SEXP latency_periodSEXP, SEXP establishment_probabilitySEXP, SEXP dispersal_percentageSEXP, SEXP overpopulation_configSEXP) { +RcppExport SEXP _PoPS_pops_model_cpp(SEXP random_seedSEXP, SEXP lethal_temperatureSEXP, SEXP lethal_temperature_monthSEXP, SEXP infectedSEXP, SEXP total_exposedSEXP, SEXP exposedSEXP, SEXP susceptibleSEXP, SEXP total_populationsSEXP, SEXP total_hostsSEXP, SEXP mortality_trackerSEXP, SEXP mortalitySEXP, SEXP quarantine_areasSEXP, SEXP treatment_mapsSEXP, SEXP treatment_datesSEXP, SEXP pesticide_durationSEXP, SEXP resistantSEXP, SEXP movementsSEXP, SEXP movements_datesSEXP, SEXP temperatureSEXP, SEXP weather_coefficientSEXP, SEXP bboxSEXP, SEXP resSEXP, SEXP rows_colsSEXP, SEXP reproductive_rateSEXP, SEXP spatial_indicesSEXP, SEXP season_month_start_endSEXP, SEXP frequency_configSEXP, SEXP bool_configSEXP, SEXP mortality_rateSEXP, SEXP mortality_time_lagSEXP, SEXP start_dateSEXP, SEXP end_dateSEXP, SEXP treatment_methodSEXP, SEXP natural_kernel_typeSEXP, SEXP anthropogenic_kernel_typeSEXP, SEXP percent_natural_dispersalSEXP, SEXP natural_distance_scaleSEXP, SEXP anthropogenic_distance_scaleSEXP, SEXP natural_dirSEXP, SEXP natural_kappaSEXP, SEXP anthropogenic_dirSEXP, SEXP anthropogenic_kappaSEXP, SEXP frequencies_n_configSEXP, SEXP model_type_SEXP, SEXP latency_periodSEXP, SEXP establishment_probabilitySEXP, SEXP dispersal_percentageSEXP, SEXP overpopulation_configSEXP, SEXP network_configSEXP, SEXP network_data_configSEXP) { SEXP rcpp_result_gen; { Rcpp::RNGScope rcpp_rngScope_gen; - rcpp_result_gen = PROTECT(_PoPS_pops_model_cpp_try(random_seedSEXP, lethal_temperatureSEXP, lethal_temperature_monthSEXP, infectedSEXP, total_exposedSEXP, exposedSEXP, susceptibleSEXP, total_populationsSEXP, total_hostsSEXP, mortality_trackerSEXP, mortalitySEXP, quarantine_areasSEXP, treatment_mapsSEXP, treatment_datesSEXP, pesticide_durationSEXP, resistantSEXP, movementsSEXP, movements_datesSEXP, temperatureSEXP, weather_coefficientSEXP, resSEXP, rows_colsSEXP, reproductive_rateSEXP, spatial_indicesSEXP, season_month_start_endSEXP, frequency_configSEXP, bool_configSEXP, mortality_rateSEXP, mortality_time_lagSEXP, start_dateSEXP, end_dateSEXP, treatment_methodSEXP, natural_kernel_typeSEXP, anthropogenic_kernel_typeSEXP, percent_natural_dispersalSEXP, natural_distance_scaleSEXP, anthropogenic_distance_scaleSEXP, natural_dirSEXP, natural_kappaSEXP, anthropogenic_dirSEXP, anthropogenic_kappaSEXP, output_frequency_nSEXP, quarantine_frequency_nSEXP, spreadrate_frequency_nSEXP, mortality_frequency_nSEXP, model_type_SEXP, latency_periodSEXP, establishment_probabilitySEXP, dispersal_percentageSEXP, overpopulation_configSEXP)); + rcpp_result_gen = PROTECT(_PoPS_pops_model_cpp_try(random_seedSEXP, lethal_temperatureSEXP, lethal_temperature_monthSEXP, infectedSEXP, total_exposedSEXP, exposedSEXP, susceptibleSEXP, total_populationsSEXP, total_hostsSEXP, mortality_trackerSEXP, mortalitySEXP, quarantine_areasSEXP, treatment_mapsSEXP, treatment_datesSEXP, pesticide_durationSEXP, resistantSEXP, movementsSEXP, movements_datesSEXP, temperatureSEXP, weather_coefficientSEXP, bboxSEXP, resSEXP, rows_colsSEXP, reproductive_rateSEXP, spatial_indicesSEXP, season_month_start_endSEXP, frequency_configSEXP, bool_configSEXP, mortality_rateSEXP, mortality_time_lagSEXP, start_dateSEXP, end_dateSEXP, treatment_methodSEXP, natural_kernel_typeSEXP, anthropogenic_kernel_typeSEXP, percent_natural_dispersalSEXP, natural_distance_scaleSEXP, anthropogenic_distance_scaleSEXP, natural_dirSEXP, natural_kappaSEXP, anthropogenic_dirSEXP, anthropogenic_kappaSEXP, frequencies_n_configSEXP, model_type_SEXP, latency_periodSEXP, establishment_probabilitySEXP, dispersal_percentageSEXP, overpopulation_configSEXP, network_configSEXP, network_data_configSEXP)); } Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); if (rcpp_isInterrupt_gen) { @@ -101,7 +101,7 @@ RcppExport SEXP _PoPS_pops_model_cpp(SEXP random_seedSEXP, SEXP lethal_temperatu static int _PoPS_RcppExport_validate(const char* sig) { static std::set signatures; if (signatures.empty()) { - signatures.insert("List(*pops_model_cpp)(int,double,int,IntegerMatrix,IntegerMatrix,std::vector,IntegerMatrix,IntegerMatrix,IntegerMatrix,std::vector,IntegerMatrix,IntegerMatrix,std::vector,std::vector,std::vector,IntegerMatrix,std::vector>,std::vector,std::vector,std::vector,List,List,double,std::vector>,List,List,List,double,int,std::string,std::string,std::string,std::string,std::string,double,double,double,std::string,double,std::string,double,int,int,int,int,std::string,int,double,double,Nullable)"); + signatures.insert("List(*pops_model_cpp)(int,double,int,IntegerMatrix,IntegerMatrix,std::vector,IntegerMatrix,IntegerMatrix,IntegerMatrix,std::vector,IntegerMatrix,IntegerMatrix,std::vector,std::vector,std::vector,IntegerMatrix,std::vector>,std::vector,std::vector,std::vector,List,List,List,double,std::vector>,List,List,List,double,int,std::string,std::string,std::string,std::string,std::string,double,double,double,std::string,double,std::string,double,Nullable,std::string,int,double,double,Nullable,Nullable,Nullable)"); } return signatures.find(sig) != signatures.end(); } diff --git a/src/pops.cpp b/src/pops.cpp index b921726d..f31df6e2 100644 --- a/src/pops.cpp +++ b/src/pops.cpp @@ -59,6 +59,7 @@ List pops_model_cpp( std::vector movements_dates, std::vector temperature, std::vector weather_coefficient, + List bbox, List res, List rows_cols, double reproductive_rate, @@ -80,15 +81,14 @@ List pops_model_cpp( double natural_kappa = 0, std::string anthropogenic_dir = "NONE", double anthropogenic_kappa = 0, - int output_frequency_n = 1, - int quarantine_frequency_n = 1, - int spreadrate_frequency_n = 1, - int mortality_frequency_n = 1, + Nullable frequencies_n_config = R_NilValue, std::string model_type_ = "SI", int latency_period = 0, double establishment_probability = 0, double dispersal_percentage = 0.99, - Nullable overpopulation_config = R_NilValue) + Nullable overpopulation_config = R_NilValue, + Nullable network_config = R_NilValue, + Nullable network_data_config = R_NilValue) { Config config; config.random_seed = random_seed; @@ -137,15 +137,20 @@ List pops_model_cpp( // movement_schedule used later config.dispersal_percentage = dispersal_percentage; + + if (frequencies_n_config.isNotNull()) { + List freq_n_config(frequencies_n_config); + config.output_frequency_n = freq_n_config["output_frequency_n"]; + config.quarantine_frequency_n = freq_n_config["quarantine_frequency_n"]; + config.spreadrate_frequency_n = freq_n_config["spreadrate_frequency_n"]; + config.mortality_frequency_n = freq_n_config["mortality_frequency_n"]; + } + config.output_frequency = output_frequency; - config.output_frequency_n = output_frequency_n; config.quarantine_frequency = quarantine_frequency; - config.quarantine_frequency_n = quarantine_frequency_n; config.use_quarantine = bool_config["use_quarantine"]; config.spreadrate_frequency = spreadrate_frequency; - config.spreadrate_frequency_n = spreadrate_frequency_n; config.mortality_frequency = mortality_frequency; - config.mortality_frequency_n = mortality_frequency_n; config.use_spreadrates = bool_config["use_spreadrates"]; config.use_overpopulation_movements = bool_config["use_overpopulation_movements"]; if (config.use_overpopulation_movements && overpopulation_config.isNotNull()) { @@ -169,6 +174,7 @@ List pops_model_cpp( std::vector> spread_rates_vector; std::tuple spread_rates; IntegerMatrix total_dispersers(config.rows, config.cols); + IntegerMatrix established_dispersers(config.rows, config.cols); int num_infected; std::vector number_infected; @@ -244,6 +250,23 @@ List pops_model_cpp( Direction escape_direction; std::vector escape_directions; + std::unique_ptr> network{nullptr}; + if (network_config.isNotNull() && network_data_config.isNotNull()) { + // The best place for bbox handling would be with rows, cols, and + // resolution, but since it is required only for network, it is here. + config.bbox.north = bbox["north"]; + config.bbox.south = bbox["south"]; + config.bbox.east = bbox["east"]; + config.bbox.west = bbox["west"]; + List net_config(network_config); + config.network_min_distance = net_config["network_min_distance"]; + config.network_max_distance = net_config["network_max_distance"]; + network.reset(new Network(config.bbox, config.ew_res, config.ns_res)); + List net_data_config(network_data_config); + std::ifstream network_stream{Rcpp::as(net_data_config["network_filename"])}; + network->load(network_stream); + } + ModelType mt = model_type_from_string(config.model_type); Simulation simulation( config.random_seed, config.rows, config.cols, mt, config.latency_period_steps); @@ -260,6 +283,7 @@ List pops_model_cpp( total_populations, total_hosts, dispersers, + established_dispersers, total_exposed, exposed, mortality_tracker, @@ -273,9 +297,10 @@ List pops_model_cpp( quarantine, quarantine_areas, movements, + network ? *network : Network::null_network(), spatial_indices); - // keeps track of cumulative dispers or propagules from a site. + // keeps track of cumulative dispersers or propagules from a site. if (config.spread_schedule()[current_index]) { total_dispersers += dispersers; } diff --git a/tests/testthat/test-pops.r b/tests/testthat/test-pops.r index 8c8cc04e..82bba0f3 100644 --- a/tests/testthat/test-pops.r +++ b/tests/testthat/test-pops.r @@ -2622,3 +2622,41 @@ test_that( expect_gte(data$infected[[1]][[3]], test_mat[[3]]) expect_gte(data$infected[[1]][[4]], test_mat[[4]]) }) + + +test_that( + "Network dispersal works as expected", { + infected_file <- + system.file("extdata", "simple20x20", "initial_infection.tif", package = "PoPS") + host_file <- system.file("extdata", "simple20x20", "all_plants.tif", package = "PoPS") + start_date <- "2008-01-01" + end_date <- "2008-03-31" + parameter_means <- c(2, 21, 0.5, 500, 0, 0) + parameter_cov_matrix <- matrix(0, nrow = 6, ncol = 6) + network_filename <- system.file("extdata", "simple20x20", "segments.csv", package = "PoPS") + network_min_distance <- 10 + network_max_distance <- 100 + anthropogenic_kernel_type <- "network" + # anthropogenic_kernel_type <- "cauchy" + # node_filename <- "" + # segment_filename <- "" + + data <- + pops(infected_file = infected_file, + host_file = host_file, + total_populations_file = host_file, + parameter_means = parameter_means, + parameter_cov_matrix = parameter_cov_matrix, + start_date = start_date, + end_date = end_date, + anthropogenic_kernel_type = anthropogenic_kernel_type, + network_min_distance = network_min_distance, + network_max_distance = network_max_distance, + network_filename = network_filename) + + test_mat <- terra::as.matrix(terra::rast(infected_file), wide = TRUE) + expect_gte(data$infected[[1]][[1]], test_mat[[1]]) + expect_gte(data$infected[[1]][[2]], test_mat[[2]]) + expect_gte(data$infected[[1]][[3]], test_mat[[3]]) + expect_gte(data$infected[[1]][[4]], test_mat[[4]]) + })