diff --git a/epiworld.hpp b/epiworld.hpp index e32dd42c..9f6eda64 100644 --- a/epiworld.hpp +++ b/epiworld.hpp @@ -6082,9 +6082,9 @@ class Model { * AgentsSample::AgentsSample(Model) these vectors are allocated. */ ///@{ - std::vector< Agent * > sampled_population; + std::shared_ptr< std::vector< Agent * > > sampled_population; size_t sampled_population_n = 0u; - std::vector< size_t > population_left; + std::shared_ptr< std::vector< size_t > > population_left; size_t population_left_n = 0u; ///@} @@ -14621,10 +14621,10 @@ class AgentsSample { size_t sample_size = 0u; - std::vector< Agent* > * agents = nullptr; ///< Pointer to sample of agents + std::shared_ptr* > > agents = nullptr; ///< Pointer to sample of agents size_t * agents_n = nullptr; ///< Size of sample of agents - std::vector< size_t > * agents_left = nullptr; ///< Pointer to agents left (iota) + std::shared_ptr > agents_left = nullptr; ///< Pointer to agents left (iota) size_t * agents_left_n = nullptr; ///< Size of agents left Model * model = nullptr; ///< Extracts runif() and (if the case) population. @@ -14651,13 +14651,17 @@ class AgentsSample { ); AgentsSample( - Model * model, Entity & entity_, size_t n, + Model * model, + Entity & entity_, + size_t n, std::vector< size_t > states_ = {}, bool truncate = false ); AgentsSample( - Model * model, Agent & agent_, size_t n, + Model * model, + Agent & agent_, + size_t n, std::vector< size_t > states_ = {}, bool truncate = false ); diff --git a/examples/11-entities/Makefile b/examples/11-entities/Makefile new file mode 100644 index 00000000..86dbde3c --- /dev/null +++ b/examples/11-entities/Makefile @@ -0,0 +1,2 @@ +main.o: main.cpp + g++ -std=c++17 -g main.cpp -o main.o \ No newline at end of file diff --git a/examples/11-entities/main.cpp b/examples/11-entities/main.cpp new file mode 100644 index 00000000..e6ab094d --- /dev/null +++ b/examples/11-entities/main.cpp @@ -0,0 +1,26 @@ +#define EPI_DEBUG +#include "../../include/epiworld/epiworld.hpp" + +using namespace epiworld; + +int main() { + + epimodels::ModelSEIREntitiesConn model( + "Flu", // std::string vname, + 10000, // epiworld_fast_uint n, + 0.01,// epiworld_double prevalence, + 4.0,// epiworld_double contact_rate, + 0.1,// epiworld_double transmission_rate, + 4.0,// epiworld_double avg_incubation_days, + 1.0/7.0,// epiworld_double recovery_rate, + {.1, .1, .8},// std::vector< epiworld_double > entities, + {"A", "B", "C"}// std::vector< std::string > entities_names + ); + + // Running and checking the results + model.run(50, 123); + model.print(); + + return 0; + +} diff --git a/include/epiworld/agent-bones.hpp b/include/epiworld/agent-bones.hpp index 1e83c069..cd703a32 100644 --- a/include/epiworld/agent-bones.hpp +++ b/include/epiworld/agent-bones.hpp @@ -102,9 +102,9 @@ class Agent { std::vector< ToolPtr > tools; epiworld_fast_uint n_tools = 0u; - std::vector< Agent * > sampled_agents; + std::vector< Agent * > sampled_agents = {}; size_t sampled_agents_n = 0u; - std::vector< size_t > sampled_agents_left; + std::vector< size_t > sampled_agents_left = {}; size_t sampled_agents_left_n = 0u; int date_last_build_sample = -99; diff --git a/include/epiworld/agentssample-bones.hpp b/include/epiworld/agentssample-bones.hpp index c8462c58..3026e3aa 100644 --- a/include/epiworld/agentssample-bones.hpp +++ b/include/epiworld/agentssample-bones.hpp @@ -30,10 +30,10 @@ class AgentsSample { size_t sample_size = 0u; - std::vector< Agent* > * agents = nullptr; ///< Pointer to sample of agents + std::vector< Agent* >* agents = nullptr; ///< Pointer to sample of agents size_t * agents_n = nullptr; ///< Size of sample of agents - std::vector< size_t > * agents_left = nullptr; ///< Pointer to agents left (iota) + std::vector< size_t >* agents_left = nullptr; ///< Pointer to agents left (iota) size_t * agents_left_n = nullptr; ///< Size of agents left Model * model = nullptr; ///< Extracts runif() and (if the case) population. @@ -49,7 +49,7 @@ class AgentsSample { public: // Not available (for now) - AgentsSample() = delete; ///< Default constructor + AgentsSample() = delete; ///< Default constructor AgentsSample(const AgentsSample & a) = delete; ///< Copy constructor AgentsSample(AgentsSample && a) = delete; ///< Move constructor @@ -60,13 +60,17 @@ class AgentsSample { ); AgentsSample( - Model * model, Entity & entity_, size_t n, + Model * model, + Entity & entity_, + size_t n, std::vector< size_t > states_ = {}, bool truncate = false ); AgentsSample( - Model * model, Agent & agent_, size_t n, + Model * model, + Agent & agent_, + size_t n, std::vector< size_t > states_ = {}, bool truncate = false ); @@ -228,19 +232,22 @@ inline AgentsSample::AgentsSample( agents->resize(n); size_t i_obs = 0u; - for (size_t i = 0u; i < agents_in_entities; ++i) + for (size_t i = 0u; i < sample_size; ++i) { + + // Sampling a single agent from the set of entities int jth = std::floor(model->runif() * agents_in_entities); for (size_t e = 0u; e < cum_agents_count.size(); ++e) { + // Are we in the limit? if (jth <= cum_agents_count[e]) { size_t agent_idx = 0u; if (e == 0) // From the first group - agent_idx = entities_a[e][jth]; + agent_idx = entities_a[e][jth]->get_id(); else - agent_idx = entities_a[e][jth - cum_agents_count[e - 1]]; + agent_idx = entities_a[e][jth - cum_agents_count[e - 1]]->get_id(); // Getting the state size_t state = model->population[agent_idx].get_state(); @@ -252,7 +259,7 @@ inline AgentsSample::AgentsSample( continue; } - agents->operator[](i_obs++) = agent_idx; + agents->operator[](i_obs++) = &(model->population[agent_idx]); break; } diff --git a/include/epiworld/model-meat.hpp b/include/epiworld/model-meat.hpp index 030e2638..dbe397b2 100644 --- a/include/epiworld/model-meat.hpp +++ b/include/epiworld/model-meat.hpp @@ -2115,6 +2115,7 @@ inline void Model::reset() { // Re distributing tools and virus dist_virus(); dist_tools(); + dist_entities(); // Distributing initial state, if specified initial_states_fun(this); diff --git a/include/epiworld/models/models.hpp b/include/epiworld/models/models.hpp index 72cfc860..e3674f39 100644 --- a/include/epiworld/models/models.hpp +++ b/include/epiworld/models/models.hpp @@ -19,6 +19,7 @@ namespace epimodels { #include "seirdconnected.hpp" #include "sirlogit.hpp" #include "diffnet.hpp" + #include "seirentitiesconnected.hpp" } diff --git a/include/epiworld/models/seirentitiesconnected.hpp b/include/epiworld/models/seirentitiesconnected.hpp new file mode 100644 index 00000000..c3e77e8c --- /dev/null +++ b/include/epiworld/models/seirentitiesconnected.hpp @@ -0,0 +1,390 @@ +#ifndef EPIWORLD_MODELS_SEIRENTITIESCONNECTED_HPP +#define EPIWORLD_MODELS_SEIRENTITIESCONNECTED_HPP + +/** + * @file seirentitiesconnected.hpp + * @brief Template for a Susceptible-Exposed-Infected-Removed (SEIR) model with entities + */ +template +class ModelSEIREntitiesConn : public epiworld::Model +{ +public: + + static const int SUSCEPTIBLE = 0; + static const int EXPOSED = 1; + static const int INFECTED = 2; + static const int RECOVERED = 3; + + + ModelSEIREntitiesConn() {}; + + + /** + * @brief Constructs a ModelSEIREntitiesConn object. + * + * @param model A reference to an existing ModelSEIREntitiesConn object. + * @param vname The name of the ModelSEIREntitiesConn object. + * @param n The number of entities in the model. + * @param prevalence The initial prevalence of the disease in the model. + * @param contact_rate The contact rate between entities in the model. + * @param transmission_rate The transmission rate of the disease in the model. + * @param avg_incubation_days The average incubation period of the disease in the model. + * @param recovery_rate The recovery rate of the disease in the model. + * @param entities A vector of entity values. + * @param entities_names A vector of entity names. + */ + ModelSEIREntitiesConn( + ModelSEIREntitiesConn & model, + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double avg_incubation_days, + epiworld_double recovery_rate, + std::vector< epiworld_double > entities, + std::vector< std::string > entities_names + ); + + /** + * @brief Constructs a ModelSEIREntitiesConn object. + * + * @param vname The name of the ModelSEIREntitiesConn object. + * @param n The number of entities in the model. + * @param prevalence The initial prevalence of the disease in the model. + * @param contact_rate The contact rate between entities in the model. + * @param transmission_rate The transmission rate of the disease in the model. + * @param avg_incubation_days The average incubation period of the disease in the model. + * @param recovery_rate The recovery rate of the disease in the model. + * @param entities A vector of entity values. + */ + ModelSEIREntitiesConn( + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double avg_incubation_days, + epiworld_double recovery_rate, + std::vector< epiworld_double > entities, + std::vector< std::string > entities_names + ); + + ModelSEIREntitiesConn & run( + epiworld_fast_uint ndays, + int seed = -1 + ); + + void reset(); + + Model * clone_ptr(); + + /** + * @brief Set the initial states of the model + * @param proportions_ Double vector with a single element: + * - The proportion of non-infected individuals who have recovered. + */ + ModelSEIREntitiesConn & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); + +}; + +// Global event that moves agents between states + + +template +inline ModelSEIREntitiesConn & ModelSEIREntitiesConn::run( + epiworld_fast_uint ndays, + int seed +) +{ + + Model::run(ndays, seed); + + return *this; + +} + +template +inline void ModelSEIREntitiesConn::reset() +{ + + Model::reset(); + + Model::set_rand_binom( + Model::size(), + static_cast( + Model::par("Contact rate"))/ + static_cast(Model::size()) + ); + + return; + +} + +template +inline Model * ModelSEIREntitiesConn::clone_ptr() +{ + + ModelSEIREntitiesConn * ptr = new ModelSEIREntitiesConn( + *dynamic_cast*>(this) + ); + + return dynamic_cast< Model *>(ptr); + +} + + +/** + * @brief Template for a Susceptible-Exposed-Infected-Removed (SEIR) model + * + * @param model A Model object where to set up the SIR. + * @param vname std::string Name of the virus + * @param prevalence Initial prevalence (proportion) + * @param contact_rate Average number of contacts (interactions) per step. + * @param transmission_rate Probability of transmission + * @param recovery_rate Probability of recovery + */ +template +inline ModelSEIREntitiesConn::ModelSEIREntitiesConn( + ModelSEIREntitiesConn & model, + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double avg_incubation_days, + epiworld_double recovery_rate, + std::vector< epiworld_double > entities, + std::vector< std::string > entities_names + ) +{ + + // Instantiating entities + int entity_id = 0; + for (auto & e : entities) + model.add_entity_n( + epiworld::Entity(entities_names[entity_id++]), + e * n + ); + + epiworld::UpdateFun update_susceptible = []( + epiworld::Agent * p, epiworld::Model * m + ) -> void + { + + // Sampling how many individuals + int ndraw = m->rbinom(); + + if (ndraw == 0) + return; + + // Sampling from the agent's entities + auto sample = epiworld::AgentsSample(m, *p, ndraw, {}, true); + + // Drawing from the set + int nviruses_tmp = 0; + for (const auto & neighbor: sample) + { + + // Can't sample itself + if (neighbor->get_id() == static_cast(p->get_id())) + continue; + + // If the neighbor is infected, then proceed + if (neighbor->get_state() == ModelSEIREntitiesConn::INFECTED) + { + + 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 + if (nviruses_tmp == 0u) + return; + + // Running the roulette + int which = roulette(nviruses_tmp, m); + + if (which < 0) + return; + + p->set_virus( + *m->array_virus_tmp[which], + m, + ModelSEIREntitiesConn::EXPOSED + ); + + return; + + }; + + epiworld::UpdateFun update_infected = []( + epiworld::Agent * p, epiworld::Model * m + ) -> void { + + auto state = p->get_state(); + + if (state == ModelSEIREntitiesConn::EXPOSED) + { + + // Getting the virus + auto & v = p->get_virus(); + + // Does the agent become infected? + if (m->runif() < 1.0/(v->get_incubation(m))) + { + + p->change_state(m, ModelSEIREntitiesConn::INFECTED); + return; + + } + + + } else if (state == ModelSEIREntitiesConn::INFECTED) + { + + + // Odd: Die, Even: Recover + epiworld_fast_uint n_events = 0u; + const auto & v = p->get_virus(); + + // Recover + m->array_double_tmp[n_events++] = + 1.0 - (1.0 - v->get_prob_recovery(m)) * (1.0 - p->get_recovery_enhancer(v, m)); + + #ifdef EPI_DEBUG + if (n_events == 0u) + { + printf_epiworld( + "[epi-debug] agent %i has 0 possible events!!\n", + static_cast(p->get_id()) + ); + throw std::logic_error("Zero events in exposed."); + } + #else + if (n_events == 0u) + return; + #endif + + + // Running the roulette + int which = roulette(n_events, m); + + if (which < 0) + return; + + // Which roulette happen? + p->rm_virus(m); + + return ; + + } else + throw std::logic_error("This function can only be applied to exposed or infected individuals. (SEIR)") ; + + return; + + }; + + // Setting up parameters + model.add_param(contact_rate, "Contact rate"); + model.add_param(transmission_rate, "Prob. Transmission"); + model.add_param(recovery_rate, "Prob. Recovery"); + model.add_param(avg_incubation_days, "Avg. Incubation days"); + + // state + model.add_state("Susceptible", update_susceptible); + model.add_state("Exposed", update_infected); + model.add_state("Infected", update_infected); + model.add_state("Recovered"); + + + // Preparing the virus ------------------------------------------- + epiworld::Virus virus(vname); + virus.set_state( + ModelSEIREntitiesConn::EXPOSED, + ModelSEIREntitiesConn::RECOVERED, + ModelSEIREntitiesConn::RECOVERED + ); + + virus.set_prob_infecting(&model("Prob. Transmission")); + virus.set_prob_recovery(&model("Prob. Recovery")); + virus.set_incubation(&model("Avg. Incubation days")); + + model.add_virus(virus, prevalence); + + model.queuing_off(); // No queuing need + + // Adding the empty population + model.agents_empty_graph(n); + + model.set_name("Susceptible-Exposed-Infected-Removed (SEIR) (connected)"); + + return; + +} + +template +inline ModelSEIREntitiesConn::ModelSEIREntitiesConn( + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double avg_incubation_days, + epiworld_double recovery_rate, + std::vector< epiworld_double > entities, + std::vector< std::string > entity_names + ) +{ + + ModelSEIREntitiesConn( + *this, + vname, + n, + prevalence, + contact_rate, + transmission_rate, + avg_incubation_days, + recovery_rate, + entities, + entity_names + ); + + return; + +} + +template +inline ModelSEIREntitiesConn & ModelSEIREntitiesConn::initial_states( + std::vector< double > proportions_, + std::vector< int > /* queue_ */ +) +{ + + Model::initial_states_fun = + create_init_function_seir(proportions_) + ; + + return *this; + +} + +#endif