From bb8d6e8bc3436ad28839da5affabfb6975365c2c Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Mon, 25 Nov 2024 13:02:41 -0700 Subject: [PATCH] Multiple patches to lfmcmc (#34) * Adding makefile and main * Updating example notes * Adding examples to makefile * Bumping version * Multiple patches * Update README.md * Update main.cpp * Adding a print to the lfmcmc test * Adding a note about the stored data in print * Add test of correct and error print burnin in tests/00-lfmcmc.cpp * Fix lfmcmc printing error and add descriptive comments --------- Co-authored-by: Andrew Pulsipher --- epiworld.hpp | 96 +++++++++---- examples/14-community-hosp/Makefile | 2 + examples/14-community-hosp/README.md | 99 +++++++++++++ examples/14-community-hosp/main.cpp | 133 ++++++++++++++++++ examples/Makefile | 3 +- .../epiworld/agent-meat-virus-sampling.hpp | 6 +- include/epiworld/epiworld.hpp | 6 +- include/epiworld/math/lfmcmc/lfmcmc-bones.hpp | 8 +- .../math/lfmcmc/lfmcmc-meat-print.hpp | 49 +++++-- include/epiworld/math/lfmcmc/lfmcmc-meat.hpp | 17 ++- include/epiworld/misc.hpp | 2 +- include/epiworld/model-meat.hpp | 2 +- tests/00-lfmcmc.cpp | 3 + 13 files changed, 364 insertions(+), 62 deletions(-) create mode 100644 examples/14-community-hosp/Makefile create mode 100644 examples/14-community-hosp/README.md create mode 100644 examples/14-community-hosp/main.cpp diff --git a/epiworld.hpp b/epiworld.hpp index f6a3022b5..0c71e9cd2 100644 --- a/epiworld.hpp +++ b/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; @@ -595,7 +595,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 @@ -1239,7 +1239,7 @@ class LFMCMC { epiworld_double * last_elapsed, std::string * unit_abbr, bool print - ); + ) const; void chrono_start(); void chrono_end(); @@ -1297,13 +1297,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; }; @@ -1517,7 +1521,7 @@ class LFMCMC { epiworld_double * last_elapsed, std::string * unit_abbr, bool print - ); + ) const; void chrono_start(); void chrono_end(); @@ -1575,13 +1579,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; }; @@ -1842,14 +1850,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) { @@ -1867,7 +1877,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 @@ -1875,7 +1888,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 @@ -2011,7 +2024,7 @@ inline void LFMCMC::get_elapsed( epiworld_double * last_elapsed, std::string * unit_abbr, bool print -) { +) const { // Preparing the result epiworld_double elapsed; @@ -2076,50 +2089,71 @@ inline void LFMCMC::get_elapsed( #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)]; } @@ -7661,7 +7695,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]; @@ -13078,7 +13112,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 = {} ) @@ -13237,7 +13271,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 = {} ) @@ -13405,7 +13439,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/examples/14-community-hosp/Makefile b/examples/14-community-hosp/Makefile new file mode 100644 index 000000000..0c4d9f21d --- /dev/null +++ b/examples/14-community-hosp/Makefile @@ -0,0 +1,2 @@ +main.o: main.cpp + g++ -std=c++17 -O3 -g main.cpp -o main.o \ No newline at end of file diff --git a/examples/14-community-hosp/README.md b/examples/14-community-hosp/README.md new file mode 100644 index 000000000..088aa54c3 --- /dev/null +++ b/examples/14-community-hosp/README.md @@ -0,0 +1,99 @@ +# Community-hospital model + +In this model, we have three states: + +- Susceptible individuals (in the community), +- Infected individuals (in the community), and +- Infected hospitalized individuals. + +Susceptible individuals may become hospitalized or not (so they have two possible transitions), and infected individuals may recover or (if hospitalized) stay infected but be discharged (so they go back as infected but in the community). + +The model has the following parameters: + +- Prob hospitalization: 0.1 +- Prob recovery: 0.33 +- Discharge infected: 0.1 + +Here is an example of the run + +```bash +root ➜ /workspaces/epiworld/examples/14-community-hosp $ make +g++ -std=c++17 -O3 -g main.cpp -o main.o +root ➜ /workspaces/epiworld/examples/14-community-hosp $ ./main.o +_________________________________________________________________________ +Running the model... +||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done. + done. +________________________________________________________________________________ +________________________________________________________________________________ +SIMULATION STUDY + +Name of the model : (none) +Population size : 1000 +Agents' data : (none) +Number of entities : 0 +Days (duration) : 100 (of 100) +Number of viruses : 1 +Last run elapsed t : 4.00ms +Last run speed : 20.81 million agents x day / second +Rewiring : off + +Global events: + (none) + +Virus(es): + - MRSA + +Tool(s): + (none) + +Model parameters: + - Discharge infected : 0.1000 + - Prob hospitalization : 0.1000 + - Prob recovery : 0.3300 + +Distribution of the population at time 100: + - (0) Susceptible : 990 -> 937 + - (1) Infected : 10 -> 49 + - (2) Infected (hospitalized) : 0 -> 14 + +Transition Probabilities: + - Susceptible 0.98 0.02 0.00 + - Infected 0.32 0.61 0.07 + - Infected (hospitalized) 0.35 0.07 0.58 +``` + +## Notes + +A few key observations from this example. + +1. **We have a sampling function that exclude population**. These two functions are used in the update functions so the sampling excludes hospitalized individuals: + + + ```cpp + // A sampler that excludes infected from the hospital + auto sampler_suscept = sampler::make_sample_virus_neighbors<>( + {States::Infected_Hospitalized} + ); + ``` + +2. **All update functions are built from scratch**. For instance, susceptible individuals are updated according to the following function: + + ```cpp + inline void update_susceptible(Agent * p, Model * m) + { + + auto virus = sampler_suscept(p, m); + if (virus != nullptr) + { + if (m->par("Prob hospitalization") > m->runif()) + p->set_virus(*virus, m, States::Infected_Hospitalized); + else + p->set_virus(*virus, m, States::Infected); + } + + + return; + + } + ``` diff --git a/examples/14-community-hosp/main.cpp b/examples/14-community-hosp/main.cpp new file mode 100644 index 000000000..bc9d6b3e7 --- /dev/null +++ b/examples/14-community-hosp/main.cpp @@ -0,0 +1,133 @@ +#include "../../epiworld.hpp" + +/*** + * A concrete example with a model that includes two populations. + * - One, a ‘community’ population and the other, + * - a ‘hospital’ population. + * In this model, individuals can move from the community to the hospital (admission) and move from the hospital back into the community (discharge). In both settings, there can be an infectious disease process (e.g. SIS), but we would assume that transmission does not occur between the community and the hospital (of course, this could be relaxed in the future). But through admission and discharge, infections in the community impact dynamics in the hospital and the reverse is true as well. + */ + +using namespace epiworld; + +enum States : size_t { + Susceptible = 0u, + Infected, + Infected_Hospitalized +}; + +// A sampler that excludes infected from the hospital +auto sampler_suscept = sampler::make_sample_virus_neighbors<>( + {States::Infected_Hospitalized} + ); + +/** + * - Susceptibles only live in the community. + * - Infection from individuals in the community only. + * - Once they become infected, they may be hospitalized or not. + */ + +inline void update_susceptible(Agent * p, Model * m) +{ + + auto virus = sampler_suscept(p, m); + if (virus != nullptr) + { + if (m->par("Prob hospitalization") > m->runif()) + p->set_virus(*virus, m, States::Infected_Hospitalized); + else + p->set_virus(*virus, m, States::Infected); + } + + + return; + +} + +/** + * Infected individuals may: + * + * - Stay the same + * - Recover + * - Be hospitalized + * + * Notice that the roulette makes the probabilities to sum to 1. + */ +inline void update_infected(Agent * p, Model * m) +{ + + // Vector of probabilities + std::vector< epiworld_double > probs = { + m->par("Prob hospitalization"), + m->par("Prob recovery") + }; + + // Sampling: + // - (-1) Nothing happens + // - (0) Hospitalization + // - (1) Recovery + int res = roulette<>(probs, m); + + if (res == 0) + p->change_state(m, States::Infected_Hospitalized); + else if (res == 1) + p->rm_virus(m, States::Susceptible); + + return; + +} + +/** + * Infected individuals who are hospitalized may: + * - Stay infected. + * - Recover (and then be discharged) + * - Stay the same and be discharged. + */ +inline void update_infected_hospitalized(Agent * p, Model * m) +{ + + if (m->par("Prob recovery") > m->runif()) { + p->rm_virus(m, States::Susceptible); + } else if (m->par("Discharge infected") > m->runif()) { + p->change_state(m, States::Infected); + } + + return; + +} + +int main() { + + // Using the mixing model as a baseline + Model<> model; + + // model.add_state("Susceptible", default_update_susceptible<>); // State 0 + model.add_state("Susceptible", update_susceptible); // State 0 + model.add_state("Infected", update_infected); // State 1 + model.add_state( + "Infected (hospitalized)", + update_infected_hospitalized + ); // State 2 + + // Adding a new virus + Virus<> mrsa("MRSA"); + mrsa.set_state(1, 0, 0); + mrsa.set_prob_infecting(.1); + mrsa.set_prob_recovery(.33); + mrsa.set_distribution(distribute_virus_randomly<>(0.01)); + + model.add_virus(mrsa); + + // Add a population + model.agents_smallworld(1000, 4, 0.1, false); + + model.add_param(0.1, "Prob hospitalization"); + model.add_param(0.33, "Prob recovery"); + model.add_param(0.1, "Discharge infected"); + + // Adding a new population + model.run(100, 1231); + model.print(); + + return 0; + +} diff --git a/examples/Makefile b/examples/Makefile index e9ec42df2..8b39894ee 100755 --- a/examples/Makefile +++ b/examples/Makefile @@ -8,7 +8,8 @@ INTELFLAGS=-g -O2 ALL_EXAMPLES=00-hello-world 01-sir 02-sir_multiple_runs 02b-sir_multiple_runs \ 03-simple-sir 04-advanced-usage 05-user-data 06-sir-omp 06b-sir-omp \ - 07-surveillance 10-likelihood-free-mcmc + 07-surveillance 10-likelihood-free-mcmc 11-entities 13-genint \ + 14-community-hosp all-examples: for f in $(ALL_EXAMPLES); do \ diff --git a/include/epiworld/agent-meat-virus-sampling.hpp b/include/epiworld/agent-meat-virus-sampling.hpp index d72899377..bfe668eeb 100644 --- a/include/epiworld/agent-meat-virus-sampling.hpp +++ b/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/include/epiworld/epiworld.hpp b/include/epiworld/epiworld.hpp index 6a060f7eb..12801802b 100644 --- a/include/epiworld/epiworld.hpp +++ b/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/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp b/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp index 7002aa247..ab86172b8 100755 --- a/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp +++ b/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/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp b/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp index 363fb766f..b2860a08c 100755 --- a/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp +++ b/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/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp b/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp index c7565121e..c515d50ca 100755 --- a/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp +++ b/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/include/epiworld/misc.hpp b/include/epiworld/misc.hpp index dbeb73d14..6e0376b1c 100644 --- a/include/epiworld/misc.hpp +++ b/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/include/epiworld/model-meat.hpp b/include/epiworld/model-meat.hpp index 5df3befb5..49160165d 100644 --- a/include/epiworld/model-meat.hpp +++ b/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/tests/00-lfmcmc.cpp b/tests/00-lfmcmc.cpp index 4f1dfcea6..9a9dac24d 100644 --- a/tests/00-lfmcmc.cpp +++ b/tests/00-lfmcmc.cpp @@ -64,6 +64,7 @@ EPIWORLD_TEST_CASE("LFMCMC", "[Basic example]") { model.print(); + model.print(50000); auto params_means = model.get_params_mean(); auto stats_means = model.get_stats_mean(); @@ -72,6 +73,8 @@ EPIWORLD_TEST_CASE("LFMCMC", "[Basic example]") { std::vector expected = {5.0, 1.5}; REQUIRE_THAT(params_means, Catch::Approx(expected).margin(0.5)); REQUIRE_THAT(stats_means, Catch::Approx(expected).margin(0.5)); + REQUIRE_THROWS(model.print(200000)); + REQUIRE_NOTHROW(model.print(50000)); #endif #ifndef CATCH_CONFIG_MAIN