diff --git a/DESCRIPTION b/DESCRIPTION index 1bdb0678..d599f8c8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: epiworldR Type: Package Title: Fast Agent-Based Epi Models -Version: 0.4-3 +Version: 0.5-0 Authors@R: c( person(given="George", family="Vega Yon", role=c("aut","cre"), email="g.vegayon@gmail.com", comment = c(ORCID = "0000-0002-3171-0844")), diff --git a/NAMESPACE b/NAMESPACE index 82bf0308..aac1efb8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -27,10 +27,8 @@ S3method(get_n_viruses,epiworld_model) S3method(get_name,epiworld_model) S3method(get_ndays,epiworld_model) S3method(get_param,epiworld_model) -S3method(get_params_mean,epiworld_lfmcmc) S3method(get_reproductive_number,epiworld_model) S3method(get_states,epiworld_model) -S3method(get_stats_mean,epiworld_lfmcmc) S3method(get_today_total,epiworld_model) S3method(get_transition_probability,epiworld_model) S3method(get_transmissions,epiworld_diffnet) @@ -77,17 +75,9 @@ S3method(queuing_on,epiworld_model) S3method(queuing_on,epiworld_seirconn) S3method(queuing_on,epiworld_sirconn) S3method(run,epiworld_model) -S3method(run_lfmcmc,epiworld_lfmcmc) S3method(run_multiple,epiworld_model) -S3method(set_kernel_fun,epiworld_lfmcmc) S3method(set_name,epiworld_model) -S3method(set_observed_data,epiworld_lfmcmc) -S3method(set_par_names,epiworld_lfmcmc) S3method(set_param,epiworld_model) -S3method(set_proposal_fun,epiworld_lfmcmc) -S3method(set_simulation_fun,epiworld_lfmcmc) -S3method(set_stats_names,epiworld_lfmcmc) -S3method(set_summary_fun,epiworld_lfmcmc) S3method(size,epiworld_model) S3method(summary,epiworld_model) S3method(today,epiworld_model) @@ -131,6 +121,8 @@ export(distribute_virus_set) export(entity) export(entity_add_agent) export(entity_get_agents) +export(get_accepted_params) +export(get_accepted_stats) export(get_agents) export(get_agents_data_ncols) export(get_agents_states) @@ -143,7 +135,10 @@ export(get_hist_tool) export(get_hist_total) export(get_hist_transition_matrix) export(get_hist_virus) +export(get_n_parameters) export(get_n_replicates) +export(get_n_samples) +export(get_n_statistics) export(get_n_tools) export(get_n_viruses) export(get_name) @@ -156,6 +151,7 @@ export(get_params_mean) export(get_reproductive_number) export(get_state) export(get_states) +export(get_statistics_hist) export(get_stats_mean) export(get_today_total) export(get_tool) diff --git a/R/LFMCMC.R b/R/LFMCMC.R index 6b8da1fe..f0ada406 100644 --- a/R/LFMCMC.R +++ b/R/LFMCMC.R @@ -1,3 +1,13 @@ +stopifnot_lfmcmc <- function(x) { + # Catching the value of x + nam <- match.call()$x + + if (!inherits(x, "epiworld_lfmcmc")) + stop(nam, " must be an object of class epiworld_lfmcmc") + +} + + #' Likelihood-Free Markhov Chain Monte Carlo (LFMCMC) #' #' @aliases epiworld_lfmcmc @@ -102,14 +112,8 @@ run_lfmcmc <- function( lfmcmc, params_init_, n_samples_, epsilon_, seed = NULL ) { - UseMethod("run_lfmcmc") -} -#' @export -run_lfmcmc.epiworld_lfmcmc <- function( - lfmcmc, params_init_, n_samples_, epsilon_, - seed = NULL - ) { + stopifnot_lfmcmc(lfmcmc) if (length(seed)) set.seed(seed) @@ -131,32 +135,29 @@ run_lfmcmc.epiworld_lfmcmc <- function( #' @param observed_data_ Observed data, treated as double #' @returns The lfmcmc model with the observed data added #' @export -set_observed_data <- function( - lfmcmc, observed_data_ - ) { - UseMethod("set_observed_data") -} +set_observed_data <- function(lfmcmc, observed_data_) { + + stopifnot_lfmcmc(lfmcmc) -#' @export -set_observed_data.epiworld_lfmcmc <- function(lfmcmc, observed_data_) { set_observed_data_cpp( lfmcmc, as.double(observed_data_) ) + invisible(lfmcmc) } #' @rdname LFMCMC #' @param lfmcmc LFMCMC model #' @param fun The LFMCMC proposal function -#' @returns The lfmcmc model with the proposal function added #' @export -set_proposal_fun <- function(lfmcmc, fun) UseMethod("set_proposal_fun") +#' @returns The lfmcmc model with the proposal function added +set_proposal_fun <- function(lfmcmc, fun) { -#' @export -set_proposal_fun.epiworld_lfmcmc <- function(lfmcmc, fun) { + stopifnot_lfmcmc(lfmcmc) set_proposal_fun_cpp(lfmcmc, fun) invisible(lfmcmc) + } #' @rdname LFMCMC @@ -164,8 +165,11 @@ set_proposal_fun.epiworld_lfmcmc <- function(lfmcmc, fun) { #' @returns The LFMCMC model with proposal function set to norm reflective #' @export use_proposal_norm_reflective <- function(lfmcmc) { + + stopifnot_lfmcmc(lfmcmc) use_proposal_norm_reflective_cpp(lfmcmc) invisible(lfmcmc) + } #' @rdname LFMCMC @@ -173,12 +177,12 @@ use_proposal_norm_reflective <- function(lfmcmc) { #' @param fun The LFMCMC simulation function #' @returns The lfmcmc model with the simulation function added #' @export -set_simulation_fun <- function(lfmcmc, fun) UseMethod("set_simulation_fun") +set_simulation_fun <- function(lfmcmc, fun) { -#' @export -set_simulation_fun.epiworld_lfmcmc <- function(lfmcmc, fun) { + stopifnot_lfmcmc(lfmcmc) set_simulation_fun_cpp(lfmcmc, fun) invisible(lfmcmc) + } #' @rdname LFMCMC @@ -186,12 +190,12 @@ set_simulation_fun.epiworld_lfmcmc <- function(lfmcmc, fun) { #' @param fun The LFMCMC sumamry function #' @returns The lfmcmc model with the summary function added #' @export -set_summary_fun <- function(lfmcmc, fun) UseMethod("set_summary_fun") +set_summary_fun <- function(lfmcmc, fun) { -#' @export -set_summary_fun.epiworld_lfmcmc <- function(lfmcmc, fun) { + stopifnot_lfmcmc(lfmcmc) set_summary_fun_cpp(lfmcmc, fun) invisible(lfmcmc) + } #' @rdname LFMCMC @@ -199,12 +203,12 @@ set_summary_fun.epiworld_lfmcmc <- function(lfmcmc, fun) { #' @param fun The LFMCMC kernel function #' @returns The lfmcmc model with the kernel function added #' @export -set_kernel_fun <- function(lfmcmc, fun) UseMethod("set_kernel_fun") +set_kernel_fun <- function(lfmcmc, fun) { -#' @export -set_kernel_fun.epiworld_lfmcmc <- function(lfmcmc, fun) { + stopifnot_lfmcmc(lfmcmc) set_kernel_fun_cpp(lfmcmc, fun) invisible(lfmcmc) + } #' @rdname LFMCMC @@ -212,8 +216,11 @@ set_kernel_fun.epiworld_lfmcmc <- function(lfmcmc, fun) { #' @returns The LFMCMC model with kernel function set to gaussian #' @export use_kernel_fun_gaussian <- function(lfmcmc) { + + stopifnot_lfmcmc(lfmcmc) use_kernel_fun_gaussian_cpp(lfmcmc) invisible(lfmcmc) + } #' @rdname LFMCMC @@ -221,12 +228,12 @@ use_kernel_fun_gaussian <- function(lfmcmc) { #' @param names The model parameter names #' @returns The lfmcmc model with the parameter names added #' @export -set_par_names <- function(lfmcmc, names) UseMethod("set_par_names") +set_par_names <- function(lfmcmc, names) { -#' @export -set_par_names.epiworld_lfmcmc <- function(lfmcmc, names) { + stopifnot_lfmcmc(lfmcmc) set_par_names_cpp(lfmcmc, names) invisible(lfmcmc) + } #' @rdname LFMCMC @@ -234,42 +241,144 @@ set_par_names.epiworld_lfmcmc <- function(lfmcmc, names) { #' @param names The model stats names #' @returns The lfmcmc model with the stats names added #' @export -set_stats_names <- function(lfmcmc, names) UseMethod("set_stats_names") +set_stats_names <- function(lfmcmc, names) { -#' @export -set_stats_names.epiworld_lfmcmc <- function(lfmcmc, names) { + stopifnot_lfmcmc(lfmcmc) set_stats_names_cpp(lfmcmc, names) invisible(lfmcmc) + } #' @rdname LFMCMC #' @param lfmcmc LFMCMC model #' @returns The param means for the given lfmcmc model #' @export -get_params_mean <- function(lfmcmc) UseMethod("get_params_mean") +get_params_mean <- function(lfmcmc) { -#' @export -get_params_mean.epiworld_lfmcmc <- function(lfmcmc) { + stopifnot_lfmcmc(lfmcmc) get_params_mean_cpp(lfmcmc) + } #' @rdname LFMCMC #' @param lfmcmc LFMCMC model #' @returns The stats means for the given lfmcmc model #' @export -get_stats_mean <- function(lfmcmc) UseMethod("get_stats_mean") +get_stats_mean <- function(lfmcmc) { -#' @export -get_stats_mean.epiworld_lfmcmc <- function(lfmcmc) { + stopifnot_lfmcmc(lfmcmc) get_stats_mean_cpp(lfmcmc) + } #' @rdname LFMCMC #' @param x LFMCMC model to print #' @param ... Ignored +#' @param burnin Integer. Number of samples to discard as burnin before computing the summary. #' @returns The lfmcmc model #' @export -print.epiworld_lfmcmc <- function(x, ...) { - print_lfmcmc_cpp(x) +print.epiworld_lfmcmc <- function(x, burnin = 0, ...) { + + if (!is.numeric(burnin)) + stop("The 'burnin' argument must be an integer.") + + if (burnin < 0) + stop("The 'burnin' argument must be a non-negative integer.") + + print_lfmcmc_cpp(x, burnin = burnin) invisible(x) + +} + +#' @rdname LFMCMC +#' @export +#' @returns +#' - The function `get_accepted_params` returns a matrix of accepted +#' parameters for the given LFMCMC model. with the number of rows equal to the +#' number of samples and the number of columns equal to the number of +#' parameters. +get_accepted_params <- function(lfmcmc) { + + stopifnot_lfmcmc(lfmcmc) + a_params <- get_accepted_params_cpp(lfmcmc) + n_params <- get_n_parameters(lfmcmc) + + matrix( + a_params, + ncol = n_params, + byrow = TRUE + ) + +} + + +#' @rdname LFMCMC +#' @export +#' @returns +#' - The function `get_accepted_stats` returns a matrix of accepted statistics +#' for the given LFMCMC model. with the number of rows equal to the number of +#' samples and the number of columns equal to the number of statistics. +get_accepted_stats <- function(lfmcmc) { + + stopifnot_lfmcmc(lfmcmc) + a_stats <- get_accepted_stats_cpp(lfmcmc) + n_stats <- get_n_statistics(lfmcmc) + + matrix( + a_stats, + ncol = n_stats, + byrow = TRUE + ) + +} + +#' @export +#' @rdname LFMCMC +#' @returns +#' - The function `get_statistics_hist` returns a matrix of statistics +#' for the given LFMCMC model. with the number of rows equal to the number of +#' samples and the number of columns equal to the number of statistics. +get_statistics_hist <- function(lfmcmc) { + + stopifnot_lfmcmc(lfmcmc) + stats <- get_statistics_hist_cpp(lfmcmc) + n_stats <- get_n_statistics(lfmcmc) + + matrix( + stats, + ncol = n_stats, + byrow = TRUE + ) + +} + +#' @export +#' @rdname LFMCMC +#' @returns +#' - The functions `get_n_parameters`, `get_n_statistics`, and `get_n_samples` +#' return the number of parameters, statistics, and samples for the given +#' LFMCMC model, respectively. +get_n_parameters <- function(lfmcmc) { + + stopifnot_lfmcmc(lfmcmc) + get_n_parameters_cpp(lfmcmc) + +} + +#' @export +#' @rdname LFMCMC +get_n_statistics <- function(lfmcmc) { + + stopifnot_lfmcmc(lfmcmc) + get_n_statistics_cpp(lfmcmc) + +} + +#' @export +#' @rdname LFMCMC +get_n_samples <- function(lfmcmc) { + + stopifnot_lfmcmc(lfmcmc) + get_n_samples_cpp(lfmcmc) + } diff --git a/R/cpp11.R b/R/cpp11.R index 7b58bc8f..4bb2accf 100644 --- a/R/cpp11.R +++ b/R/cpp11.R @@ -272,8 +272,32 @@ get_stats_mean_cpp <- function(lfmcmc) { .Call(`_epiworldR_get_stats_mean_cpp`, lfmcmc) } -print_lfmcmc_cpp <- function(lfmcmc) { - .Call(`_epiworldR_print_lfmcmc_cpp`, lfmcmc) +print_lfmcmc_cpp <- function(lfmcmc, burnin) { + .Call(`_epiworldR_print_lfmcmc_cpp`, lfmcmc, burnin) +} + +get_statistics_hist_cpp <- function(lfmcmc) { + .Call(`_epiworldR_get_statistics_hist_cpp`, lfmcmc) +} + +get_accepted_params_cpp <- function(lfmcmc) { + .Call(`_epiworldR_get_accepted_params_cpp`, lfmcmc) +} + +get_accepted_stats_cpp <- function(lfmcmc) { + .Call(`_epiworldR_get_accepted_stats_cpp`, lfmcmc) +} + +get_n_samples_cpp <- function(lfmcmc) { + .Call(`_epiworldR_get_n_samples_cpp`, lfmcmc) +} + +get_n_statistics_cpp <- function(lfmcmc) { + .Call(`_epiworldR_get_n_statistics_cpp`, lfmcmc) +} + +get_n_parameters_cpp <- function(lfmcmc) { + .Call(`_epiworldR_get_n_parameters_cpp`, lfmcmc) } print_cpp <- function(m, lite) { diff --git a/inst/include/epiworld/agent-meat-virus-sampling.hpp b/inst/include/epiworld/agent-meat-virus-sampling.hpp index d7289937..bfe668ee 100644 --- a/inst/include/epiworld/agent-meat-virus-sampling.hpp +++ b/inst/include/epiworld/agent-meat-virus-sampling.hpp @@ -20,7 +20,7 @@ namespace sampler { * @return Virus* of the selected virus. If none selected (or none * available,) returns a nullptr; */ -template +template inline std::function*,Model*)> make_update_susceptible( std::vector< epiworld_fast_uint > exclude = {} ) @@ -179,7 +179,7 @@ inline std::function*,Model*)> make_update_susceptible( * @return Virus* of the selected virus. If none selected (or none * available,) returns a nullptr; */ -template +template inline std::function*(Agent*,Model*)> make_sample_virus_neighbors( std::vector< epiworld_fast_uint > exclude = {} ) @@ -347,7 +347,7 @@ inline std::function*(Agent*,Model*)> make_sample_virus_ * @return Virus* of the selected virus. If none selected (or none * available,) returns a nullptr; */ -template +template inline Virus * sample_virus_single(Agent * p, Model * m) { diff --git a/inst/include/epiworld/epiworld.hpp b/inst/include/epiworld/epiworld.hpp index 6a060f7e..12801802 100644 --- a/inst/include/epiworld/epiworld.hpp +++ b/inst/include/epiworld/epiworld.hpp @@ -18,8 +18,8 @@ /* Versioning */ #define EPIWORLD_VERSION_MAJOR 0 -#define EPIWORLD_VERSION_MINOR 4 -#define EPIWORLD_VERSION_PATCH 3 +#define EPIWORLD_VERSION_MINOR 5 +#define EPIWORLD_VERSION_PATCH 0 static const int epiworld_version_major = EPIWORLD_VERSION_MAJOR; static const int epiworld_version_minor = EPIWORLD_VERSION_MINOR; @@ -88,4 +88,4 @@ namespace epiworld { } -#endif \ No newline at end of file +#endif diff --git a/inst/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp b/inst/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp index 7002aa24..ab86172b 100755 --- a/inst/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp +++ b/inst/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp @@ -177,7 +177,7 @@ class LFMCMC { epiworld_double * last_elapsed, std::string * unit_abbr, bool print - ); + ) const; void chrono_start(); void chrono_end(); @@ -235,13 +235,17 @@ class LFMCMC { const std::vector< epiworld_double > & get_drawn_prob() {return drawn_prob;}; std::vector< TData > * get_sampled_data() {return sampled_data;}; + const std::vector< epiworld_double > & get_accepted_params() {return accepted_params;}; + const std::vector< epiworld_double > & get_accepted_stats() {return accepted_stats;}; + + void set_par_names(std::vector< std::string > names); void set_stats_names(std::vector< std::string > names); std::vector< epiworld_double > get_params_mean(); std::vector< epiworld_double > get_stats_mean(); - void print() ; + void print(size_t burnin = 0u) const; }; diff --git a/inst/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp b/inst/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp index 363fb766..b2860a08 100755 --- a/inst/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp +++ b/inst/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp @@ -2,50 +2,71 @@ #define LFMCMC_MEAT_PRINT_HPP template -inline void LFMCMC::print() +inline void LFMCMC::print(size_t burnin) const { + + // For each statistic or parameter in the model, we print three values: + // - mean, the 2.5% quantile, and the 97.5% quantile std::vector< epiworld_double > summ_params(n_parameters * 3, 0.0); std::vector< epiworld_double > summ_stats(n_statistics * 3, 0.0); + // Compute the number of samples to use based on burnin rate + size_t n_samples_print = n_samples; + if (burnin > 0) + { + if (burnin >= n_samples) + throw std::length_error( + "The burnin is greater than or equal to the number of samples." + ); + + n_samples_print = n_samples - burnin; + + } + + epiworld_double n_samples_dbl = static_cast< epiworld_double >( + n_samples_print + ); + + // Compute parameter summary values for (size_t k = 0u; k < n_parameters; ++k) { // Retrieving the relevant parameter - std::vector< epiworld_double > par_i(n_samples); - for (size_t i = 0u; i < n_samples; ++i) + std::vector< epiworld_double > par_i(n_samples_print); + for (size_t i = burnin; i < n_samples; ++i) { - par_i[i] = accepted_params[i * n_parameters + k]; - summ_params[k * 3] += par_i[i]/n_samples; + par_i[i-burnin] = accepted_params[i * n_parameters + k]; + summ_params[k * 3] += par_i[i-burnin]/n_samples_dbl; } // Computing the 95% Credible interval std::sort(par_i.begin(), par_i.end()); summ_params[k * 3 + 1u] = - par_i[std::floor(.025 * static_cast(n_samples))]; + par_i[std::floor(.025 * n_samples_dbl)]; summ_params[k * 3 + 2u] = - par_i[std::floor(.975 * static_cast(n_samples))]; + par_i[std::floor(.975 * n_samples_dbl)]; } + // Compute statistics summary values for (size_t k = 0u; k < n_statistics; ++k) { // Retrieving the relevant parameter - std::vector< epiworld_double > stat_k(n_samples); - for (size_t i = 0u; i < n_samples; ++i) + std::vector< epiworld_double > stat_k(n_samples_print); + for (size_t i = burnin; i < n_samples; ++i) { - stat_k[i] = accepted_stats[i * n_statistics + k]; - summ_stats[k * 3] += stat_k[i]/n_samples; + stat_k[i-burnin] = accepted_stats[i * n_statistics + k]; + summ_stats[k * 3] += stat_k[i-burnin]/n_samples_dbl; } // Computing the 95% Credible interval std::sort(stat_k.begin(), stat_k.end()); - summ_stats[k * 3 + 1u] = - stat_k[std::floor(.025 * static_cast(n_samples))]; + stat_k[std::floor(.025 * n_samples_dbl)]; summ_stats[k * 3 + 2u] = - stat_k[std::floor(.975 * static_cast(n_samples))]; + stat_k[std::floor(.975 * n_samples_dbl)]; } diff --git a/inst/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp b/inst/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp index c7565121..c515d50c 100755 --- a/inst/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp +++ b/inst/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp @@ -249,14 +249,16 @@ inline void LFMCMC::run( std::vector< epiworld_double > proposed_stats_i; summary_fun(proposed_stats_i, data_i, this); - accepted_params_prob[0u] = kernel_fun(proposed_stats_i, observed_stats, epsilon, this); + accepted_params_prob[0u] = kernel_fun( + proposed_stats_i, observed_stats, epsilon, this + ); // Recording statistics for (size_t i = 0u; i < n_statistics; ++i) sampled_stats[i] = proposed_stats_i[i]; - for (size_t k = 0u; k < n_statistics; ++k) - accepted_params[k] = proposed_stats_i[k]; + for (size_t k = 0u; k < n_parameters; ++k) + accepted_params[k] = params_init[k]; for (size_t i = 1u; i < n_samples; ++i) { @@ -274,7 +276,10 @@ inline void LFMCMC::run( summary_fun(proposed_stats_i, data_i, this); // Step 4: Compute the hastings ratio using the kernel function - epiworld_double hr = kernel_fun(proposed_stats_i, observed_stats, epsilon, this); + epiworld_double hr = kernel_fun( + proposed_stats_i, observed_stats, epsilon, this + ); + sampled_stats_prob[i] = hr; // Storing data @@ -282,7 +287,7 @@ inline void LFMCMC::run( sampled_stats[i * n_statistics + k] = proposed_stats_i[k]; // Running Hastings ratio - epiworld_double r = runif(); + epiworld_double r = runif(); drawn_prob[i] = r; // Step 5: Update if likely @@ -418,7 +423,7 @@ inline void LFMCMC::get_elapsed( epiworld_double * last_elapsed, std::string * unit_abbr, bool print -) { +) const { // Preparing the result epiworld_double elapsed; diff --git a/inst/include/epiworld/misc.hpp b/inst/include/epiworld/misc.hpp index dbeb73d1..6e0376b1 100644 --- a/inst/include/epiworld/misc.hpp +++ b/inst/include/epiworld/misc.hpp @@ -124,7 +124,7 @@ inline bool IN(const Ta & a, const std::vector< Ta > & b) noexcept * @return int If -1 then it means that none got sampled, otherwise the index * of the entry that got drawn. */ -template +template inline int roulette( const std::vector< TDbl > & probs, Model * m diff --git a/inst/include/epiworld/model-meat.hpp b/inst/include/epiworld/model-meat.hpp index 5df3befb..49160165 100644 --- a/inst/include/epiworld/model-meat.hpp +++ b/inst/include/epiworld/model-meat.hpp @@ -714,7 +714,7 @@ template inline epiworld_double & Model::operator()(std::string pname) { if (parameters.find(pname) == parameters.end()) - throw std::range_error("The parameter "+ pname + "is not in the model."); + throw std::range_error("The parameter '"+ pname + "' is not in the model."); return parameters[pname]; diff --git a/inst/tinytest/test-lfmcmc.R b/inst/tinytest/test-lfmcmc.R index 1425b42d..8eeba0a7 100644 --- a/inst/tinytest/test-lfmcmc.R +++ b/inst/tinytest/test-lfmcmc.R @@ -71,10 +71,27 @@ expect_silent(run_lfmcmc( expect_silent(set_stats_names(lfmcmc_model, get_states(model_sir))) expect_silent(set_par_names(lfmcmc_model, c("Immune recovery", "Infectiousness"))) +# Check printing LFMCMC -------------------------------------------------------- expect_stdout(print(lfmcmc_model)) +expect_stdout(print(lfmcmc_model, burnin = n_samp / 2)) +expect_error(print(lfmcmc_model, burnin = n_samp), "burnin is greater than or equal to the number of samples") +expect_error(print(lfmcmc_model, burnin = n_samp + 50), "burnin is greater than or equal to the number of samples") +expect_error(print(lfmcmc_model, burnin = -n_samp / 2), "argument must be a non-negative integer") +expect_error(print(lfmcmc_model, burnin = "n_samp"), "argument must be an integer") -expect_equal(get_stats_mean(lfmcmc_model), c(284.7140, 0.8485, 713.9375)) -expect_equal(get_params_mean(lfmcmc_model), c(0.3132901, 0.2782186), tolerance = 0.00001) +# Check LFMCMC getters --------------------------------------------------------- +expect_equal(get_n_samples(lfmcmc_model), n_samp) + +expected_stats_mean <- c(284.7140, 0.8485, 713.9375) +expect_equal(get_stats_mean(lfmcmc_model), expected_stats_mean) +expect_equal(get_n_statistics(lfmcmc_model), length(expected_stats_mean)) + +expected_params_mean <- c(0.3133401, 0.2749686) +expect_equal(get_params_mean(lfmcmc_model), expected_params_mean, tolerance = 0.0001) +expect_equal(get_n_parameters(lfmcmc_model), length(expected_params_mean)) + +expect_equal(dim(get_accepted_params(lfmcmc_model)), c(n_samp, length(expected_params_mean))) +expect_equal(dim(get_statistics_hist(lfmcmc_model)), c(n_samp, length(expected_stats_mean))) # Check LFMCMC using factory functions ----------------------------------------- expect_silent(use_proposal_norm_reflective(lfmcmc_model)) @@ -194,3 +211,36 @@ expect_equivalent( c(-5, 2.5), tolerance = 0.1 ) + +# Check functions fail when not passing an LFMCMC object ----------------------- +# Target is 56 tests +expected_error_msg <- "must be an object of class epiworld_lfmcmc" +not_lfmcmc <- c("NOT LFMCMC") + +expect_error(run_lfmcmc(not_lfmcmc, + params_init_ = par0_int, + n_samples_ = n_samp_double, + epsilon_ = epsil_int, + seed = model_seed + ), expected_error_msg) + +expect_error(set_observed_data(not_lfmcmc, obs_data), expected_error_msg) + +expect_error(set_proposal_fun(not_lfmcmc, propfun), expected_error_msg) +expect_error(use_proposal_norm_reflective(not_lfmcmc), expected_error_msg) +expect_error(set_simulation_fun(not_lfmcmc, simfun), expected_error_msg) +expect_error(set_summary_fun(not_lfmcmc, sumfun), expected_error_msg) +expect_error(set_kernel_fun(not_lfmcmc, kernelfun), expected_error_msg) +expect_error(use_kernel_fun_gaussian(not_lfmcmc), expected_error_msg) + +expect_error(set_par_names(not_lfmcmc, c("Par 1", "Par 2")), expected_error_msg) +expect_error(set_stats_names(not_lfmcmc, get_states(model_sir)), expected_error_msg) + +expect_error(get_params_mean(not_lfmcmc), expected_error_msg) +expect_error(get_stats_mean(not_lfmcmc), expected_error_msg) +expect_error(get_accepted_params(not_lfmcmc), expected_error_msg) +expect_error(get_accepted_stats(not_lfmcmc), expected_error_msg) +expect_error(get_statistics_hist(not_lfmcmc), expected_error_msg) +expect_error(get_n_parameters(not_lfmcmc), expected_error_msg) +expect_error(get_n_statistics(not_lfmcmc), expected_error_msg) +expect_error(get_n_samples(not_lfmcmc), expected_error_msg) \ No newline at end of file diff --git a/man/LFMCMC.Rd b/man/LFMCMC.Rd index be54c14a..24623b17 100644 --- a/man/LFMCMC.Rd +++ b/man/LFMCMC.Rd @@ -16,6 +16,12 @@ \alias{get_params_mean} \alias{get_stats_mean} \alias{print.epiworld_lfmcmc} +\alias{get_accepted_params} +\alias{get_accepted_stats} +\alias{get_statistics_hist} +\alias{get_n_parameters} +\alias{get_n_statistics} +\alias{get_n_samples} \title{Likelihood-Free Markhov Chain Monte Carlo (LFMCMC)} \usage{ LFMCMC(model = NULL) @@ -44,7 +50,19 @@ get_params_mean(lfmcmc) get_stats_mean(lfmcmc) -\method{print}{epiworld_lfmcmc}(x, ...) +\method{print}{epiworld_lfmcmc}(x, burnin = 0, ...) + +get_accepted_params(lfmcmc) + +get_accepted_stats(lfmcmc) + +get_statistics_hist(lfmcmc) + +get_n_parameters(lfmcmc) + +get_n_statistics(lfmcmc) + +get_n_samples(lfmcmc) } \arguments{ \item{model}{A model of class \link{epiworld_model} or \code{NULL} (see details).} @@ -67,6 +85,8 @@ get_stats_mean(lfmcmc) \item{x}{LFMCMC model to print} +\item{burnin}{Integer. Number of samples to discard as burnin before computing the summary.} + \item{...}{Ignored} } \value{ @@ -97,6 +117,31 @@ The param means for the given lfmcmc model The stats means for the given lfmcmc model The lfmcmc model + +\itemize{ +\item The function \code{get_accepted_params} returns a matrix of accepted +parameters for the given LFMCMC model. with the number of rows equal to the +number of samples and the number of columns equal to the number of +parameters. +} + +\itemize{ +\item The function \code{get_accepted_stats} returns a matrix of accepted statistics +for the given LFMCMC model. with the number of rows equal to the number of +samples and the number of columns equal to the number of statistics. +} + +\itemize{ +\item The function \code{get_statistics_hist} returns a matrix of statistics +for the given LFMCMC model. with the number of rows equal to the number of +samples and the number of columns equal to the number of statistics. +} + +\itemize{ +\item The functions \code{get_n_parameters}, \code{get_n_statistics}, and \code{get_n_samples} +return the number of parameters, statistics, and samples for the given +LFMCMC model, respectively. +} } \description{ Likelihood-Free Markhov Chain Monte Carlo (LFMCMC) diff --git a/src/cpp11.cpp b/src/cpp11.cpp index 499185ec..24afd8f0 100644 --- a/src/cpp11.cpp +++ b/src/cpp11.cpp @@ -482,10 +482,52 @@ extern "C" SEXP _epiworldR_get_stats_mean_cpp(SEXP lfmcmc) { END_CPP11 } // lfmcmc.cpp -SEXP print_lfmcmc_cpp(SEXP lfmcmc); -extern "C" SEXP _epiworldR_print_lfmcmc_cpp(SEXP lfmcmc) { +SEXP print_lfmcmc_cpp(SEXP lfmcmc, int burnin); +extern "C" SEXP _epiworldR_print_lfmcmc_cpp(SEXP lfmcmc, SEXP burnin) { BEGIN_CPP11 - return cpp11::as_sexp(print_lfmcmc_cpp(cpp11::as_cpp>(lfmcmc))); + return cpp11::as_sexp(print_lfmcmc_cpp(cpp11::as_cpp>(lfmcmc), cpp11::as_cpp>(burnin))); + END_CPP11 +} +// lfmcmc.cpp +SEXP get_statistics_hist_cpp(SEXP lfmcmc); +extern "C" SEXP _epiworldR_get_statistics_hist_cpp(SEXP lfmcmc) { + BEGIN_CPP11 + return cpp11::as_sexp(get_statistics_hist_cpp(cpp11::as_cpp>(lfmcmc))); + END_CPP11 +} +// lfmcmc.cpp +SEXP get_accepted_params_cpp(SEXP lfmcmc); +extern "C" SEXP _epiworldR_get_accepted_params_cpp(SEXP lfmcmc) { + BEGIN_CPP11 + return cpp11::as_sexp(get_accepted_params_cpp(cpp11::as_cpp>(lfmcmc))); + END_CPP11 +} +// lfmcmc.cpp +SEXP get_accepted_stats_cpp(SEXP lfmcmc); +extern "C" SEXP _epiworldR_get_accepted_stats_cpp(SEXP lfmcmc) { + BEGIN_CPP11 + return cpp11::as_sexp(get_accepted_stats_cpp(cpp11::as_cpp>(lfmcmc))); + END_CPP11 +} +// lfmcmc.cpp +int get_n_samples_cpp(SEXP lfmcmc); +extern "C" SEXP _epiworldR_get_n_samples_cpp(SEXP lfmcmc) { + BEGIN_CPP11 + return cpp11::as_sexp(get_n_samples_cpp(cpp11::as_cpp>(lfmcmc))); + END_CPP11 +} +// lfmcmc.cpp +int get_n_statistics_cpp(SEXP lfmcmc); +extern "C" SEXP _epiworldR_get_n_statistics_cpp(SEXP lfmcmc) { + BEGIN_CPP11 + return cpp11::as_sexp(get_n_statistics_cpp(cpp11::as_cpp>(lfmcmc))); + END_CPP11 +} +// lfmcmc.cpp +int get_n_parameters_cpp(SEXP lfmcmc); +extern "C" SEXP _epiworldR_get_n_parameters_cpp(SEXP lfmcmc) { + BEGIN_CPP11 + return cpp11::as_sexp(get_n_parameters_cpp(cpp11::as_cpp>(lfmcmc))); END_CPP11 } // model.cpp @@ -1059,6 +1101,8 @@ static const R_CallMethodDef CallEntries[] = { {"_epiworldR_entity_add_agent_cpp", (DL_FUNC) &_epiworldR_entity_add_agent_cpp, 3}, {"_epiworldR_entity_cpp", (DL_FUNC) &_epiworldR_entity_cpp, 4}, {"_epiworldR_entity_get_agents_cpp", (DL_FUNC) &_epiworldR_entity_get_agents_cpp, 1}, + {"_epiworldR_get_accepted_params_cpp", (DL_FUNC) &_epiworldR_get_accepted_params_cpp, 1}, + {"_epiworldR_get_accepted_stats_cpp", (DL_FUNC) &_epiworldR_get_accepted_stats_cpp, 1}, {"_epiworldR_get_agent_cpp", (DL_FUNC) &_epiworldR_get_agent_cpp, 2}, {"_epiworldR_get_agents_cpp", (DL_FUNC) &_epiworldR_get_agents_cpp, 1}, {"_epiworldR_get_agents_data_ncols_cpp", (DL_FUNC) &_epiworldR_get_agents_data_ncols_cpp, 1}, @@ -1073,7 +1117,10 @@ static const R_CallMethodDef CallEntries[] = { {"_epiworldR_get_hist_total_cpp", (DL_FUNC) &_epiworldR_get_hist_total_cpp, 1}, {"_epiworldR_get_hist_transition_matrix_cpp", (DL_FUNC) &_epiworldR_get_hist_transition_matrix_cpp, 2}, {"_epiworldR_get_hist_virus_cpp", (DL_FUNC) &_epiworldR_get_hist_virus_cpp, 1}, + {"_epiworldR_get_n_parameters_cpp", (DL_FUNC) &_epiworldR_get_n_parameters_cpp, 1}, {"_epiworldR_get_n_replicates_cpp", (DL_FUNC) &_epiworldR_get_n_replicates_cpp, 1}, + {"_epiworldR_get_n_samples_cpp", (DL_FUNC) &_epiworldR_get_n_samples_cpp, 1}, + {"_epiworldR_get_n_statistics_cpp", (DL_FUNC) &_epiworldR_get_n_statistics_cpp, 1}, {"_epiworldR_get_n_tools_cpp", (DL_FUNC) &_epiworldR_get_n_tools_cpp, 1}, {"_epiworldR_get_n_viruses_cpp", (DL_FUNC) &_epiworldR_get_n_viruses_cpp, 1}, {"_epiworldR_get_name_cpp", (DL_FUNC) &_epiworldR_get_name_cpp, 1}, @@ -1086,6 +1133,7 @@ static const R_CallMethodDef CallEntries[] = { {"_epiworldR_get_reproductive_number_cpp", (DL_FUNC) &_epiworldR_get_reproductive_number_cpp, 1}, {"_epiworldR_get_state_agent_cpp", (DL_FUNC) &_epiworldR_get_state_agent_cpp, 1}, {"_epiworldR_get_states_cpp", (DL_FUNC) &_epiworldR_get_states_cpp, 1}, + {"_epiworldR_get_statistics_hist_cpp", (DL_FUNC) &_epiworldR_get_statistics_hist_cpp, 1}, {"_epiworldR_get_stats_mean_cpp", (DL_FUNC) &_epiworldR_get_stats_mean_cpp, 1}, {"_epiworldR_get_today_total_cpp", (DL_FUNC) &_epiworldR_get_today_total_cpp, 1}, {"_epiworldR_get_tool_model_cpp", (DL_FUNC) &_epiworldR_get_tool_model_cpp, 2}, @@ -1106,7 +1154,7 @@ static const R_CallMethodDef CallEntries[] = { {"_epiworldR_print_cpp", (DL_FUNC) &_epiworldR_print_cpp, 2}, {"_epiworldR_print_entity_cpp", (DL_FUNC) &_epiworldR_print_entity_cpp, 1}, {"_epiworldR_print_global_action_cpp", (DL_FUNC) &_epiworldR_print_global_action_cpp, 1}, - {"_epiworldR_print_lfmcmc_cpp", (DL_FUNC) &_epiworldR_print_lfmcmc_cpp, 1}, + {"_epiworldR_print_lfmcmc_cpp", (DL_FUNC) &_epiworldR_print_lfmcmc_cpp, 2}, {"_epiworldR_print_tool_cpp", (DL_FUNC) &_epiworldR_print_tool_cpp, 1}, {"_epiworldR_print_virus_cpp", (DL_FUNC) &_epiworldR_print_virus_cpp, 1}, {"_epiworldR_queuing_off_cpp", (DL_FUNC) &_epiworldR_queuing_off_cpp, 1}, diff --git a/src/lfmcmc.cpp b/src/lfmcmc.cpp index 98453140..eb420e13 100644 --- a/src/lfmcmc.cpp +++ b/src/lfmcmc.cpp @@ -249,11 +249,60 @@ cpp11::writable::doubles get_stats_mean_cpp( [[cpp11::register]] SEXP print_lfmcmc_cpp( - SEXP lfmcmc + SEXP lfmcmc, + int burnin ) { WrapLFMCMC(lfmcmc_ptr)(lfmcmc); - lfmcmc_ptr->print(); + lfmcmc_ptr->print(static_cast(burnin)); return lfmcmc; } +[[cpp11::register]] +SEXP get_statistics_hist_cpp(SEXP lfmcmc) { + + WrapLFMCMC(lfmcmc_ptr)(lfmcmc); + return cpp11::doubles(lfmcmc_ptr->get_statistics_hist()); + +} + +[[cpp11::register]] +SEXP get_accepted_params_cpp(SEXP lfmcmc) { + + WrapLFMCMC(lfmcmc_ptr)(lfmcmc); + return cpp11::doubles(lfmcmc_ptr->get_accepted_params()); + +} + +[[cpp11::register]] +SEXP get_accepted_stats_cpp(SEXP lfmcmc) { + + WrapLFMCMC(lfmcmc_ptr)(lfmcmc); + return cpp11::doubles(lfmcmc_ptr->get_accepted_stats()); + +} + +[[cpp11::register]] +int get_n_samples_cpp(SEXP lfmcmc) { + + WrapLFMCMC(lfmcmc_ptr)(lfmcmc); + return static_cast(lfmcmc_ptr->get_n_samples()); + +} + +[[cpp11::register]] +int get_n_statistics_cpp(SEXP lfmcmc) { + + WrapLFMCMC(lfmcmc_ptr)(lfmcmc); + return static_cast(lfmcmc_ptr->get_n_statistics()); + +} + +[[cpp11::register]] +int get_n_parameters_cpp(SEXP lfmcmc) { + + WrapLFMCMC(lfmcmc_ptr)(lfmcmc); + return static_cast(lfmcmc_ptr->get_n_parameters()); + +} + #undef WrapLFMCMC