From b3cfbbd92b72ea34096e38761a73b718235d1750 Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Mon, 18 Nov 2024 16:20:11 -0700 Subject: [PATCH 01/11] Adding makefile and main --- examples/14-community-hosp/Makefile | 2 + examples/14-community-hosp/main.cpp | 75 +++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+) create mode 100644 examples/14-community-hosp/Makefile create mode 100644 examples/14-community-hosp/main.cpp diff --git a/examples/14-community-hosp/Makefile b/examples/14-community-hosp/Makefile new file mode 100644 index 00000000..0c4d9f21 --- /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/main.cpp b/examples/14-community-hosp/main.cpp new file mode 100644 index 00000000..7a020551 --- /dev/null +++ b/examples/14-community-hosp/main.cpp @@ -0,0 +1,75 @@ +#include "../../epiworld.hpp" + +/*** + * A concrete example would be 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; + +template +inline void my_suscept( + Agent * p, + Model * m + ) +{ + + // Counting how many entities the agent has + auto n_entities = p->get_n_entities(); + + if (n_entities == 0) + { + // You can move agents to an entity + // p->add_entity(m->get_entity("Hospital")); + return; + + } + + Virus * virus = sampler::sample_virus_single(p, m); + + // You can remove them from the entity + // p->rm_entity(m->get_entity("Hospital")); + + if (virus == nullptr) + return; + + p->set_virus(*virus, m); + + 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", my_suscept<>); // State 0 + model.add_state("Infected", default_update_exposed<>); // State 1 + + Entity<> hospital("Hospital"); + hospital.set_distribution(distribute_entity_randomly<>(.5, true, true)); + model.add_entity(hospital); + + // 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); + + // Adding a new population + model.run(100, 1231); + model.print(); + + return 0; + +} \ No newline at end of file From 19980a77cc65572c27c892cd9c2d66e928de4cc8 Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Mon, 18 Nov 2024 17:40:18 -0700 Subject: [PATCH 02/11] Updating example notes --- epiworld.hpp | 10 +- examples/14-community-hosp/README.md | 99 ++++++++++++++++ examples/14-community-hosp/main.cpp | 110 +++++++++++++----- .../epiworld/agent-meat-virus-sampling.hpp | 6 +- include/epiworld/misc.hpp | 2 +- include/epiworld/model-meat.hpp | 2 +- 6 files changed, 193 insertions(+), 36 deletions(-) create mode 100644 examples/14-community-hosp/README.md diff --git a/epiworld.hpp b/epiworld.hpp index f6a3022b..01073f31 100644 --- a/epiworld.hpp +++ b/epiworld.hpp @@ -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 @@ -7661,7 +7661,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 +13078,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 +13237,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 +13405,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/README.md b/examples/14-community-hosp/README.md new file mode 100644 index 00000000..0f4bae93 --- /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 (example-karim) $ make +g++ -std=c++17 -O3 -g main.cpp -o main.o +root ➜ /workspaces/epiworld/examples/14-community-hosp (example-karim) $ ./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 + +An 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, when susceptible (in the community) sample, the sampling excludes individuals who are hospitalized. Likewise, hospitalized + + + ```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 index 7a020551..de3bc136 100644 --- a/examples/14-community-hosp/main.cpp +++ b/examples/14-community-hosp/main.cpp @@ -9,34 +9,88 @@ using namespace epiworld; -template -inline void my_suscept( - Agent * p, - Model * m - ) -{ +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. + */ - // Counting how many entities the agent has - auto n_entities = p->get_n_entities(); +inline void update_susceptible(Agent * p, Model * m) +{ - if (n_entities == 0) + auto virus = sampler_suscept(p, m); + if (virus != nullptr) { - // You can move agents to an entity - // p->add_entity(m->get_entity("Hospital")); - return; - + if (m->par("Prob hospitalization") > m->runif()) + p->set_virus(*virus, m, States::Infected_Hospitalized); + else + p->set_virus(*virus, m, States::Infected); } - Virus * virus = sampler::sample_virus_single(p, m); - // You can remove them from the entity - // p->rm_entity(m->get_entity("Hospital")); - - if (virus == nullptr) - return; + 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") + }; - p->set_virus(*virus, m); + // 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; } @@ -47,12 +101,12 @@ int main() { Model<> model; // model.add_state("Susceptible", default_update_susceptible<>); // State 0 - model.add_state("Susceptible", my_suscept<>); // State 0 - model.add_state("Infected", default_update_exposed<>); // State 1 - - Entity<> hospital("Hospital"); - hospital.set_distribution(distribute_entity_randomly<>(.5, true, true)); - model.add_entity(hospital); + 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"); @@ -66,6 +120,10 @@ int main() { // 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(); diff --git a/include/epiworld/agent-meat-virus-sampling.hpp b/include/epiworld/agent-meat-virus-sampling.hpp index d7289937..bfe668ee 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/misc.hpp b/include/epiworld/misc.hpp index dbeb73d1..6e0376b1 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 5df3befb..49160165 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]; From 82b914bd7c0d0c085a1e98b5121ea880cfce0495 Mon Sep 17 00:00:00 2001 From: George Vega Yon Date: Tue, 19 Nov 2024 03:20:54 +0000 Subject: [PATCH 03/11] Adding examples to makefile --- examples/Makefile | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/Makefile b/examples/Makefile index e9ec42df..8b39894e 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 \ From 31d8724fc88871aa0a6216a50253ea267f867d51 Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Mon, 18 Nov 2024 21:43:15 -0700 Subject: [PATCH 04/11] Bumping version --- epiworld.hpp | 2 +- include/epiworld/epiworld.hpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/epiworld.hpp b/epiworld.hpp index 01073f31..2c1d3620 100644 --- a/epiworld.hpp +++ b/epiworld.hpp @@ -19,7 +19,7 @@ /* Versioning */ #define EPIWORLD_VERSION_MAJOR 0 #define EPIWORLD_VERSION_MINOR 4 -#define EPIWORLD_VERSION_PATCH 3 +#define EPIWORLD_VERSION_PATCH 4 static const int epiworld_version_major = EPIWORLD_VERSION_MAJOR; static const int epiworld_version_minor = EPIWORLD_VERSION_MINOR; diff --git a/include/epiworld/epiworld.hpp b/include/epiworld/epiworld.hpp index 6a060f7e..b4c0782d 100644 --- a/include/epiworld/epiworld.hpp +++ b/include/epiworld/epiworld.hpp @@ -19,7 +19,7 @@ /* Versioning */ #define EPIWORLD_VERSION_MAJOR 0 #define EPIWORLD_VERSION_MINOR 4 -#define EPIWORLD_VERSION_PATCH 3 +#define EPIWORLD_VERSION_PATCH 4 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 From d7e04a4a5d232d254b21511690bb0caabc1e1d84 Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Thu, 21 Nov 2024 17:13:21 -0700 Subject: [PATCH 05/11] Multiple patches --- epiworld.hpp | 76 +++++++++++++------ include/epiworld/epiworld.hpp | 4 +- include/epiworld/math/lfmcmc/lfmcmc-bones.hpp | 8 +- .../math/lfmcmc/lfmcmc-meat-print.hpp | 39 +++++++--- include/epiworld/math/lfmcmc/lfmcmc-meat.hpp | 17 +++-- 5 files changed, 100 insertions(+), 44 deletions(-) diff --git a/epiworld.hpp b/epiworld.hpp index 2c1d3620..ee13cbb6 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 4 +#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; @@ -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,29 +2089,46 @@ inline void LFMCMC::get_elapsed( #define LFMCMC_MEAT_PRINT_HPP template -inline void LFMCMC::print() +inline void LFMCMC::print(size_t burnin) const { + std::vector< epiworld_double > summ_params(n_parameters * 3, 0.0); std::vector< epiworld_double > summ_stats(n_statistics * 3, 0.0); + size_t n_samples_print = n_samples; + if (burnin > 0) + { + if (burnin >= n_samples) + throw std::length_error( + "The burnin is greater than the number of samples." + ); + + n_samples_print = n_samples - burnin; + + } + + epiworld_double n_samples_dbl = static_cast< epiworld_double >( + n_samples_print + ); + 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; + summ_params[k * 3] += par_i[i]/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)]; } @@ -2106,20 +2136,20 @@ inline void LFMCMC::print() { // 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; + summ_stats[k * 3] += stat_k[i]/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/epiworld.hpp b/include/epiworld/epiworld.hpp index b4c0782d..12801802 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 4 +#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; diff --git a/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp b/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp index 7002aa24..ab86172b 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 363fb766..1bd5879e 100755 --- a/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp +++ b/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp @@ -2,29 +2,46 @@ #define LFMCMC_MEAT_PRINT_HPP template -inline void LFMCMC::print() +inline void LFMCMC::print(size_t burnin) const { + std::vector< epiworld_double > summ_params(n_parameters * 3, 0.0); std::vector< epiworld_double > summ_stats(n_statistics * 3, 0.0); + size_t n_samples_print = n_samples; + if (burnin > 0) + { + if (burnin >= n_samples) + throw std::length_error( + "The burnin is greater than the number of samples." + ); + + n_samples_print = n_samples - burnin; + + } + + epiworld_double n_samples_dbl = static_cast< epiworld_double >( + n_samples_print + ); + 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; + summ_params[k * 3] += par_i[i]/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)]; } @@ -32,20 +49,20 @@ inline void LFMCMC::print() { // 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; + summ_stats[k * 3] += stat_k[i]/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 c7565121..c515d50c 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; From 94887e0ab23961394c84b6aace2fd72b220ba726 Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Fri, 22 Nov 2024 16:07:53 -0700 Subject: [PATCH 06/11] Update README.md --- examples/14-community-hosp/README.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/14-community-hosp/README.md b/examples/14-community-hosp/README.md index 0f4bae93..088aa54c 100644 --- a/examples/14-community-hosp/README.md +++ b/examples/14-community-hosp/README.md @@ -17,9 +17,9 @@ The model has the following parameters: Here is an example of the run ```bash -root ➜ /workspaces/epiworld/examples/14-community-hosp (example-karim) $ make +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 (example-karim) $ ./main.o +root ➜ /workspaces/epiworld/examples/14-community-hosp $ ./main.o _________________________________________________________________________ Running the model... ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done. @@ -65,9 +65,9 @@ Transition Probabilities: ## Notes -An few key observations from this example. +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, when susceptible (in the community) sample, the sampling excludes individuals who are hospitalized. Likewise, hospitalized +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 From 8fe73e953d18710404bdade069cc6b53b7df771a Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Fri, 22 Nov 2024 16:08:22 -0700 Subject: [PATCH 07/11] Update main.cpp --- examples/14-community-hosp/main.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/14-community-hosp/main.cpp b/examples/14-community-hosp/main.cpp index de3bc136..bc9d6b3e 100644 --- a/examples/14-community-hosp/main.cpp +++ b/examples/14-community-hosp/main.cpp @@ -1,7 +1,7 @@ #include "../../epiworld.hpp" /*** - * A concrete example would be a model that includes two populations. + * 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. @@ -130,4 +130,4 @@ int main() { return 0; -} \ No newline at end of file +} From 561912751ebd3ea7bc02045853ff821ba8c793ac Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Mon, 25 Nov 2024 09:11:01 -0700 Subject: [PATCH 08/11] Adding a print to the lfmcmc test --- tests/00-lfmcmc.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/00-lfmcmc.cpp b/tests/00-lfmcmc.cpp index 4f1dfcea..b793f02e 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(); From af9056e3318ed9b1933f7f816e6c60ad59d75184 Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Mon, 25 Nov 2024 09:17:49 -0700 Subject: [PATCH 09/11] Adding a note about the stored data in print --- include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp b/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp index 1bd5879e..6b7c2fe5 100755 --- a/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp +++ b/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp @@ -5,6 +5,8 @@ template inline void LFMCMC::print(size_t burnin) const { + // The summary statistics will have 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); @@ -58,7 +60,6 @@ inline void LFMCMC::print(size_t burnin) const // Computing the 95% Credible interval std::sort(stat_k.begin(), stat_k.end()); - summ_stats[k * 3 + 1u] = stat_k[std::floor(.025 * n_samples_dbl)]; summ_stats[k * 3 + 2u] = From 0f853d9611aed7122c5871995ef8640437accaa1 Mon Sep 17 00:00:00 2001 From: Andrew Pulsipher Date: Mon, 25 Nov 2024 11:10:42 -0700 Subject: [PATCH 10/11] Add test of correct and error print burning in tests/00-lfmcmc.cpp --- tests/00-lfmcmc.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/00-lfmcmc.cpp b/tests/00-lfmcmc.cpp index b793f02e..9a9dac24 100644 --- a/tests/00-lfmcmc.cpp +++ b/tests/00-lfmcmc.cpp @@ -73,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 From 310652dc94d67005a1fb33874f27dc4d56fe7630 Mon Sep 17 00:00:00 2001 From: Andrew Pulsipher Date: Mon, 25 Nov 2024 12:13:46 -0700 Subject: [PATCH 11/11] Fix lfmcmc printing error and add descriptive comments --- epiworld.hpp | 16 ++++++++++------ .../epiworld/math/lfmcmc/lfmcmc-meat-print.hpp | 17 ++++++++++------- 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/epiworld.hpp b/epiworld.hpp index ee13cbb6..0c71e9cd 100644 --- a/epiworld.hpp +++ b/epiworld.hpp @@ -2092,15 +2092,18 @@ template 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 the number of samples." + "The burnin is greater than or equal to the number of samples." ); n_samples_print = n_samples - burnin; @@ -2111,6 +2114,7 @@ inline void LFMCMC::print(size_t burnin) const n_samples_print ); + // Compute parameter summary values for (size_t k = 0u; k < n_parameters; ++k) { @@ -2118,8 +2122,8 @@ inline void LFMCMC::print(size_t burnin) const 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_dbl; + 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 @@ -2132,6 +2136,7 @@ inline void LFMCMC::print(size_t burnin) const } + // Compute statistics summary values for (size_t k = 0u; k < n_statistics; ++k) { @@ -2139,13 +2144,12 @@ inline void LFMCMC::print(size_t burnin) const 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_dbl; + 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 * n_samples_dbl)]; summ_stats[k * 3 + 2u] = diff --git a/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp b/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp index 6b7c2fe5..b2860a08 100755 --- a/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp +++ b/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp @@ -5,17 +5,18 @@ template inline void LFMCMC::print(size_t burnin) const { - // The summary statistics will have three values: mean, the 2.5% - // quantile, and the 97.5% quantile + // 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 the number of samples." + "The burnin is greater than or equal to the number of samples." ); n_samples_print = n_samples - burnin; @@ -26,6 +27,7 @@ inline void LFMCMC::print(size_t burnin) const n_samples_print ); + // Compute parameter summary values for (size_t k = 0u; k < n_parameters; ++k) { @@ -33,8 +35,8 @@ inline void LFMCMC::print(size_t burnin) const 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_dbl; + 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 @@ -47,6 +49,7 @@ inline void LFMCMC::print(size_t burnin) const } + // Compute statistics summary values for (size_t k = 0u; k < n_statistics; ++k) { @@ -54,8 +57,8 @@ inline void LFMCMC::print(size_t burnin) const 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_dbl; + 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