From 3e116504fe9ff0833d260857691404b8ed10ad41 Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Thu, 25 Apr 2024 13:47:09 -0600 Subject: [PATCH] Adding the compiled version --- epiworld.hpp | 337 ++++++++++++++++++++++++---------- include/epiworld/epiworld.hpp | 2 +- 2 files changed, 239 insertions(+), 100 deletions(-) diff --git a/epiworld.hpp b/epiworld.hpp index e32dd42c..e18df010 100644 --- a/epiworld.hpp +++ b/epiworld.hpp @@ -16,6 +16,16 @@ #ifndef EPIWORLD_HPP #define EPIWORLD_HPP + +/* Versioning */ +#define EPIWORLD_VERSION_MAJOR 0 +#define EPIWORLD_VERSION_MINOR 1 +#define EPIWORLD_VERSION_PATCH 0 + +static const int epiworld_version_major = EPIWORLD_VERSION_MAJOR; +static const int epiworld_version_minor = EPIWORLD_VERSION_MINOR; +static const int epiworld_version_patch = EPIWORLD_VERSION_PATCH; + namespace epiworld { /*////////////////////////////////////////////////////////////////////////////// @@ -8295,7 +8305,11 @@ inline void Model::run_multiple( std::function*)> fun, bool reset, bool verbose, + #ifdef _OPENMP int nthreads + #else + int + #endif ) { @@ -16444,13 +16458,19 @@ inline ModelSURV::ModelSURV( template class ModelSIRCONN : public epiworld::Model { + private: + + std::vector< epiworld::Agent * > infected; + void update_infected(); + +public: + static const int SUSCEPTIBLE = 0; static const int INFECTED = 1; static const int RECOVERED = 2; -public: - + ModelSIRCONN() {}; ModelSIRCONN( @@ -16491,9 +16511,43 @@ class ModelSIRCONN : public epiworld::Model std::vector< int > queue_ = {} ); + /** + * @brief Get the infected individuals + * @return std::vector< epiworld::Agent * > + */ + size_t get_n_infected() const + { + return infected.size(); + } + }; +template +inline void ModelSIRCONN::update_infected() +{ + + infected.clear(); + infected.reserve(this->size()); + + for (auto & p : this->get_agents()) + { + if (p.get_state() == ModelSIRCONN::INFECTED) + { + infected.push_back(&p); + } + } + + Model::set_rand_binom( + this->get_n_infected(), + static_cast(Model::par("Contact rate"))/ + static_cast(Model::size()) + ); + + return; + +} + template inline ModelSIRCONN & ModelSIRCONN::run( epiworld_fast_uint ndays, @@ -16511,6 +16565,9 @@ inline void ModelSIRCONN::reset() { Model::reset(); + + this->update_infected(); + return; } @@ -16555,26 +16612,21 @@ inline ModelSIRCONN::ModelSIRCONN( ) -> void { - // Sampling how many individuals - m->set_rand_binom( - m->size(), - static_cast( - m->par("Contact rate"))/ - static_cast(m->size()) - ); - int ndraw = m->rbinom(); if (ndraw == 0) return; + ModelSIRCONN * model = dynamic_cast *>(m); + size_t ninfected = model->get_n_infected(); + // Drawing from the set int nviruses_tmp = 0; for (int i = 0; i < ndraw; ++i) { // Now selecting who is transmitting the disease int which = static_cast( - std::floor(m->size() * m->runif()) + std::floor(ninfected * m->runif()) ); /* There is a bug in which runif() returns 1.0. It is rare, but @@ -16584,39 +16636,35 @@ inline ModelSIRCONN::ModelSIRCONN( * https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63176 * */ - if (which == static_cast(m->size())) + if (which == static_cast(ninfected)) --which; + epiworld::Agent & neighbor = *model->infected[which]; + // Can't sample itself - if (which == static_cast(p->get_id())) + if (neighbor.get_id() == p->get_id()) continue; - // If the neighbor is infected, then proceed - auto & neighbor = m->get_agents()[which]; - if (neighbor.get_state() == ModelSIRCONN::INFECTED) - { - - if (neighbor.get_virus() == nullptr) - continue; - - auto & v = neighbor.get_virus(); + // The neighbor is infected because it is on the list! + if (neighbor.get_virus() == nullptr) + continue; - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor.get_transmission_reduction(v, m)) - ; - - m->array_virus_tmp[nviruses_tmp++] = &(*v); - + auto & v = neighbor.get_virus(); - } + #ifdef EPI_DEBUG + if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif + + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor.get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); + } // No virus to compute @@ -16698,6 +16746,17 @@ inline ModelSIRCONN::ModelSIRCONN( model.add_param(transmission_rate, "Transmission rate"); model.add_param(recovery_rate, "Recovery rate"); // model.add_param(prob_reinfection, "Prob. Reinfection"); + + // Adding update function + epiworld::GlobalFun update = [](epiworld::Model * m) -> void + { + ModelSIRCONN * model = dynamic_cast *>(m); + model->update_infected(); + + return; + }; + + model.add_globalevent(update, "Update infected individuals"); // Preparing the virus ------------------------------------------- epiworld::Virus virus(vname); @@ -16782,6 +16841,10 @@ inline ModelSIRCONN & ModelSIRCONN::initial_states( template class ModelSEIRCONN : public epiworld::Model { +private: + std::vector< epiworld::Agent * > infected; + void update_infected(); + public: static const int SUSCEPTIBLE = 0; @@ -16832,8 +16895,35 @@ class ModelSEIRCONN : public epiworld::Model std::vector< int > queue_ = {} ); + size_t get_n_infected() const { return infected.size(); } + }; +template +inline void ModelSEIRCONN::update_infected() +{ + + infected.clear(); + infected.reserve(this->size()); + + for (auto & p : this->get_agents()) + { + if (p.get_state() == ModelSEIRCONN::INFECTED) + { + infected.push_back(&p); + } + } + + Model::set_rand_binom( + this->get_n_infected(), + static_cast(Model::par("Contact rate"))/ + static_cast(Model::size()) + ); + + return; + +} + template inline ModelSEIRCONN & ModelSEIRCONN::run( epiworld_fast_uint ndays, @@ -16852,13 +16942,7 @@ inline void ModelSEIRCONN::reset() { Model::reset(); - - Model::set_rand_binom( - Model::size(), - static_cast( - Model::par("Contact rate"))/ - static_cast(Model::size()) - ); + this->update_infected(); return; @@ -16911,13 +16995,16 @@ inline ModelSEIRCONN::ModelSEIRCONN( if (ndraw == 0) return; + ModelSEIRCONN * model = dynamic_cast *>(m); + size_t ninfected = model->get_n_infected(); + // Drawing from the set int nviruses_tmp = 0; for (int i = 0; i < ndraw; ++i) { // Now selecting who is transmitting the disease int which = static_cast( - std::floor(m->size() * m->runif()) + std::floor(ninfected * m->runif()) ); /* There is a bug in which runif() returns 1.0. It is rare, but @@ -16927,35 +17014,32 @@ inline ModelSEIRCONN::ModelSEIRCONN( * https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63176 * */ - if (which == static_cast(m->size())) + if (which == static_cast(ninfected)) --which; + epiworld::Agent & neighbor = *model->infected[which]; + // Can't sample itself - if (which == static_cast(p->get_id())) + if (neighbor.get_id() == p->get_id()) continue; - // If the neighbor is infected, then proceed - auto & neighbor = m->get_agents()[which]; - if (neighbor.get_state() == ModelSEIRCONN::INFECTED) - { - - auto & v = neighbor.get_virus(); + // The neighbor is infected by construction + auto & v = neighbor.get_virus(); - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor.get_transmission_reduction(v, m)) - ; - - m->array_virus_tmp[nviruses_tmp++] = &(*v); + #ifdef EPI_DEBUG + if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif + + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor.get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); - } } // No virus to compute @@ -17057,6 +17141,22 @@ inline ModelSEIRCONN::ModelSEIRCONN( model.add_state("Infected", update_infected); model.add_state("Recovered"); + // Adding update function + epiworld::GlobalFun update = []( + epiworld::Model * m + ) -> void + { + + ModelSEIRCONN * model = dynamic_cast *>(m); + + model->update_infected(); + + return; + + }; + + model.add_globalevent(update, "Update infected individuals"); + // Preparing the virus ------------------------------------------- epiworld::Virus virus(vname); @@ -18021,6 +18121,10 @@ inline ModelSIRDCONN::ModelSIRDCONN( template class ModelSEIRDCONN : public epiworld::Model { +private: + std::vector< epiworld::Agent * > infected; + void update_infected(); + public: static const int SUSCEPTIBLE = 0; @@ -18029,7 +18133,6 @@ class ModelSEIRDCONN : public epiworld::Model static const int REMOVED = 3; static const int DECEASED = 4; - ModelSEIRDCONN() {}; ModelSEIRDCONN( @@ -18075,8 +18178,36 @@ class ModelSEIRDCONN : public epiworld::Model std::vector< int > queue_ = {} ); + size_t get_n_infected() const + { + return infected.size(); + } + }; +template +inline void ModelSEIRDCONN::update_infected() +{ + infected.clear(); + infected.reserve(this->size()); + + for (auto & p : this->get_agents()) + { + if (p.get_state() == ModelSEIRDCONN::INFECTED) + { + infected.push_back(&p); + } + } + + Model::set_rand_binom( + this->get_n_infected(), + static_cast(Model::par("Contact rate"))/ + static_cast(Model::size()) + ); + + return; +} + template inline ModelSEIRDCONN & ModelSEIRDCONN::run( epiworld_fast_uint ndays, @@ -18096,12 +18227,7 @@ inline void ModelSEIRDCONN::reset() Model::reset(); - Model::set_rand_binom( - Model::size(), - static_cast( - Model::par("Contact rate"))/ - static_cast(Model::size()) - ); + this->update_infected(); return; @@ -18156,13 +18282,19 @@ inline ModelSEIRDCONN::ModelSEIRDCONN( if (ndraw == 0) return; + ModelSEIRDCONN * model = dynamic_cast *>( + m + ); + + size_t ninfected = model->get_n_infected(); + // Drawing from the set int nviruses_tmp = 0; for (int i = 0; i < ndraw; ++i) { // Now selecting who is transmitting the disease int which = static_cast( - std::floor(m->size() * m->runif()) + std::floor(ninfected * m->runif()) ); /* There is a bug in which runif() returns 1.0. It is rare, but @@ -18172,36 +18304,31 @@ inline ModelSEIRDCONN::ModelSEIRDCONN( * https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63176 * */ - if (which == static_cast(m->size())) + if (which == static_cast(ninfected)) --which; + epiworld::Agent & neighbor = *model->infected[which]; + // Can't sample itself - if (which == static_cast(p->get_id())) + if (neighbor.get_id() == p->get_id()) continue; - // If the neighbor is infected, then proceed - auto & neighbor = m->get_agents()[which]; - if (neighbor.get_state() == ModelSEIRDCONN::INFECTED) - { - - const auto & v = neighbor.get_virus(); - - - #ifdef EPI_DEBUG - if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) - throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); - #endif - - /* And it is a function of susceptibility_reduction as well */ - m->array_double_tmp[nviruses_tmp] = - (1.0 - p->get_susceptibility_reduction(v, m)) * - v->get_prob_infecting(m) * - (1.0 - neighbor.get_transmission_reduction(v, m)) - ; - - m->array_virus_tmp[nviruses_tmp++] = &(*v); + // All neighbors in this set are infected by construction + const auto & v = neighbor.get_virus(); + + #ifdef EPI_DEBUG + if (nviruses_tmp >= static_cast(m->array_virus_tmp.size())) + throw std::logic_error("Trying to add an extra element to a temporal array outside of the range."); + #endif - } + /* And it is a function of susceptibility_reduction as well */ + m->array_double_tmp[nviruses_tmp] = + (1.0 - p->get_susceptibility_reduction(v, m)) * + v->get_prob_infecting(m) * + (1.0 - neighbor.get_transmission_reduction(v, m)) + ; + + m->array_virus_tmp[nviruses_tmp++] = &(*v); } // No virus to compute @@ -18318,6 +18445,18 @@ inline ModelSEIRDCONN::ModelSEIRDCONN( model.add_state("Deceased"); + // Adding update function + epiworld::GlobalFun update = [](epiworld::Model * m) -> void + { + ModelSEIRDCONN * model = dynamic_cast *>(m); + model->update_infected(); + + return; + }; + + model.add_globalevent(update, "Update infected individuals"); + + // Preparing the virus ------------------------------------------- epiworld::Virus virus(vname); virus.set_state( diff --git a/include/epiworld/epiworld.hpp b/include/epiworld/epiworld.hpp index 33ecbe54..5d7d7bdc 100644 --- a/include/epiworld/epiworld.hpp +++ b/include/epiworld/epiworld.hpp @@ -20,7 +20,7 @@ /* Versioning */ #define EPIWORLD_VERSION_MAJOR 0 #define EPIWORLD_VERSION_MINOR 1 -#define EPIWORLD_VERSION_PATCH 0 +#define EPIWORLD_VERSION_PATCH 1 static const int epiworld_version_major = EPIWORLD_VERSION_MAJOR; static const int epiworld_version_minor = EPIWORLD_VERSION_MINOR;