Multiple patches to lfmcmc

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<typename TSeq, typename TDbl>
+template<typename TSeq = EPI_DEFAULT_TSEQ, typename TDbl = epiworld_double >
 inline int roulette(
     const std::vector< TDbl > & probs,
     Model<TSeq> * 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<TData>::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,15 +1877,18 @@ inline void LFMCMC<TData>::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
         for (size_t k = 0u; k < n_statistics; ++k)
             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<TData>::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<TData>::get_elapsed(
 #define LFMCMC_MEAT_PRINT_HPP
 
 template<typename TData>
-inline void LFMCMC<TData>::print()
+inline void LFMCMC<TData>::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<epiworld_double>(n_samples))];
+            par_i[std::floor(.025 * n_samples_dbl)];
         summ_params[k * 3 + 2u] = 
-            par_i[std::floor(.975 * static_cast<epiworld_double>(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<epiworld_double>(n_samples))];
+            stat_k[std::floor(.025 * n_samples_dbl)];
         summ_stats[k * 3 + 2u] = 
-            stat_k[std::floor(.975 * static_cast<epiworld_double>(n_samples))];
+            stat_k[std::floor(.975 * n_samples_dbl)];
 
     }
 
@@ -7661,7 +7695,7 @@ template<typename TSeq>
 inline epiworld_double & Model<TSeq>::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<TSeq>* of the selected virus. If none selected (or none
  * available,) returns a nullptr;
  */
-template<typename TSeq>
+template<typename TSeq = EPI_DEFAULT_TSEQ>
 inline std::function<void(Agent<TSeq>*,Model<TSeq>*)> make_update_susceptible(
     std::vector< epiworld_fast_uint > exclude = {}
     )
@@ -13237,7 +13271,7 @@ inline std::function<void(Agent<TSeq>*,Model<TSeq>*)> make_update_susceptible(
  * @return Virus<TSeq>* of the selected virus. If none selected (or none
  * available,) returns a nullptr;
  */
-template<typename TSeq = int>
+template<typename TSeq = EPI_DEFAULT_TSEQ>
 inline std::function<Virus<TSeq>*(Agent<TSeq>*,Model<TSeq>*)> make_sample_virus_neighbors(
     std::vector< epiworld_fast_uint > exclude = {}
 )
@@ -13405,7 +13439,7 @@ inline std::function<Virus<TSeq>*(Agent<TSeq>*,Model<TSeq>*)> make_sample_virus_
  * @return Virus<TSeq>* of the selected virus. If none selected (or none
  * available,) returns a nullptr;
  */
-template<typename TSeq = int>
+template<typename TSeq = EPI_DEFAULT_TSEQ>
 inline Virus<TSeq> * sample_virus_single(Agent<TSeq> * p, Model<TSeq> * m)
 {
 

examples/14-community-hosp/Makefile

@@ -0,0 +1,2 @@
+main.o: main.cpp
+	g++ -std=c++17 -O3 -g main.cpp -o main.o

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<int> * p, Model<int> * 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;
+
+    }
+    ```