Skip to content

Commit

Permalink
Multiple patches to lfmcmc (#34)
Browse files Browse the repository at this point in the history
* 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 <[email protected]>
  • Loading branch information
gvegayon and apulsipher authored Nov 25, 2024
1 parent 27c32a5 commit bb8d6e8
Show file tree
Hide file tree
Showing 13 changed files with 364 additions and 62 deletions.
96 changes: 65 additions & 31 deletions epiworld.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -1239,7 +1239,7 @@ class LFMCMC {
epiworld_double * last_elapsed,
std::string * unit_abbr,
bool print
);
) const;

void chrono_start();
void chrono_end();
Expand Down Expand Up @@ -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;

};

Expand Down Expand Up @@ -1517,7 +1521,7 @@ class LFMCMC {
epiworld_double * last_elapsed,
std::string * unit_abbr,
bool print
);
) const;

void chrono_start();
void chrono_end();
Expand Down Expand Up @@ -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;

};

Expand Down Expand Up @@ -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)
{
Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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)];

}

Expand Down Expand Up @@ -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];

Expand Down Expand Up @@ -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 = {}
)
Expand Down Expand Up @@ -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 = {}
)
Expand Down Expand Up @@ -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)
{

Expand Down
2 changes: 2 additions & 0 deletions examples/14-community-hosp/Makefile
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
main.o: main.cpp
g++ -std=c++17 -O3 -g main.cpp -o main.o
99 changes: 99 additions & 0 deletions examples/14-community-hosp/README.md
Original file line number Diff line number Diff line change
@@ -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;
}
```
Loading

0 comments on commit bb8d6e8

Please sign in to comment.