diff --git a/.gitignore b/.gitignore index 37b16013..fbadc30d 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,6 @@ examples/*/* !*.cpp .Rproj.user .Rhistory -tests +tests/* +!tests/ +!tests/*.cpp diff --git a/.vscode/settings.json b/.vscode/settings.json index 639f0a55..4976af68 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -70,5 +70,11 @@ "cfenv": "cpp", "cinttypes": "cpp" }, - "intel-corporation.oneapi-analysis-configurator.binary-path": "/home/george/Documents/development/epiworld/tests/01c.o" + "intel-corporation.oneapi-analysis-configurator.binary-path": "/home/george/Documents/development/epiworld/tests/01c.o", + "grammarly.selectors": [ + { + "language": "quarto", + "scheme": "file" + } + ] } \ No newline at end of file diff --git a/epiworld.hpp b/epiworld.hpp index e18df010..927837e4 100644 --- a/epiworld.hpp +++ b/epiworld.hpp @@ -19,7 +19,7 @@ /* Versioning */ #define EPIWORLD_VERSION_MAJOR 0 -#define EPIWORLD_VERSION_MINOR 1 +#define EPIWORLD_VERSION_MINOR 2 #define EPIWORLD_VERSION_PATCH 0 static const int epiworld_version_major = EPIWORLD_VERSION_MAJOR; @@ -6124,6 +6124,9 @@ class Model { std::vector< ToolToAgentFun > tools_dist_funs = {}; std::vector< Entity > entities = {}; + std::vector< epiworld_double > prevalence_entity = {}; + std::vector< bool > prevalence_entity_as_proportion = {}; + std::vector< EntityToAgentFun > entities_dist_funs = {}; std::vector< Entity > entities_backup = {}; std::mt19937 engine; @@ -6162,7 +6165,7 @@ class Model { void dist_tools(); void dist_virus(); - // void dist_entities(); + void dist_entities(); std::chrono::time_point time_start; std::chrono::time_point time_end; @@ -6237,6 +6240,7 @@ class Model { std::vector array_double_tmp; std::vector * > array_virus_tmp; + std::vector< int > array_int_tmp; Model(); Model(const Model & m); @@ -6323,6 +6327,8 @@ class Model { void add_tool_n(Tool & t, epiworld_fast_uint preval); void add_tool_fun(Tool & t, ToolToAgentFun fun); void add_entity(Entity e); + void add_entity_n(Entity e, epiworld_fast_uint preval); + void add_entity_fun(Entity e, EntityToAgentFun fun); void rm_virus(size_t virus_pos); void rm_tool(size_t tool_pos); void rm_entity(size_t entity_pos); @@ -6339,6 +6345,14 @@ class Model { * @param skip How many rows to skip. */ void load_agents_entities_ties(std::string fn, int skip); + + /** + * @brief Associate agents-entities from data + */ + void load_agents_entities_ties( + const std::vector & agents_ids, + const std::vector & entities_ids + ); /** * @name Accessing population of the model @@ -7115,10 +7129,10 @@ inline Model::Model(const Model & model) : prevalence_tool_as_proportion(model.prevalence_tool_as_proportion), tools_dist_funs(model.tools_dist_funs), entities(model.entities), + prevalence_entity(model.prevalence_entity), + prevalence_entity_as_proportion(model.prevalence_entity_as_proportion), + entities_dist_funs(model.entities_dist_funs), entities_backup(model.entities_backup), - // prevalence_entity(model.prevalence_entity), - // prevalence_entity_as_proportion(model.prevalence_entity_as_proportion), - // entities_dist_funs(model.entities_dist_funs), rewire_fun(model.rewire_fun), rewire_prop(model.rewire_prop), parameters(model.parameters), @@ -7134,7 +7148,8 @@ inline Model::Model(const Model & model) : queue(model.queue), use_queuing(model.use_queuing), array_double_tmp(model.array_double_tmp.size()), - array_virus_tmp(model.array_virus_tmp.size()) + array_virus_tmp(model.array_virus_tmp.size()), + array_int_tmp(model.array_int_tmp.size()) { @@ -7191,10 +7206,10 @@ inline Model::Model(Model && model) : tools_dist_funs(std::move(model.tools_dist_funs)), // Entities entities(std::move(model.entities)), + prevalence_entity(std::move(model.prevalence_entity)), + prevalence_entity_as_proportion(std::move(model.prevalence_entity_as_proportion)), + entities_dist_funs(std::move(model.entities_dist_funs)), entities_backup(std::move(model.entities_backup)), - // prevalence_entity(std::move(model.prevalence_entity)), - // prevalence_entity_as_proportion(std::move(model.prevalence_entity_as_proportion)), - // entities_dist_funs(std::move(model.entities_dist_funs)), // Pseudo-RNG engine(std::move(model.engine)), runifd(std::move(model.runifd)), @@ -7219,7 +7234,8 @@ inline Model::Model(Model && model) : queue(std::move(model.queue)), use_queuing(model.use_queuing), array_double_tmp(model.array_double_tmp.size()), - array_virus_tmp(model.array_virus_tmp.size()) + array_virus_tmp(model.array_virus_tmp.size()), + array_int_tmp(model.array_int_tmp.size()) { db.model = this; @@ -7269,10 +7285,10 @@ inline Model & Model::operator=(const Model & m) tools_dist_funs = m.tools_dist_funs; entities = m.entities; + prevalence_entity = m.prevalence_entity; + prevalence_entity_as_proportion = m.prevalence_entity_as_proportion; + entities_dist_funs = m.entities_dist_funs; entities_backup = m.entities_backup; - // prevalence_entity = m.prevalence_entity; - // prevalence_entity_as_proportion = m.prevalence_entity_as_proportion; - // entities_dist_funs = m.entities_dist_funs; rewire_fun = m.rewire_fun; rewire_prop = m.rewire_prop; @@ -7313,6 +7329,7 @@ inline Model & Model::operator=(const Model & m) )); array_virus_tmp.resize(1024u); + array_int_tmp.resize(1024u * 1024); return *this; @@ -7592,62 +7609,62 @@ inline void Model::dist_tools() } -// template -// inline void Model::dist_entities() -// { +template +inline void Model::dist_entities() +{ + + // Starting first infection + int n = size(); + std::vector< size_t > idx(n); + for (epiworld_fast_uint e = 0; e < entities.size(); ++e) + { + + if (entities_dist_funs[e]) + { + + entities_dist_funs[e](entities[e], this); -// // Starting first infection -// int n = size(); -// std::vector< size_t > idx(n); -// for (epiworld_fast_uint e = 0; e < entities.size(); ++e) -// { - -// if (entities_dist_funs[e]) -// { - -// entities_dist_funs[e](entities[e], this); - -// } else { - -// // Picking how many -// int nsampled; -// if (prevalence_entity_as_proportion[e]) -// { -// nsampled = static_cast(std::floor(prevalence_entity[e] * size())); -// } -// else -// { -// nsampled = static_cast(prevalence_entity[e]); -// } - -// if (nsampled > static_cast(size())) -// throw std::range_error("There are only " + std::to_string(size()) + -// " individuals in the population. Cannot add the entity to " + std::to_string(nsampled)); + } else { + + // Picking how many + int nsampled; + if (prevalence_entity_as_proportion[e]) + { + nsampled = static_cast(std::floor(prevalence_entity[e] * size())); + } + else + { + nsampled = static_cast(prevalence_entity[e]); + } + + if (nsampled > static_cast(size())) + throw std::range_error("There are only " + std::to_string(size()) + + " individuals in the population. Cannot add the entity to " + std::to_string(nsampled)); -// Entity & entity = entities[e]; + Entity & entity = entities[e]; -// int n_left = n; -// std::iota(idx.begin(), idx.end(), 0); -// while (nsampled > 0) -// { -// int loc = static_cast(floor(runif() * n_left--)); + int n_left = n; + std::iota(idx.begin(), idx.end(), 0); + while (nsampled > 0) + { + int loc = static_cast(floor(runif() * n_left--)); -// population[idx[loc]].add_entity(entity, this, entity.state_init, entity.queue_init); + population[idx[loc]].add_entity(entity, this, entity.state_init, entity.queue_init); -// nsampled--; + nsampled--; -// std::swap(idx[loc], idx[n_left]); + std::swap(idx[loc], idx[n_left]); -// } + } -// } + } -// // Apply the events -// events_run(); + // Apply the events + events_run(); -// } + } -// } +} template inline void Model::chrono_start() { @@ -7928,6 +7945,32 @@ inline void Model::add_entity(Entity e) } +template +inline void Model::add_entity_n(Entity e, epiworld_fast_uint preval) +{ + + e.model = this; + e.id = entities.size(); + entities.push_back(e); + prevalence_entity.push_back(preval); + prevalence_entity_as_proportion.push_back(false); + entities_dist_funs.push_back(nullptr); + +} + +template +inline void Model::add_entity_fun(Entity e, EntityToAgentFun fun) +{ + + e.model = this; + e.id = entities.size(); + entities.push_back(e); + prevalence_entity.push_back(0.0); + prevalence_entity_as_proportion.push_back(false); + entities_dist_funs.push_back(fun); + +} + template inline void Model::rm_virus(size_t virus_pos) { @@ -8018,7 +8061,6 @@ inline void Model::load_agents_entities_ties( throw std::logic_error("The file " + fn + " was not found."); int linenum = 0; - std::vector< epiworld_fast_uint > source_; std::vector< std::vector< epiworld_fast_uint > > target_(entities.size()); target_.reserve(1e5); @@ -8059,36 +8101,64 @@ inline void Model::load_agents_entities_ties( } - // // Iterating over entities - // for (size_t e = 0u; e < entities.size(); ++e) - // { + return; + +} + +template +inline void Model::load_agents_entities_ties( + const std::vector< int > & agents_ids, + const std::vector< int > & entities_ids +) { - // // This entity will have individuals assigned to it, so we add it - // if (target_[e].size() > 0u) - // { + if (agents_ids.size() != entities_ids.size()) + throw std::length_error( + std::string("agents_ids (") + + std::to_string(agents_ids.size()) + + std::string(") and entities_ids (") + + std::to_string(entities_ids.size()) + + std::string(") should match.") + ); - // // Filling in the gaps - // prevalence_entity[e] = static_cast(target_[e].size()); - // prevalence_entity_as_proportion[e] = false; - // // Generating the assignment function - // auto who = target_[e]; - // entities_dist_funs[e] = - // [who](Entity & e, Model* m) -> void { + size_t n_entries = agents_ids.size(); + for (size_t i = 0u; i < n_entries; ++i) + { - // for (auto w : who) - // m->population[w].add_entity(e, m, e.state_init, e.queue_init); - - // return; - - // }; + if (agents_ids[i] >= this->population.size()) + throw std::length_error( + std::string("agents_ids[") + + std::to_string(i) + + std::string("] = ") + + std::to_string(agents_ids[i]) + + std::string(" is out of range (population size: ") + + std::to_string(this->population.size()) + + std::string(").") + ); + + + if (entities_ids[i] >= this->entities.size()) + throw std::length_error( + std::string("entities_ids[") + + std::to_string(i) + + std::string("] = ") + + std::to_string(entities_ids[i]) + + std::string(" is out of range (entities size: ") + + std::to_string(this->entities.size()) + + std::string(").") + ); - // } + // Adding the entity to the agent + this->population[agents_ids[i]].add_entity( + this->entities[entities_ids[i]], + nullptr /* Immediately add it to the agent */ + ); - // } + } return; + } template @@ -8215,6 +8285,7 @@ inline Model & Model::run( array_virus_tmp.resize(1024); + array_int_tmp.resize(1024 * 1024); // Checking whether the proposed state in/out/removed // are valid @@ -8799,11 +8870,9 @@ inline void Model::reset() { #endif } - else - { - for (auto & e: entities) - e.reset(); - } + + for (auto & e: entities) + e.reset(); current_date = 0; @@ -8816,6 +8885,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); @@ -12076,7 +12146,7 @@ inline void Entity::rm_agent(size_t idx) " out of " + std::to_string(n_agents) ); - model->population[idx].rm_entity(*this); + model->get_agents()[agents[idx]].rm_entity(*this, model); return; } @@ -12137,7 +12207,10 @@ template inline Agent * Entity::operator[](size_t i) { if (n_agents <= i) - throw std::logic_error("There are not that many agents in this entity."); + throw std::logic_error( + "There are not that many agents in this entity. " + + std::to_string(n_agents) + " <= " + std::to_string(i) + ); return &model->get_agents()[i]; } @@ -12209,6 +12282,13 @@ inline void Entity::reset() sampled_agents_n = 0u; sampled_agents_left.clear(); sampled_agents_left_n = 0u; + + // Removing agents from entities + for (size_t i = 0u; i < n_agents; ++i) + this->rm_agent(i); + + return; + } template @@ -13160,9 +13240,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; @@ -13652,11 +13732,15 @@ inline void default_rm_entity(Event & a, Model * m) // When we move the end entity to the new location, the // moved entity needs to reflect the change, i.e., where the // entity will now be located in the agent - size_t agent_location_in_last_entity = p->entities_locations[p->n_entities]; - Entity * last_entity = &m->get_entities()[p->entities[p->n_entities]]; ///< Last entity of the agent + size_t agent_location_in_last_entity = + p->entities_locations[p->n_entities]; + + Entity * last_entity = + &m->get_entities()[p->entities[p->n_entities]]; ///< Last entity of the agent // The end entity will be located where the removed was - last_entity->agents_location[agent_location_in_last_entity] = idx_entity_in_agent; + last_entity->agents_location[agent_location_in_last_entity] = + idx_entity_in_agent; // We now make the swap std::swap( @@ -13673,10 +13757,13 @@ inline void default_rm_entity(Event & a, Model * m) // moved agent needs to reflect the change, i.e., where the // agent will now be located in the entity size_t entity_location_in_last_agent = e->agents_location[e->n_agents]; - Agent * last_agent = &m->get_agents()[e->agents[e->n_agents]]; ///< Last agent of the entity + + Agent * last_agent = + &m->get_agents()[e->agents[e->n_agents]]; ///< Last agent of the entity // The end entity will be located where the removed was - last_agent->entities_locations[entity_location_in_last_agent] = idx_agent_in_entity; + last_agent->entities_locations[entity_location_in_last_agent] = + idx_agent_in_entity; // We now make the swap std::swap( @@ -13935,7 +14022,7 @@ inline void Agent::add_entity( -1, -1 ); - // default_add_entity(a, model); /* passing model makes nothing */ + default_add_entity(a, model); /* passing model makes nothing */ } @@ -14029,8 +14116,15 @@ inline void Agent::rm_entity( CHECK_COALESCE_(queue, model->entities[entity_idx].queue_post, Queue::NoOne); model->events_add( - this, nullptr, nullptr, model->entities[entity_idx], state_new, queue, - default_rm_entity, entities_locations[entity_idx], entity_idx + this, + nullptr, + nullptr, + &model->entities[entities[entity_idx]], + state_new, + queue, + default_rm_entity, + entities_locations[entity_idx], + entity_idx ); } @@ -14047,22 +14141,35 @@ inline void Agent::rm_entity( int entity_idx = -1; for (size_t i = 0u; i < n_entities; ++i) { - if (entities[i] == entity->get_id()) + if (static_cast(entities[i]) == entity.get_id()) + { entity_idx = i; + break; + } } if (entity_idx == -1) throw std::logic_error( - "The agent " + std::to_string(id) + " is not associated with entity \"" + - entity.get_name() + "\"." + std::string("The agent ") + + std::to_string(id) + + std::string(" is not associated with entity \"") + + entity.get_name() + + std::string("\".") ); CHECK_COALESCE_(state_new, entity.state_post, state); CHECK_COALESCE_(queue, entity.queue_post, Queue::NoOne); model->events_add( - this, nullptr, nullptr, entities[entity_idx], state_new, queue, - default_rm_entity, entities_locations[entity_idx], entity_idx + this, + nullptr, + nullptr, + &model->entities[entity.get_id()], + state_new, + queue, + default_rm_entity, + entities_locations[entity_idx], + entity_idx ); } @@ -14635,10 +14742,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. @@ -14654,7 +14761,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 @@ -14665,13 +14772,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 ); @@ -14792,9 +14903,6 @@ inline AgentsSample::AgentsSample( agents = &agent_.sampled_agents; agents_n = &agent_.sampled_agents_n; - agents_left = &agent_.sampled_agents_left; - agents_left_n = &agent_.sampled_agents_left_n; - // Computing the cumulative sum of counts across entities size_t agents_in_entities = 0; Entities entities_a = agent->get_entities(); @@ -14833,31 +14941,37 @@ 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(); // Checking if states was specified if (states.size()) { + + // Getting the state + size_t state = model->population[agent_idx].get_state(); + if (std::find(states.begin(), states.end(), state) != states.end()) continue; + } - agents->operator[](i_obs++) = agent_idx; + agents->operator[](i_obs++) = &(model->population[agent_idx]); break; } @@ -15077,6 +15191,186 @@ inline void AgentsSample::sample_n(size_t n) +/*////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + Start of -include/epiworld/groupsampler-bones.hpp- + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////*/ + + +#ifndef GROUPSAMPLER_BONES_HPP +#define GROUPSAMPLER_BONES_HPP + +/** + * @brief Weighted sampling of groups + */ +template +class GroupSampler { + +private: + + std::vector< double > contact_matrix; ///< Contact matrix between groups + std::vector< size_t > group_sizes; ///< Sizes of the groups + std::vector< double > cumulate; ///< Cumulative sum of the contact matrix (row-major for faster access) + + /** + * @brief Get the index of the contact matrix + * + * The matrix is a vector stored in column-major order. + * + * @param i Index of the row + * @param j Index of the column + * @return Index of the contact matrix + */ + inline int idx(const int i, const int j, bool rowmajor = false) const + { + + if (rowmajor) + return i * group_sizes.size() + j; + + return j * group_sizes.size() + i; + + } + +public: + + GroupSampler() {}; + + GroupSampler( + const std::vector< double > & contact_matrix_, + const std::vector< size_t > & group_sizes_, + bool normalize = true + ); + + int sample_1( + Model * model, + const int origin_group + ); + + void sample_n( + Model * model, + std::vector< int > & sample, + const int origin_group, + const int nsamples + ); + +}; + +#endif +/*////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + End of -include/epiworld/groupsampler-bones.hpp- + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////*/ + + +/*////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + Start of -include/epiworld/groupsampler-meat.hpp- + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////*/ + + +#ifndef GROUPSAMPLER_MEAT_HPP +#define GROUPSAMPLER_MEAT_HPP + +template +inline GroupSampler::GroupSampler( + const std::vector< double > & contact_matrix_, + const std::vector< size_t > & group_sizes_, + bool normalize + ): contact_matrix(contact_matrix_), group_sizes(group_sizes_) { + + + this->cumulate.resize(contact_matrix.size()); + std::fill(cumulate.begin(), cumulate.end(), 0.0); + + // Cumulative sum + for (size_t j = 0; j < group_sizes.size(); ++j) + { + for (size_t i = 0; i < group_sizes.size(); ++i) + cumulate[idx(i, j, true)] += + cumulate[idx(i, j - 1, true)] + + contact_matrix[idx(i, j)]; + } + + if (normalize) + { + for (size_t i = 0; i < group_sizes.size(); ++i) + { + double sum = 0.0; + for (size_t j = 0; j < group_sizes.size(); ++j) + sum += contact_matrix[idx(i, j, true)]; + for (size_t j = 0; j < group_sizes.size(); ++j) + contact_matrix[idx(i, j, true)] /= sum; + } + } + + }; + +template +int GroupSampler::sample_1( + Model * model, + const int origin_group + ) +{ + + // Random number + double r = model->runif(); + + // Finding the group + size_t j = 0; + while (r > cumulate[idx(origin_group, j, true)]) + ++j; + + // Adjusting the prob + r = r - (j == 0 ? 0.0 : cumulate[idx(origin_group, j - 1, true)]); + + int res = static_cast( + std::floor(r * group_sizes[j]) + ); + + // Making sure we are not picking outside of the group + if (res >= static_cast(group_sizes[j])) + res = static_cast(group_sizes[j]) - 1; + + return model->get_entities()[j][res]->get_id(); + +} + +template +void GroupSampler::sample_n( + Model * model, + std::vector< int > & sample, + const int origin_group, + const int nsamples +) +{ + + for (int i = 0; i < nsamples; ++i) + sample[i] = sample_1(model, origin_group); + + return; + +} + +#endif +/*////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + End of -include/epiworld/groupsampler-meat.hpp- + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////*/ + + + /*////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// @@ -19137,6 +19431,1055 @@ inline ModelDiffNet::ModelDiffNet( //////////////////////////////////////////////////////////////////////////////*/ +/*////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + Start of -include/epiworld//models/seirmixing.hpp- + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////*/ + + +#ifndef EPIWORLD_MODELS_SEIRMIXING_HPP +#define EPIWORLD_MODELS_SEIRMIXING_HPP + +/** + * @file seirentitiesconnected.hpp + * @brief Template for a Susceptible-Exposed-Infected-Removed (SEIR) model with mixing + */ +template +class ModelSEIRMixing : public epiworld::Model +{ +private: + std::vector< std::vector< epiworld::Agent * > > infected; + void update_infected(); + std::vector< epiworld::Agent * > sampled_agents; + size_t sample_agents( + epiworld::Agent * agent, + std::vector< epiworld::Agent * > & sampled_agents + ); + double adjusted_contact_rate; + std::vector< double > contact_matrix; + + size_t index(size_t i, size_t j, size_t n) { + return j * n + i; + } + +public: + + static const int SUSCEPTIBLE = 0; + static const int EXPOSED = 1; + static const int INFECTED = 2; + static const int RECOVERED = 3; + + ModelSEIRMixing() {}; + + /** + * @brief Constructs a ModelSEIRMixing object. + * + * @param model A reference to an existing ModelSEIRMixing object. + * @param vname The name of the ModelSEIRMixing 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 contact_matrix The contact matrix between entities in the model. + */ + ModelSEIRMixing( + ModelSEIRMixing & 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< double > contact_matrix + ); + + /** + * @brief Constructs a ModelSEIRMixing object. + * + * @param vname The name of the ModelSEIRMixing 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 contact_matrix The contact matrix between entities in the model. + */ + ModelSEIRMixing( + 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< double > contact_matrix + ); + + ModelSEIRMixing & 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. + */ + ModelSEIRMixing & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); + + size_t get_n_infected(size_t group) const + { + return infected[group].size(); + } + + void set_contact_matrix(std::vector< double > cmat) + { + contact_matrix = cmat; + return; + }; + +}; + +template +inline void ModelSEIRMixing::update_infected() +{ + + auto & agents = Model::get_agents(); + auto & entities = Model::get_entities(); + + infected.resize(entities.size()); + sampled_agents.resize(agents.size()); + + // Checking contact matrix's rows add to one + size_t nentities = entities.size(); + if (this->contact_matrix.size() != nentities*nentities) + throw std::length_error( + std::string("The contact matrix must be a square matrix of size ") + + std::string("nentities x nentities. ") + + std::to_string(this->contact_matrix.size()) + + std::string(" != ") + std::to_string(nentities*nentities) + + std::string(".") + ); + + for (size_t i = 0u; i < entities.size(); ++i) + { + double sum = 0.0; + for (size_t j = 0u; j < entities.size(); ++j) + { + if (this->contact_matrix[index(i, j, nentities)] < 0.0) + throw std::range_error( + std::string("The contact matrix must be non-negative. ") + + std::to_string(this->contact_matrix[index(i, j, nentities)]) + + std::string(" < 0.") + ); + sum += this->contact_matrix[index(i, j, nentities)]; + } + if (sum < 0.999 || sum > 1.001) + throw std::range_error( + std::string("The contact matrix must have rows that add to one. ") + + std::to_string(sum) + + std::string(" != 1.") + ); + } + + for (size_t i = 0; i < entities.size(); ++i) + { + infected[i].clear(); + infected[i].reserve(agents.size()); + } + + for (auto & a : agents) + { + + if (a.get_state() == ModelSEIRMixing::INFECTED) + infected[a.get_entity(0u).get_id()].push_back(&a); + + } + + // Adjusting contact rate + adjusted_contact_rate = Model::get_param("Contact rate") / + agents.size(); + + return; + +} + +template +inline size_t ModelSEIRMixing::sample_agents( + epiworld::Agent * agent, + std::vector< epiworld::Agent * > & sampled_agents + ) +{ + + size_t agent_group_id = agent->get_entity(0u).get_id(); + size_t ngroups = infected.size(); + + int samp_id = 0; + for (size_t g = 0; g < infected.size(); ++g) + { + + // How many from this entity? + int nsamples = epiworld::Model::rbinom( + infected[g].size(), + adjusted_contact_rate * contact_matrix[ + index(agent_group_id, g, ngroups) + ] + ); + + if (nsamples == 0) + continue; + + // Sampling from the entity + for (int s = 0; s < nsamples; ++s) + { + + // Randomly selecting an agent + int which = epiworld::Model::runif() * infected[g].size(); + + // Correcting overflow error + if (which >= static_cast(infected[g].size())) + which = static_cast(infected[g].size()) - 1; + + auto & a = infected[g][which]; + + // Can't sample itself + if (a->get_id() == agent->get_id()) + continue; + + sampled_agents[samp_id++] = a; + + } + + } + + return samp_id; + +} + +template +inline ModelSEIRMixing & ModelSEIRMixing::run( + epiworld_fast_uint ndays, + int seed +) +{ + + Model::run(ndays, seed); + return *this; + +} + +template +inline void ModelSEIRMixing::reset() +{ + + Model::reset(); + this->update_infected(); + + return; + +} + +template +inline Model * ModelSEIRMixing::clone_ptr() +{ + + ModelSEIRMixing * ptr = new ModelSEIRMixing( + *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 ModelSEIRMixing::ModelSEIRMixing( + ModelSEIRMixing & 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< double > contact_matrix + ) +{ + + // Setting up the contact matrix + this->contact_matrix = contact_matrix; + + epiworld::UpdateFun update_susceptible = []( + epiworld::Agent * p, epiworld::Model * m + ) -> void + { + + if (p->get_n_entities() == 0) + return; + + // Downcasting to retrieve the sampler attached to the + // class + ModelSEIRMixing * m_down = + dynamic_cast *>(m); + + size_t ndraws = m_down->sample_agents(p, m_down->sampled_agents); + + if (ndraws == 0u) + return; + + + // Drawing from the set + int nviruses_tmp = 0; + for (size_t n = 0u; n < ndraws; ++n) + { + + auto & neighbor = m_down->sampled_agents[n]; + + 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); + + } + + // Running the roulette + int which = roulette(nviruses_tmp, m); + + if (which < 0) + return; + + p->set_virus( + *m->array_virus_tmp[which], + m, + ModelSEIRMixing::EXPOSED + ); + + return; + + }; + + epiworld::UpdateFun update_infected = []( + epiworld::Agent * p, epiworld::Model * m + ) -> void { + + auto state = p->get_state(); + + if (state == ModelSEIRMixing::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, ModelSEIRMixing::INFECTED); + return; + + } + + + } else if (state == ModelSEIRMixing::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"); + + // Global function + epiworld::GlobalFun update = [](epiworld::Model * m) -> void + { + + ModelSEIRMixing * m_down = + dynamic_cast *>(m); + + m_down->update_infected(); + + return; + + }; + + model.add_globalevent(update, "Update infected individuals"); + + + // Preparing the virus ------------------------------------------- + epiworld::Virus virus(vname); + virus.set_state( + ModelSEIRMixing::EXPOSED, + ModelSEIRMixing::RECOVERED, + ModelSEIRMixing::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) with Mixing"); + + return; + +} + +template +inline ModelSEIRMixing::ModelSEIRMixing( + 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< double > contact_matrix + ) +{ + + this->contact_matrix = contact_matrix; + + ModelSEIRMixing( + *this, + vname, + n, + prevalence, + contact_rate, + transmission_rate, + avg_incubation_days, + recovery_rate, + contact_matrix + ); + + return; + +} + +template +inline ModelSEIRMixing & ModelSEIRMixing::initial_states( + std::vector< double > proportions_, + std::vector< int > /* queue_ */ +) +{ + + Model::initial_states_fun = + create_init_function_seir(proportions_) + ; + + return *this; + +} + +#endif +/*////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + End of -include/epiworld//models/seirmixing.hpp- + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////*/ + + +/*////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + Start of -include/epiworld//models/sirmixing.hpp- + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////*/ + + +#ifndef EPIWORLD_MODELS_SIRMIXING_HPP +#define EPIWORLD_MODELS_SIRMIXING_HPP + +/** + * @file seirentitiesconnected.hpp + * @brief Template for a Susceptible-Exposed-Infected-Removed (SEIR) model with mixing + */ +template +class ModelSIRMixing : public epiworld::Model +{ +private: + std::vector< std::vector< epiworld::Agent * > > infected; + void update_infected_list(); + std::vector< epiworld::Agent * > sampled_agents; + size_t sample_agents( + epiworld::Agent * agent, + std::vector< epiworld::Agent * > & sampled_agents + ); + double adjusted_contact_rate; + std::vector< double > contact_matrix; + + size_t index(size_t i, size_t j, size_t n) { + return j * n + i; + } + +public: + + static const int SUSCEPTIBLE = 0; + static const int INFECTED = 1; + static const int RECOVERED = 2; + + ModelSIRMixing() {}; + + /** + * @brief Constructs a ModelSIRMixing object. + * + * @param model A reference to an existing ModelSIRMixing object. + * @param vname The name of the ModelSIRMixing 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 recovery_rate The recovery rate of the disease in the model. + * @param contact_matrix The contact matrix between entities in the model. + */ + ModelSIRMixing( + ModelSIRMixing & model, + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double recovery_rate, + std::vector< double > contact_matrix + ); + + /** + * @brief Constructs a ModelSIRMixing object. + * + * @param vname The name of the ModelSIRMixing 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 recovery_rate The recovery rate of the disease in the model. + * @param contact_matrix The contact matrix between entities in the model. + */ + ModelSIRMixing( + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double recovery_rate, + std::vector< double > contact_matrix + ); + + ModelSIRMixing & 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. + */ + ModelSIRMixing & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); + + size_t get_n_infected(size_t group) const + { + return infected[group].size(); + } + + void set_contact_matrix(std::vector< double > cmat) + { + contact_matrix = cmat; + return; + }; + +}; + +template +inline void ModelSIRMixing::update_infected_list() +{ + + auto & agents = Model::get_agents(); + auto & entities = Model::get_entities(); + + infected.resize(entities.size()); + sampled_agents.resize(agents.size()); + + // Checking contact matrix's rows add to one + size_t nentities = entities.size(); + if (this->contact_matrix.size() != nentities*nentities) + throw std::length_error( + std::string("The contact matrix must be a square matrix of size ") + + std::string("nentities x nentities. ") + + std::to_string(this->contact_matrix.size()) + + std::string(" != ") + std::to_string(nentities*nentities) + + std::string(".") + ); + + for (size_t i = 0u; i < entities.size(); ++i) + { + double sum = 0.0; + for (size_t j = 0u; j < entities.size(); ++j) + { + if (this->contact_matrix[index(i, j, nentities)] < 0.0) + throw std::range_error( + std::string("The contact matrix must be non-negative. ") + + std::to_string(this->contact_matrix[index(i, j, nentities)]) + + std::string(" < 0.") + ); + sum += this->contact_matrix[index(i, j, nentities)]; + } + if (sum < 0.999 || sum > 1.001) + throw std::range_error( + std::string("The contact matrix must have rows that add to one. ") + + std::to_string(sum) + + std::string(" != 1.") + ); + } + + for (size_t i = 0; i < entities.size(); ++i) + { + infected[i].clear(); + infected[i].reserve(agents.size()); + } + + for (auto & a : agents) + { + + if (a.get_state() == ModelSIRMixing::INFECTED) + infected[a.get_entity(0u).get_id()].push_back(&a); + + } + + // Adjusting contact rate + adjusted_contact_rate = Model::get_param("Contact rate") / + agents.size(); + + return; + +} + +template +inline size_t ModelSIRMixing::sample_agents( + epiworld::Agent * agent, + std::vector< epiworld::Agent * > & sampled_agents + ) +{ + + size_t agent_group_id = agent->get_entity(0u).get_id(); + size_t ngroups = infected.size(); + + int samp_id = 0; + for (size_t g = 0; g < infected.size(); ++g) + { + + // How many from this entity? + int nsamples = epiworld::Model::rbinom( + infected[g].size(), + adjusted_contact_rate * contact_matrix[ + index(agent_group_id, g, ngroups) + ] + ); + + if (nsamples == 0) + continue; + + // Sampling from the entity + for (int s = 0; s < nsamples; ++s) + { + + // Randomly selecting an agent + int which = epiworld::Model::runif() * infected[g].size(); + + // Correcting overflow error + if (which >= static_cast(infected[g].size())) + which = static_cast(infected[g].size()) - 1; + + auto & a = infected[g][which]; + + // Can't sample itself + if (a->get_id() == agent->get_id()) + continue; + + sampled_agents[samp_id++] = a; + + } + + } + + return samp_id; + +} + +template +inline ModelSIRMixing & ModelSIRMixing::run( + epiworld_fast_uint ndays, + int seed +) +{ + + Model::run(ndays, seed); + return *this; + +} + +template +inline void ModelSIRMixing::reset() +{ + + Model::reset(); + this->update_infected_list(); + + return; + +} + +template +inline Model * ModelSIRMixing::clone_ptr() +{ + + ModelSIRMixing * ptr = new ModelSIRMixing( + *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 ModelSIRMixing::ModelSIRMixing( + ModelSIRMixing & model, + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double recovery_rate, + std::vector< double > contact_matrix + ) +{ + + // Setting up the contact matrix + this->contact_matrix = contact_matrix; + + epiworld::UpdateFun update_susceptible = []( + epiworld::Agent * p, epiworld::Model * m + ) -> void + { + + if (p->get_n_entities() == 0) + return; + + // Downcasting to retrieve the sampler attached to the + // class + ModelSIRMixing * m_down = + dynamic_cast *>(m); + + size_t ndraws = m_down->sample_agents(p, m_down->sampled_agents); + + if (ndraws == 0u) + return; + + + // Drawing from the set + int nviruses_tmp = 0; + for (size_t n = 0u; n < ndraws; ++n) + { + + auto & neighbor = m_down->sampled_agents[n]; + + 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); + + } + + // Running the roulette + int which = roulette(nviruses_tmp, m); + + if (which < 0) + return; + + p->set_virus( + *m->array_virus_tmp[which], + m, + ModelSIRMixing::INFECTED + ); + + return; + + }; + + epiworld::UpdateFun update_infected = []( + epiworld::Agent * p, epiworld::Model * m + ) -> void { + + auto state = p->get_state(); + + if (state == ModelSIRMixing::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 infected."); + } + #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 infected individuals. (SIR)") ; + + 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"); + + // state + model.add_state("Susceptible", update_susceptible); + model.add_state("Infected", update_infected); + model.add_state("Recovered"); + + // Global function + epiworld::GlobalFun update = [](epiworld::Model * m) -> void + { + + ModelSIRMixing * m_down = + dynamic_cast *>(m); + + m_down->update_infected_list(); + + return; + + }; + + model.add_globalevent(update, "Update infected individuals"); + + + // Preparing the virus ------------------------------------------- + epiworld::Virus virus(vname); + virus.set_state( + ModelSIRMixing::INFECTED, + ModelSIRMixing::RECOVERED, + ModelSIRMixing::RECOVERED + ); + + virus.set_prob_infecting(&model("Prob. Transmission")); + virus.set_prob_recovery(&model("Prob. Recovery")); + + model.add_virus(virus, prevalence); + + model.queuing_off(); // No queuing need + + // Adding the empty population + model.agents_empty_graph(n); + + model.set_name("Susceptible-Infected-Removed (SIR) with Mixing"); + + return; + +} + +template +inline ModelSIRMixing::ModelSIRMixing( + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double recovery_rate, + std::vector< double > contact_matrix + ) +{ + + this->contact_matrix = contact_matrix; + + ModelSIRMixing( + *this, + vname, + n, + prevalence, + contact_rate, + transmission_rate, + recovery_rate, + contact_matrix + ); + + return; + +} + +template +inline ModelSIRMixing & ModelSIRMixing::initial_states( + std::vector< double > proportions_, + std::vector< int > /* queue_ */ +) +{ + + Model::initial_states_fun = + create_init_function_sir(proportions_) + ; + + return *this; + +} + +#endif +/*////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + End of -include/epiworld//models/sirmixing.hpp- + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////*/ + + } diff --git a/examples/11-entities/Makefile b/examples/11-entities/Makefile new file mode 100644 index 00000000..0c4d9f21 --- /dev/null +++ b/examples/11-entities/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/11-entities/main.cpp b/examples/11-entities/main.cpp new file mode 100644 index 00000000..6db469bd --- /dev/null +++ b/examples/11-entities/main.cpp @@ -0,0 +1,55 @@ +// #define EPI_DEBUG +#include "../../include/epiworld/epiworld.hpp" + +using namespace epiworld; + +template +EntityToAgentFun dist_factory(int from, int to) { + return [from, to](Entity<> & e, Model<> * m) -> void { + + auto & agents = m->get_agents(); + for (size_t i = from; i < to; ++i) + { + e.add_agent(&agents[i], m); + } + + return; + + }; +} + +int main() { + + std::vector< double > contact_matrix = { + 0.9, 0.1, 0.1, + 0.05, 0.8, .2, + 0.05, 0.1, 0.7 + }; + + epimodels::ModelSEIRMixing<> model( + "Flu", // std::string vname, + 10000, // epiworld_fast_uint n, + 0.01,// epiworld_double prevalence, + 10.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, + contact_matrix + ); + + // Creating three groups + Entity<> e1("Entity 1"); + Entity<> e2("Entity 2"); + Entity<> e3("Entity 3"); + + model.add_entity_fun(e1, dist_factory<>(0, 3000)); + model.add_entity_fun(e2, dist_factory<>(3000, 6000)); + model.add_entity_fun(e3, dist_factory<>(6000, 10000)); + + // 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/agent-events-meat.hpp b/include/epiworld/agent-events-meat.hpp index 8da043fb..275896e6 100644 --- a/include/epiworld/agent-events-meat.hpp +++ b/include/epiworld/agent-events-meat.hpp @@ -253,11 +253,15 @@ inline void default_rm_entity(Event & a, Model * m) // When we move the end entity to the new location, the // moved entity needs to reflect the change, i.e., where the // entity will now be located in the agent - size_t agent_location_in_last_entity = p->entities_locations[p->n_entities]; - Entity * last_entity = &m->get_entities()[p->entities[p->n_entities]]; ///< Last entity of the agent + size_t agent_location_in_last_entity = + p->entities_locations[p->n_entities]; + + Entity * last_entity = + &m->get_entities()[p->entities[p->n_entities]]; ///< Last entity of the agent // The end entity will be located where the removed was - last_entity->agents_location[agent_location_in_last_entity] = idx_entity_in_agent; + last_entity->agents_location[agent_location_in_last_entity] = + idx_entity_in_agent; // We now make the swap std::swap( @@ -274,10 +278,13 @@ inline void default_rm_entity(Event & a, Model * m) // moved agent needs to reflect the change, i.e., where the // agent will now be located in the entity size_t entity_location_in_last_agent = e->agents_location[e->n_agents]; - Agent * last_agent = &m->get_agents()[e->agents[e->n_agents]]; ///< Last agent of the entity + + Agent * last_agent = + &m->get_agents()[e->agents[e->n_agents]]; ///< Last agent of the entity // The end entity will be located where the removed was - last_agent->entities_locations[entity_location_in_last_agent] = idx_agent_in_entity; + last_agent->entities_locations[entity_location_in_last_agent] = + idx_agent_in_entity; // We now make the swap std::swap( diff --git a/include/epiworld/agent-meat.hpp b/include/epiworld/agent-meat.hpp index a5006b5d..1d907ed8 100644 --- a/include/epiworld/agent-meat.hpp +++ b/include/epiworld/agent-meat.hpp @@ -242,7 +242,7 @@ inline void Agent::add_entity( -1, -1 ); - // default_add_entity(a, model); /* passing model makes nothing */ + default_add_entity(a, model); /* passing model makes nothing */ } @@ -336,8 +336,15 @@ inline void Agent::rm_entity( CHECK_COALESCE_(queue, model->entities[entity_idx].queue_post, Queue::NoOne); model->events_add( - this, nullptr, nullptr, model->entities[entity_idx], state_new, queue, - default_rm_entity, entities_locations[entity_idx], entity_idx + this, + nullptr, + nullptr, + &model->entities[entities[entity_idx]], + state_new, + queue, + default_rm_entity, + entities_locations[entity_idx], + entity_idx ); } @@ -354,22 +361,35 @@ inline void Agent::rm_entity( int entity_idx = -1; for (size_t i = 0u; i < n_entities; ++i) { - if (entities[i] == entity->get_id()) + if (static_cast(entities[i]) == entity.get_id()) + { entity_idx = i; + break; + } } if (entity_idx == -1) throw std::logic_error( - "The agent " + std::to_string(id) + " is not associated with entity \"" + - entity.get_name() + "\"." + std::string("The agent ") + + std::to_string(id) + + std::string(" is not associated with entity \"") + + entity.get_name() + + std::string("\".") ); CHECK_COALESCE_(state_new, entity.state_post, state); CHECK_COALESCE_(queue, entity.queue_post, Queue::NoOne); model->events_add( - this, nullptr, nullptr, entities[entity_idx], state_new, queue, - default_rm_entity, entities_locations[entity_idx], entity_idx + this, + nullptr, + nullptr, + &model->entities[entity.get_id()], + state_new, + queue, + default_rm_entity, + entities_locations[entity_idx], + entity_idx ); } diff --git a/include/epiworld/agentssample-bones.hpp b/include/epiworld/agentssample-bones.hpp index c8462c58..e0d1285a 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 ); @@ -187,9 +191,6 @@ inline AgentsSample::AgentsSample( agents = &agent_.sampled_agents; agents_n = &agent_.sampled_agents_n; - agents_left = &agent_.sampled_agents_left; - agents_left_n = &agent_.sampled_agents_left_n; - // Computing the cumulative sum of counts across entities size_t agents_in_entities = 0; Entities entities_a = agent->get_entities(); @@ -228,31 +229,37 @@ 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(); // Checking if states was specified if (states.size()) { + + // Getting the state + size_t state = model->population[agent_idx].get_state(); + if (std::find(states.begin(), states.end(), state) != states.end()) continue; + } - agents->operator[](i_obs++) = agent_idx; + agents->operator[](i_obs++) = &(model->population[agent_idx]); break; } diff --git a/include/epiworld/entity-meat.hpp b/include/epiworld/entity-meat.hpp index 4c5a30d4..25fdbad0 100644 --- a/include/epiworld/entity-meat.hpp +++ b/include/epiworld/entity-meat.hpp @@ -31,7 +31,7 @@ inline void Entity::rm_agent(size_t idx) " out of " + std::to_string(n_agents) ); - model->population[idx].rm_entity(*this); + model->get_agents()[agents[idx]].rm_entity(*this, model); return; } @@ -92,7 +92,10 @@ template inline Agent * Entity::operator[](size_t i) { if (n_agents <= i) - throw std::logic_error("There are not that many agents in this entity."); + throw std::logic_error( + "There are not that many agents in this entity. " + + std::to_string(n_agents) + " <= " + std::to_string(i) + ); return &model->get_agents()[i]; } @@ -164,6 +167,13 @@ inline void Entity::reset() sampled_agents_n = 0u; sampled_agents_left.clear(); sampled_agents_left_n = 0u; + + // Removing agents from entities + for (size_t i = 0u; i < n_agents; ++i) + this->rm_agent(i); + + return; + } template diff --git a/include/epiworld/epiworld.hpp b/include/epiworld/epiworld.hpp index 5d7d7bdc..5ed8afe7 100644 --- a/include/epiworld/epiworld.hpp +++ b/include/epiworld/epiworld.hpp @@ -19,8 +19,8 @@ /* Versioning */ #define EPIWORLD_VERSION_MAJOR 0 -#define EPIWORLD_VERSION_MINOR 1 -#define EPIWORLD_VERSION_PATCH 1 +#define EPIWORLD_VERSION_MINOR 2 +#define EPIWORLD_VERSION_PATCH 0 static const int epiworld_version_major = EPIWORLD_VERSION_MAJOR; static const int epiworld_version_minor = EPIWORLD_VERSION_MINOR; @@ -79,6 +79,9 @@ namespace epiworld { #include "agentssample-bones.hpp" + #include "groupsampler-bones.hpp" + #include "groupsampler-meat.hpp" + #include "models/models.hpp" } diff --git a/include/epiworld/groupsampler-bones.hpp b/include/epiworld/groupsampler-bones.hpp new file mode 100644 index 00000000..88ec690a --- /dev/null +++ b/include/epiworld/groupsampler-bones.hpp @@ -0,0 +1,59 @@ +#ifndef GROUPSAMPLER_BONES_HPP +#define GROUPSAMPLER_BONES_HPP + +/** + * @brief Weighted sampling of groups + */ +template +class GroupSampler { + +private: + + std::vector< double > contact_matrix; ///< Contact matrix between groups + std::vector< size_t > group_sizes; ///< Sizes of the groups + std::vector< double > cumulate; ///< Cumulative sum of the contact matrix (row-major for faster access) + + /** + * @brief Get the index of the contact matrix + * + * The matrix is a vector stored in column-major order. + * + * @param i Index of the row + * @param j Index of the column + * @return Index of the contact matrix + */ + inline int idx(const int i, const int j, bool rowmajor = false) const + { + + if (rowmajor) + return i * group_sizes.size() + j; + + return j * group_sizes.size() + i; + + } + +public: + + GroupSampler() {}; + + GroupSampler( + const std::vector< double > & contact_matrix_, + const std::vector< size_t > & group_sizes_, + bool normalize = true + ); + + int sample_1( + Model * model, + const int origin_group + ); + + void sample_n( + Model * model, + std::vector< int > & sample, + const int origin_group, + const int nsamples + ); + +}; + +#endif \ No newline at end of file diff --git a/include/epiworld/groupsampler-meat.hpp b/include/epiworld/groupsampler-meat.hpp new file mode 100644 index 00000000..6e960965 --- /dev/null +++ b/include/epiworld/groupsampler-meat.hpp @@ -0,0 +1,84 @@ +#ifndef GROUPSAMPLER_MEAT_HPP +#define GROUPSAMPLER_MEAT_HPP + +template +inline GroupSampler::GroupSampler( + const std::vector< double > & contact_matrix_, + const std::vector< size_t > & group_sizes_, + bool normalize + ): contact_matrix(contact_matrix_), group_sizes(group_sizes_) { + + + this->cumulate.resize(contact_matrix.size()); + std::fill(cumulate.begin(), cumulate.end(), 0.0); + + // Cumulative sum + for (size_t j = 0; j < group_sizes.size(); ++j) + { + for (size_t i = 0; i < group_sizes.size(); ++i) + cumulate[idx(i, j, true)] += + cumulate[idx(i, j - 1, true)] + + contact_matrix[idx(i, j)]; + } + + if (normalize) + { + for (size_t i = 0; i < group_sizes.size(); ++i) + { + double sum = 0.0; + for (size_t j = 0; j < group_sizes.size(); ++j) + sum += contact_matrix[idx(i, j, true)]; + for (size_t j = 0; j < group_sizes.size(); ++j) + contact_matrix[idx(i, j, true)] /= sum; + } + } + + }; + +template +int GroupSampler::sample_1( + Model * model, + const int origin_group + ) +{ + + // Random number + double r = model->runif(); + + // Finding the group + size_t j = 0; + while (r > cumulate[idx(origin_group, j, true)]) + ++j; + + // Adjusting the prob + r = r - (j == 0 ? 0.0 : cumulate[idx(origin_group, j - 1, true)]); + + int res = static_cast( + std::floor(r * group_sizes[j]) + ); + + // Making sure we are not picking outside of the group + if (res >= static_cast(group_sizes[j])) + res = static_cast(group_sizes[j]) - 1; + + return model->get_entities()[j][res]->get_id(); + +} + +template +void GroupSampler::sample_n( + Model * model, + std::vector< int > & sample, + const int origin_group, + const int nsamples +) +{ + + for (int i = 0; i < nsamples; ++i) + sample[i] = sample_1(model, origin_group); + + return; + +} + +#endif \ No newline at end of file diff --git a/include/epiworld/model-bones.hpp b/include/epiworld/model-bones.hpp index 3b51e2a8..7959c608 100644 --- a/include/epiworld/model-bones.hpp +++ b/include/epiworld/model-bones.hpp @@ -145,6 +145,9 @@ class Model { std::vector< ToolToAgentFun > tools_dist_funs = {}; std::vector< Entity > entities = {}; + std::vector< epiworld_double > prevalence_entity = {}; + std::vector< bool > prevalence_entity_as_proportion = {}; + std::vector< EntityToAgentFun > entities_dist_funs = {}; std::vector< Entity > entities_backup = {}; std::mt19937 engine; @@ -183,7 +186,7 @@ class Model { void dist_tools(); void dist_virus(); - // void dist_entities(); + void dist_entities(); std::chrono::time_point time_start; std::chrono::time_point time_end; @@ -258,6 +261,7 @@ class Model { std::vector array_double_tmp; std::vector * > array_virus_tmp; + std::vector< int > array_int_tmp; Model(); Model(const Model & m); @@ -344,6 +348,8 @@ class Model { void add_tool_n(Tool & t, epiworld_fast_uint preval); void add_tool_fun(Tool & t, ToolToAgentFun fun); void add_entity(Entity e); + void add_entity_n(Entity e, epiworld_fast_uint preval); + void add_entity_fun(Entity e, EntityToAgentFun fun); void rm_virus(size_t virus_pos); void rm_tool(size_t tool_pos); void rm_entity(size_t entity_pos); @@ -360,6 +366,14 @@ class Model { * @param skip How many rows to skip. */ void load_agents_entities_ties(std::string fn, int skip); + + /** + * @brief Associate agents-entities from data + */ + void load_agents_entities_ties( + const std::vector & agents_ids, + const std::vector & entities_ids + ); /** * @name Accessing population of the model diff --git a/include/epiworld/model-meat.hpp b/include/epiworld/model-meat.hpp index 63bdc3e2..37e93379 100644 --- a/include/epiworld/model-meat.hpp +++ b/include/epiworld/model-meat.hpp @@ -388,10 +388,10 @@ inline Model::Model(const Model & model) : prevalence_tool_as_proportion(model.prevalence_tool_as_proportion), tools_dist_funs(model.tools_dist_funs), entities(model.entities), + prevalence_entity(model.prevalence_entity), + prevalence_entity_as_proportion(model.prevalence_entity_as_proportion), + entities_dist_funs(model.entities_dist_funs), entities_backup(model.entities_backup), - // prevalence_entity(model.prevalence_entity), - // prevalence_entity_as_proportion(model.prevalence_entity_as_proportion), - // entities_dist_funs(model.entities_dist_funs), rewire_fun(model.rewire_fun), rewire_prop(model.rewire_prop), parameters(model.parameters), @@ -407,7 +407,8 @@ inline Model::Model(const Model & model) : queue(model.queue), use_queuing(model.use_queuing), array_double_tmp(model.array_double_tmp.size()), - array_virus_tmp(model.array_virus_tmp.size()) + array_virus_tmp(model.array_virus_tmp.size()), + array_int_tmp(model.array_int_tmp.size()) { @@ -464,10 +465,10 @@ inline Model::Model(Model && model) : tools_dist_funs(std::move(model.tools_dist_funs)), // Entities entities(std::move(model.entities)), + prevalence_entity(std::move(model.prevalence_entity)), + prevalence_entity_as_proportion(std::move(model.prevalence_entity_as_proportion)), + entities_dist_funs(std::move(model.entities_dist_funs)), entities_backup(std::move(model.entities_backup)), - // prevalence_entity(std::move(model.prevalence_entity)), - // prevalence_entity_as_proportion(std::move(model.prevalence_entity_as_proportion)), - // entities_dist_funs(std::move(model.entities_dist_funs)), // Pseudo-RNG engine(std::move(model.engine)), runifd(std::move(model.runifd)), @@ -492,7 +493,8 @@ inline Model::Model(Model && model) : queue(std::move(model.queue)), use_queuing(model.use_queuing), array_double_tmp(model.array_double_tmp.size()), - array_virus_tmp(model.array_virus_tmp.size()) + array_virus_tmp(model.array_virus_tmp.size()), + array_int_tmp(model.array_int_tmp.size()) { db.model = this; @@ -542,10 +544,10 @@ inline Model & Model::operator=(const Model & m) tools_dist_funs = m.tools_dist_funs; entities = m.entities; + prevalence_entity = m.prevalence_entity; + prevalence_entity_as_proportion = m.prevalence_entity_as_proportion; + entities_dist_funs = m.entities_dist_funs; entities_backup = m.entities_backup; - // prevalence_entity = m.prevalence_entity; - // prevalence_entity_as_proportion = m.prevalence_entity_as_proportion; - // entities_dist_funs = m.entities_dist_funs; rewire_fun = m.rewire_fun; rewire_prop = m.rewire_prop; @@ -586,6 +588,7 @@ inline Model & Model::operator=(const Model & m) )); array_virus_tmp.resize(1024u); + array_int_tmp.resize(1024u * 1024); return *this; @@ -865,62 +868,62 @@ inline void Model::dist_tools() } -// template -// inline void Model::dist_entities() -// { +template +inline void Model::dist_entities() +{ + + // Starting first infection + int n = size(); + std::vector< size_t > idx(n); + for (epiworld_fast_uint e = 0; e < entities.size(); ++e) + { + + if (entities_dist_funs[e]) + { + + entities_dist_funs[e](entities[e], this); -// // Starting first infection -// int n = size(); -// std::vector< size_t > idx(n); -// for (epiworld_fast_uint e = 0; e < entities.size(); ++e) -// { - -// if (entities_dist_funs[e]) -// { - -// entities_dist_funs[e](entities[e], this); - -// } else { - -// // Picking how many -// int nsampled; -// if (prevalence_entity_as_proportion[e]) -// { -// nsampled = static_cast(std::floor(prevalence_entity[e] * size())); -// } -// else -// { -// nsampled = static_cast(prevalence_entity[e]); -// } - -// if (nsampled > static_cast(size())) -// throw std::range_error("There are only " + std::to_string(size()) + -// " individuals in the population. Cannot add the entity to " + std::to_string(nsampled)); + } else { + + // Picking how many + int nsampled; + if (prevalence_entity_as_proportion[e]) + { + nsampled = static_cast(std::floor(prevalence_entity[e] * size())); + } + else + { + nsampled = static_cast(prevalence_entity[e]); + } + + if (nsampled > static_cast(size())) + throw std::range_error("There are only " + std::to_string(size()) + + " individuals in the population. Cannot add the entity to " + std::to_string(nsampled)); -// Entity & entity = entities[e]; + Entity & entity = entities[e]; -// int n_left = n; -// std::iota(idx.begin(), idx.end(), 0); -// while (nsampled > 0) -// { -// int loc = static_cast(floor(runif() * n_left--)); + int n_left = n; + std::iota(idx.begin(), idx.end(), 0); + while (nsampled > 0) + { + int loc = static_cast(floor(runif() * n_left--)); -// population[idx[loc]].add_entity(entity, this, entity.state_init, entity.queue_init); + population[idx[loc]].add_entity(entity, this, entity.state_init, entity.queue_init); -// nsampled--; + nsampled--; -// std::swap(idx[loc], idx[n_left]); + std::swap(idx[loc], idx[n_left]); -// } + } -// } + } -// // Apply the events -// events_run(); + // Apply the events + events_run(); -// } + } -// } +} template inline void Model::chrono_start() { @@ -1201,6 +1204,32 @@ inline void Model::add_entity(Entity e) } +template +inline void Model::add_entity_n(Entity e, epiworld_fast_uint preval) +{ + + e.model = this; + e.id = entities.size(); + entities.push_back(e); + prevalence_entity.push_back(preval); + prevalence_entity_as_proportion.push_back(false); + entities_dist_funs.push_back(nullptr); + +} + +template +inline void Model::add_entity_fun(Entity e, EntityToAgentFun fun) +{ + + e.model = this; + e.id = entities.size(); + entities.push_back(e); + prevalence_entity.push_back(0.0); + prevalence_entity_as_proportion.push_back(false); + entities_dist_funs.push_back(fun); + +} + template inline void Model::rm_virus(size_t virus_pos) { @@ -1291,7 +1320,6 @@ inline void Model::load_agents_entities_ties( throw std::logic_error("The file " + fn + " was not found."); int linenum = 0; - std::vector< epiworld_fast_uint > source_; std::vector< std::vector< epiworld_fast_uint > > target_(entities.size()); target_.reserve(1e5); @@ -1332,36 +1360,64 @@ inline void Model::load_agents_entities_ties( } - // // Iterating over entities - // for (size_t e = 0u; e < entities.size(); ++e) - // { + return; - // // This entity will have individuals assigned to it, so we add it - // if (target_[e].size() > 0u) - // { +} - // // Filling in the gaps - // prevalence_entity[e] = static_cast(target_[e].size()); - // prevalence_entity_as_proportion[e] = false; +template +inline void Model::load_agents_entities_ties( + const std::vector< int > & agents_ids, + const std::vector< int > & entities_ids +) { - // // Generating the assignment function - // auto who = target_[e]; - // entities_dist_funs[e] = - // [who](Entity & e, Model* m) -> void { + if (agents_ids.size() != entities_ids.size()) + throw std::length_error( + std::string("agents_ids (") + + std::to_string(agents_ids.size()) + + std::string(") and entities_ids (") + + std::to_string(entities_ids.size()) + + std::string(") should match.") + ); - // for (auto w : who) - // m->population[w].add_entity(e, m, e.state_init, e.queue_init); - - // return; - - // }; - // } + size_t n_entries = agents_ids.size(); + for (size_t i = 0u; i < n_entries; ++i) + { + + if (agents_ids[i] >= this->population.size()) + throw std::length_error( + std::string("agents_ids[") + + std::to_string(i) + + std::string("] = ") + + std::to_string(agents_ids[i]) + + std::string(" is out of range (population size: ") + + std::to_string(this->population.size()) + + std::string(").") + ); - // } + + if (entities_ids[i] >= this->entities.size()) + throw std::length_error( + std::string("entities_ids[") + + std::to_string(i) + + std::string("] = ") + + std::to_string(entities_ids[i]) + + std::string(" is out of range (entities size: ") + + std::to_string(this->entities.size()) + + std::string(").") + ); + + // Adding the entity to the agent + this->population[agents_ids[i]].add_entity( + this->entities[entities_ids[i]], + nullptr /* Immediately add it to the agent */ + ); + + } return; + } template @@ -1488,6 +1544,7 @@ inline Model & Model::run( array_virus_tmp.resize(1024); + array_int_tmp.resize(1024 * 1024); // Checking whether the proposed state in/out/removed // are valid @@ -2072,11 +2129,9 @@ inline void Model::reset() { #endif } - else - { - for (auto & e: entities) - e.reset(); - } + + for (auto & e: entities) + e.reset(); current_date = 0; @@ -2089,6 +2144,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..a61ae91d 100644 --- a/include/epiworld/models/models.hpp +++ b/include/epiworld/models/models.hpp @@ -19,6 +19,8 @@ namespace epimodels { #include "seirdconnected.hpp" #include "sirlogit.hpp" #include "diffnet.hpp" + #include "seirmixing.hpp" + #include "sirmixing.hpp" } diff --git a/include/epiworld/models/seirmixing.hpp b/include/epiworld/models/seirmixing.hpp new file mode 100644 index 00000000..91c86a93 --- /dev/null +++ b/include/epiworld/models/seirmixing.hpp @@ -0,0 +1,520 @@ +#ifndef EPIWORLD_MODELS_SEIRMIXING_HPP +#define EPIWORLD_MODELS_SEIRMIXING_HPP + +/** + * @file seirentitiesconnected.hpp + * @brief Template for a Susceptible-Exposed-Infected-Removed (SEIR) model with mixing + */ +template +class ModelSEIRMixing : public epiworld::Model +{ +private: + std::vector< std::vector< epiworld::Agent * > > infected; + void update_infected(); + std::vector< epiworld::Agent * > sampled_agents; + size_t sample_agents( + epiworld::Agent * agent, + std::vector< epiworld::Agent * > & sampled_agents + ); + double adjusted_contact_rate; + std::vector< double > contact_matrix; + + size_t index(size_t i, size_t j, size_t n) { + return j * n + i; + } + +public: + + static const int SUSCEPTIBLE = 0; + static const int EXPOSED = 1; + static const int INFECTED = 2; + static const int RECOVERED = 3; + + ModelSEIRMixing() {}; + + /** + * @brief Constructs a ModelSEIRMixing object. + * + * @param model A reference to an existing ModelSEIRMixing object. + * @param vname The name of the ModelSEIRMixing 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 contact_matrix The contact matrix between entities in the model. + */ + ModelSEIRMixing( + ModelSEIRMixing & 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< double > contact_matrix + ); + + /** + * @brief Constructs a ModelSEIRMixing object. + * + * @param vname The name of the ModelSEIRMixing 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 contact_matrix The contact matrix between entities in the model. + */ + ModelSEIRMixing( + 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< double > contact_matrix + ); + + ModelSEIRMixing & 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. + */ + ModelSEIRMixing & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); + + size_t get_n_infected(size_t group) const + { + return infected[group].size(); + } + + void set_contact_matrix(std::vector< double > cmat) + { + contact_matrix = cmat; + return; + }; + +}; + +template +inline void ModelSEIRMixing::update_infected() +{ + + auto & agents = Model::get_agents(); + auto & entities = Model::get_entities(); + + infected.resize(entities.size()); + sampled_agents.resize(agents.size()); + + // Checking contact matrix's rows add to one + size_t nentities = entities.size(); + if (this->contact_matrix.size() != nentities*nentities) + throw std::length_error( + std::string("The contact matrix must be a square matrix of size ") + + std::string("nentities x nentities. ") + + std::to_string(this->contact_matrix.size()) + + std::string(" != ") + std::to_string(nentities*nentities) + + std::string(".") + ); + + for (size_t i = 0u; i < entities.size(); ++i) + { + double sum = 0.0; + for (size_t j = 0u; j < entities.size(); ++j) + { + if (this->contact_matrix[index(i, j, nentities)] < 0.0) + throw std::range_error( + std::string("The contact matrix must be non-negative. ") + + std::to_string(this->contact_matrix[index(i, j, nentities)]) + + std::string(" < 0.") + ); + sum += this->contact_matrix[index(i, j, nentities)]; + } + if (sum < 0.999 || sum > 1.001) + throw std::range_error( + std::string("The contact matrix must have rows that add to one. ") + + std::to_string(sum) + + std::string(" != 1.") + ); + } + + for (size_t i = 0; i < entities.size(); ++i) + { + infected[i].clear(); + infected[i].reserve(agents.size()); + } + + for (auto & a : agents) + { + + if (a.get_state() == ModelSEIRMixing::INFECTED) + infected[a.get_entity(0u).get_id()].push_back(&a); + + } + + // Adjusting contact rate + adjusted_contact_rate = Model::get_param("Contact rate") / + agents.size(); + + return; + +} + +template +inline size_t ModelSEIRMixing::sample_agents( + epiworld::Agent * agent, + std::vector< epiworld::Agent * > & sampled_agents + ) +{ + + size_t agent_group_id = agent->get_entity(0u).get_id(); + size_t ngroups = infected.size(); + + int samp_id = 0; + for (size_t g = 0; g < infected.size(); ++g) + { + + // How many from this entity? + int nsamples = epiworld::Model::rbinom( + infected[g].size(), + adjusted_contact_rate * contact_matrix[ + index(agent_group_id, g, ngroups) + ] + ); + + if (nsamples == 0) + continue; + + // Sampling from the entity + for (int s = 0; s < nsamples; ++s) + { + + // Randomly selecting an agent + int which = epiworld::Model::runif() * infected[g].size(); + + // Correcting overflow error + if (which >= static_cast(infected[g].size())) + which = static_cast(infected[g].size()) - 1; + + auto & a = infected[g][which]; + + // Can't sample itself + if (a->get_id() == agent->get_id()) + continue; + + sampled_agents[samp_id++] = a; + + } + + } + + return samp_id; + +} + +template +inline ModelSEIRMixing & ModelSEIRMixing::run( + epiworld_fast_uint ndays, + int seed +) +{ + + Model::run(ndays, seed); + return *this; + +} + +template +inline void ModelSEIRMixing::reset() +{ + + Model::reset(); + this->update_infected(); + + return; + +} + +template +inline Model * ModelSEIRMixing::clone_ptr() +{ + + ModelSEIRMixing * ptr = new ModelSEIRMixing( + *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 ModelSEIRMixing::ModelSEIRMixing( + ModelSEIRMixing & 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< double > contact_matrix + ) +{ + + // Setting up the contact matrix + this->contact_matrix = contact_matrix; + + epiworld::UpdateFun update_susceptible = []( + epiworld::Agent * p, epiworld::Model * m + ) -> void + { + + if (p->get_n_entities() == 0) + return; + + // Downcasting to retrieve the sampler attached to the + // class + ModelSEIRMixing * m_down = + dynamic_cast *>(m); + + size_t ndraws = m_down->sample_agents(p, m_down->sampled_agents); + + if (ndraws == 0u) + return; + + + // Drawing from the set + int nviruses_tmp = 0; + for (size_t n = 0u; n < ndraws; ++n) + { + + auto & neighbor = m_down->sampled_agents[n]; + + 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); + + } + + // Running the roulette + int which = roulette(nviruses_tmp, m); + + if (which < 0) + return; + + p->set_virus( + *m->array_virus_tmp[which], + m, + ModelSEIRMixing::EXPOSED + ); + + return; + + }; + + epiworld::UpdateFun update_infected = []( + epiworld::Agent * p, epiworld::Model * m + ) -> void { + + auto state = p->get_state(); + + if (state == ModelSEIRMixing::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, ModelSEIRMixing::INFECTED); + return; + + } + + + } else if (state == ModelSEIRMixing::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"); + + // Global function + epiworld::GlobalFun update = [](epiworld::Model * m) -> void + { + + ModelSEIRMixing * m_down = + dynamic_cast *>(m); + + m_down->update_infected(); + + return; + + }; + + model.add_globalevent(update, "Update infected individuals"); + + + // Preparing the virus ------------------------------------------- + epiworld::Virus virus(vname); + virus.set_state( + ModelSEIRMixing::EXPOSED, + ModelSEIRMixing::RECOVERED, + ModelSEIRMixing::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) with Mixing"); + + return; + +} + +template +inline ModelSEIRMixing::ModelSEIRMixing( + 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< double > contact_matrix + ) +{ + + this->contact_matrix = contact_matrix; + + ModelSEIRMixing( + *this, + vname, + n, + prevalence, + contact_rate, + transmission_rate, + avg_incubation_days, + recovery_rate, + contact_matrix + ); + + return; + +} + +template +inline ModelSEIRMixing & ModelSEIRMixing::initial_states( + std::vector< double > proportions_, + std::vector< int > /* queue_ */ +) +{ + + Model::initial_states_fun = + create_init_function_seir(proportions_) + ; + + return *this; + +} + +#endif diff --git a/include/epiworld/models/sirmixing.hpp b/include/epiworld/models/sirmixing.hpp new file mode 100644 index 00000000..74c2169a --- /dev/null +++ b/include/epiworld/models/sirmixing.hpp @@ -0,0 +1,493 @@ +#ifndef EPIWORLD_MODELS_SIRMIXING_HPP +#define EPIWORLD_MODELS_SIRMIXING_HPP + +/** + * @file seirentitiesconnected.hpp + * @brief Template for a Susceptible-Exposed-Infected-Removed (SEIR) model with mixing + */ +template +class ModelSIRMixing : public epiworld::Model +{ +private: + std::vector< std::vector< epiworld::Agent * > > infected; + void update_infected_list(); + std::vector< epiworld::Agent * > sampled_agents; + size_t sample_agents( + epiworld::Agent * agent, + std::vector< epiworld::Agent * > & sampled_agents + ); + double adjusted_contact_rate; + std::vector< double > contact_matrix; + + size_t index(size_t i, size_t j, size_t n) { + return j * n + i; + } + +public: + + static const int SUSCEPTIBLE = 0; + static const int INFECTED = 1; + static const int RECOVERED = 2; + + ModelSIRMixing() {}; + + /** + * @brief Constructs a ModelSIRMixing object. + * + * @param model A reference to an existing ModelSIRMixing object. + * @param vname The name of the ModelSIRMixing 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 recovery_rate The recovery rate of the disease in the model. + * @param contact_matrix The contact matrix between entities in the model. + */ + ModelSIRMixing( + ModelSIRMixing & model, + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double recovery_rate, + std::vector< double > contact_matrix + ); + + /** + * @brief Constructs a ModelSIRMixing object. + * + * @param vname The name of the ModelSIRMixing 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 recovery_rate The recovery rate of the disease in the model. + * @param contact_matrix The contact matrix between entities in the model. + */ + ModelSIRMixing( + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double recovery_rate, + std::vector< double > contact_matrix + ); + + ModelSIRMixing & 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. + */ + ModelSIRMixing & initial_states( + std::vector< double > proportions_, + std::vector< int > queue_ = {} + ); + + size_t get_n_infected(size_t group) const + { + return infected[group].size(); + } + + void set_contact_matrix(std::vector< double > cmat) + { + contact_matrix = cmat; + return; + }; + +}; + +template +inline void ModelSIRMixing::update_infected_list() +{ + + auto & agents = Model::get_agents(); + auto & entities = Model::get_entities(); + + infected.resize(entities.size()); + sampled_agents.resize(agents.size()); + + // Checking contact matrix's rows add to one + size_t nentities = entities.size(); + if (this->contact_matrix.size() != nentities*nentities) + throw std::length_error( + std::string("The contact matrix must be a square matrix of size ") + + std::string("nentities x nentities. ") + + std::to_string(this->contact_matrix.size()) + + std::string(" != ") + std::to_string(nentities*nentities) + + std::string(".") + ); + + for (size_t i = 0u; i < entities.size(); ++i) + { + double sum = 0.0; + for (size_t j = 0u; j < entities.size(); ++j) + { + if (this->contact_matrix[index(i, j, nentities)] < 0.0) + throw std::range_error( + std::string("The contact matrix must be non-negative. ") + + std::to_string(this->contact_matrix[index(i, j, nentities)]) + + std::string(" < 0.") + ); + sum += this->contact_matrix[index(i, j, nentities)]; + } + if (sum < 0.999 || sum > 1.001) + throw std::range_error( + std::string("The contact matrix must have rows that add to one. ") + + std::to_string(sum) + + std::string(" != 1.") + ); + } + + for (size_t i = 0; i < entities.size(); ++i) + { + infected[i].clear(); + infected[i].reserve(agents.size()); + } + + for (auto & a : agents) + { + + if (a.get_state() == ModelSIRMixing::INFECTED) + infected[a.get_entity(0u).get_id()].push_back(&a); + + } + + // Adjusting contact rate + adjusted_contact_rate = Model::get_param("Contact rate") / + agents.size(); + + return; + +} + +template +inline size_t ModelSIRMixing::sample_agents( + epiworld::Agent * agent, + std::vector< epiworld::Agent * > & sampled_agents + ) +{ + + size_t agent_group_id = agent->get_entity(0u).get_id(); + size_t ngroups = infected.size(); + + int samp_id = 0; + for (size_t g = 0; g < infected.size(); ++g) + { + + // How many from this entity? + int nsamples = epiworld::Model::rbinom( + infected[g].size(), + adjusted_contact_rate * contact_matrix[ + index(agent_group_id, g, ngroups) + ] + ); + + if (nsamples == 0) + continue; + + // Sampling from the entity + for (int s = 0; s < nsamples; ++s) + { + + // Randomly selecting an agent + int which = epiworld::Model::runif() * infected[g].size(); + + // Correcting overflow error + if (which >= static_cast(infected[g].size())) + which = static_cast(infected[g].size()) - 1; + + auto & a = infected[g][which]; + + // Can't sample itself + if (a->get_id() == agent->get_id()) + continue; + + sampled_agents[samp_id++] = a; + + } + + } + + return samp_id; + +} + +template +inline ModelSIRMixing & ModelSIRMixing::run( + epiworld_fast_uint ndays, + int seed +) +{ + + Model::run(ndays, seed); + return *this; + +} + +template +inline void ModelSIRMixing::reset() +{ + + Model::reset(); + this->update_infected_list(); + + return; + +} + +template +inline Model * ModelSIRMixing::clone_ptr() +{ + + ModelSIRMixing * ptr = new ModelSIRMixing( + *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 ModelSIRMixing::ModelSIRMixing( + ModelSIRMixing & model, + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double recovery_rate, + std::vector< double > contact_matrix + ) +{ + + // Setting up the contact matrix + this->contact_matrix = contact_matrix; + + epiworld::UpdateFun update_susceptible = []( + epiworld::Agent * p, epiworld::Model * m + ) -> void + { + + if (p->get_n_entities() == 0) + return; + + // Downcasting to retrieve the sampler attached to the + // class + ModelSIRMixing * m_down = + dynamic_cast *>(m); + + size_t ndraws = m_down->sample_agents(p, m_down->sampled_agents); + + if (ndraws == 0u) + return; + + + // Drawing from the set + int nviruses_tmp = 0; + for (size_t n = 0u; n < ndraws; ++n) + { + + auto & neighbor = m_down->sampled_agents[n]; + + 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); + + } + + // Running the roulette + int which = roulette(nviruses_tmp, m); + + if (which < 0) + return; + + p->set_virus( + *m->array_virus_tmp[which], + m, + ModelSIRMixing::INFECTED + ); + + return; + + }; + + epiworld::UpdateFun update_infected = []( + epiworld::Agent * p, epiworld::Model * m + ) -> void { + + auto state = p->get_state(); + + if (state == ModelSIRMixing::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 infected."); + } + #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 infected individuals. (SIR)") ; + + 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"); + + // state + model.add_state("Susceptible", update_susceptible); + model.add_state("Infected", update_infected); + model.add_state("Recovered"); + + // Global function + epiworld::GlobalFun update = [](epiworld::Model * m) -> void + { + + ModelSIRMixing * m_down = + dynamic_cast *>(m); + + m_down->update_infected_list(); + + return; + + }; + + model.add_globalevent(update, "Update infected individuals"); + + + // Preparing the virus ------------------------------------------- + epiworld::Virus virus(vname); + virus.set_state( + ModelSIRMixing::INFECTED, + ModelSIRMixing::RECOVERED, + ModelSIRMixing::RECOVERED + ); + + virus.set_prob_infecting(&model("Prob. Transmission")); + virus.set_prob_recovery(&model("Prob. Recovery")); + + model.add_virus(virus, prevalence); + + model.queuing_off(); // No queuing need + + // Adding the empty population + model.agents_empty_graph(n); + + model.set_name("Susceptible-Infected-Removed (SIR) with Mixing"); + + return; + +} + +template +inline ModelSIRMixing::ModelSIRMixing( + std::string vname, + epiworld_fast_uint n, + epiworld_double prevalence, + epiworld_double contact_rate, + epiworld_double transmission_rate, + epiworld_double recovery_rate, + std::vector< double > contact_matrix + ) +{ + + this->contact_matrix = contact_matrix; + + ModelSIRMixing( + *this, + vname, + n, + prevalence, + contact_rate, + transmission_rate, + recovery_rate, + contact_matrix + ); + + return; + +} + +template +inline ModelSIRMixing & ModelSIRMixing::initial_states( + std::vector< double > proportions_, + std::vector< int > /* queue_ */ +) +{ + + Model::initial_states_fun = + create_init_function_sir(proportions_) + ; + + return *this; + +} + +#endif diff --git a/paper/mixing.md b/paper/mixing.md new file mode 100644 index 00000000..2ac7d373 --- /dev/null +++ b/paper/mixing.md @@ -0,0 +1,159 @@ +# Mixing probabilities in connected model +George G. Vega Yon, Ph.D. +2024-04-25 + +## Case 1: No grouping + +We will look into the probability of drawing infected individuals to +simplify the algorithm. There are $I$ infected individuals at any time +in the simulation; thus, instead of drawing from $Bern(c/N, N)$, we will +be drawing from $Bern(c/N, I)$. The next step is to check which infected +individuals should be drawn. Let’s compare the distributions using the +hypergeometric as an example: + +``` r +set.seed(132) +nsims <- 1e4 +N <- 400 +rate <- 5 +p <- rate/N +I <- 10 + +sim_complex <- parallel::mclapply(1:nsims, \(i) { + nsamples <- rbinom(N, N, p) + sum(rbinom(N, size = nsamples, prob = I/N) > 0) +}, mc.cores = 4L) |> unlist() + +sim_simple <- parallel::mclapply(1:nsims, \(i) { + sum(rbinom(N, I, p) > 0) +}, mc.cores = 4L) |> unlist() + + +op <- par(mfrow = c(1,2)) +MASS::truehist(sim_complex) +MASS::truehist(sim_simple) +``` + +![](mixing_files/figure-commonmark/Simulation-1.png) + +``` r +par(op) + +quantile(sim_complex) +``` + + 0% 25% 50% 75% 100% + 27 43 47 51 71 + +``` r +quantile(sim_simple) +``` + + 0% 25% 50% 75% 100% + 23 43 47 51 71 + +``` r +plotter(sim_complex, sim_simple) +``` + +![](mixing_files/figure-commonmark/Simulation-2.png) + +These two approaches are equivalent, but the second one is more +efficient from the computational perspective. + +## Case 2: Grouping + +This explores the case when we have mixing across groups. The question +is if we can replicate the effect at the group level. + +``` r +set.seed(123133) + +ngroups <- 3 +mixing <- matrix( + c(0.1, 0.2, 0.3, 0.2, 0.1, 0.2, 0.3, 0.2, 0.1), + nrow = ngroups, + ncol = ngroups + ) + +mixing <- mixing/rowSums(mixing) +mixing +``` + + [,1] [,2] [,3] + [1,] 0.1666667 0.3333333 0.5000000 + [2,] 0.4000000 0.2000000 0.4000000 + [3,] 0.5000000 0.3333333 0.1666667 + +``` r +N <- 500 +sizes <- c(100, 150, 250) +rate <- 5 +p <- rate/N +I <- c(10, 30, 20) + +ids <- rep.int(1:ngroups, times = sizes) + +nsims <- 1e4 + +sim_complex <- parallel::mclapply(1:nsims, \(i) { + + # Sampling group first + sapply(1:ngroups, \(g) { + + # How many each individual will sample from the groups + ans <- rbinom( + n = N, size = sizes[g], prob = mixing[ids,][,g] * p + ) |> sum() + + # Sampling with replacement + rbinom(ans, size = 1, prob = I[g]/sizes[g]) |> sum() + + }) |> sum() + +}, mc.cores = 4L) |> unlist() +``` + +Using the alternative method in which we directly weight the +probabilities: + +``` r +sim_simple <- parallel::mclapply(1:nsims, \(i) { + + # Sampling group first + sapply(1:ngroups, \(g) { + rbinom( + n = N, size = I[g], prob = mixing[cbind(ids,g)] * p + ) |> sum() + }) |> sum() + +}, mc.cores = 4L) |> unlist() + +op <- par(mfrow = c(1,2)) +MASS::truehist(sim_complex) +MASS::truehist(sim_simple) +``` + +![](mixing_files/figure-commonmark/sim-with-mixing-2-1.png) + +``` r +par(op) + +quantile(sim_complex) +``` + + 0% 25% 50% 75% 100% + 57 88 94 101 131 + +``` r +quantile(sim_simple) +``` + + 0% 25% 50% 75% 100% + 58 87 94 101 135 + +``` r +plotter(sim_complex, sim_simple) +``` + +![](mixing_files/figure-commonmark/sim-with-mixing-2-2.png) diff --git a/paper/mixing.qmd b/paper/mixing.qmd new file mode 100644 index 00000000..adf719eb --- /dev/null +++ b/paper/mixing.qmd @@ -0,0 +1,138 @@ +--- +format: gfm +title: Mixing probabilities in connected model +author: George G. Vega Yon, Ph.D. +date: 2024-04-25 +--- + +```{r} +#| label: setup +#| echo: false +tabulator <- function(x) { + x <- table(x) |> prop.table() |> cumsum() |> data.frame() + + data.frame( + x = as.integer(rownames(x)), + y = x[[1]] + ) + +} + +plotter <- function(a, b) { + a <- tabulator(a) + b <- tabulator(b) + + plot(a, col = "red") + points(b, col = "blue") +} +``` + +## Case 1: No grouping + +We will look into the probability of drawing infected individuals to simplify the algorithm. There are $I$ infected individuals at any time in the simulation; thus, instead of drawing from $Bern(c/N, N)$, we will be drawing from $Bern(c/N, I)$. The next step is to check which infected individuals should be drawn. Let's compare the distributions: + + +```{r} +#| label: Simulation +set.seed(132) +nsims <- 1e4 +N <- 400 +rate <- 5 +p <- rate/N +I <- 10 + +sim_complex <- parallel::mclapply(1:nsims, \(i) { + nsamples <- rbinom(N, N, p) + sum(rbinom(N, size = nsamples, prob = I/N) > 0) +}, mc.cores = 4L) |> unlist() + +sim_simple <- parallel::mclapply(1:nsims, \(i) { + sum(rbinom(N, I, p) > 0) +}, mc.cores = 4L) |> unlist() + + +op <- par(mfrow = c(1,2)) +MASS::truehist(sim_complex) +MASS::truehist(sim_simple) +par(op) + +quantile(sim_complex) +quantile(sim_simple) + +plotter(sim_complex, sim_simple) +``` + +These two approaches are equivalent, but the second one is more efficient from the computational perspective. + +## Case 2: Grouping + +This explores the case when we have mixing across groups. The question is if we can replicate the effect at the group level. + +```{r} +#| label: sim-with-mixing-1 +set.seed(123133) + +ngroups <- 3 +mixing <- matrix( + c(0.1, 0.2, 0.3, 0.2, 0.1, 0.2, 0.3, 0.2, 0.1), + nrow = ngroups, + ncol = ngroups + ) + +mixing <- mixing/rowSums(mixing) +mixing + +N <- 500 +sizes <- c(100, 150, 250) +rate <- 5 +p <- rate/N +I <- c(10, 30, 20) + +ids <- rep.int(1:ngroups, times = sizes) + +nsims <- 1e4 + +sim_complex <- parallel::mclapply(1:nsims, \(i) { + + # Sampling group first + sapply(1:ngroups, \(g) { + + # How many each individual will sample from the groups + ans <- rbinom( + n = N, size = sizes[g], prob = mixing[ids,][,g] * p + ) |> sum() + + # Sampling with replacement + rbinom(ans, size = 1, prob = I[g]/sizes[g]) |> sum() + + }) |> sum() + +}, mc.cores = 4L) |> unlist() +``` + +Using the alternative method in which we directly weight the probabilities: + +```{r} +#| label: sim-with-mixing-2 +sim_simple <- parallel::mclapply(1:nsims, \(i) { + + # Sampling group first + sapply(1:ngroups, \(g) { + rbinom( + n = N, size = I[g], prob = mixing[cbind(ids,g)] * p + ) |> sum() + }) |> sum() + +}, mc.cores = 4L) |> unlist() + +op <- par(mfrow = c(1,2)) +MASS::truehist(sim_complex) +MASS::truehist(sim_simple) +par(op) + +quantile(sim_complex) +quantile(sim_simple) + +plotter(sim_complex, sim_simple) + +``` \ No newline at end of file diff --git a/paper/mixing_files/figure-commonmark/Simulation-1.png b/paper/mixing_files/figure-commonmark/Simulation-1.png new file mode 100644 index 00000000..0a55ade6 Binary files /dev/null and b/paper/mixing_files/figure-commonmark/Simulation-1.png differ diff --git a/paper/mixing_files/figure-commonmark/Simulation-2.png b/paper/mixing_files/figure-commonmark/Simulation-2.png new file mode 100644 index 00000000..c40171f1 Binary files /dev/null and b/paper/mixing_files/figure-commonmark/Simulation-2.png differ diff --git a/paper/mixing_files/figure-commonmark/sim-with-mixing-2-1.png b/paper/mixing_files/figure-commonmark/sim-with-mixing-2-1.png new file mode 100644 index 00000000..abae2785 Binary files /dev/null and b/paper/mixing_files/figure-commonmark/sim-with-mixing-2-1.png differ diff --git a/paper/mixing_files/figure-commonmark/sim-with-mixing-2-2.png b/paper/mixing_files/figure-commonmark/sim-with-mixing-2-2.png new file mode 100644 index 00000000..ad53d2ce Binary files /dev/null and b/paper/mixing_files/figure-commonmark/sim-with-mixing-2-2.png differ diff --git a/tests/05-mixing.cpp b/tests/05-mixing.cpp new file mode 100644 index 00000000..67396a25 --- /dev/null +++ b/tests/05-mixing.cpp @@ -0,0 +1,137 @@ +#ifndef CATCH_CONFIG_MAIN +#define EPI_DEBUG +#endif + +#include "tests.hpp" + +using namespace epiworld; + +EPIWORLD_TEST_CASE("SEIRMixing", "[SEIR-mixing]") { + + std::vector< double > contact_matrix = { + 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0 + }; + + epimodels::ModelSEIRMixing<> model( + "Flu", // std::string vname, + 10000, // epiworld_fast_uint n, + 0.01,// epiworld_double prevalence, + 40.0,// epiworld_double contact_rate, + 1.0,// epiworld_double transmission_rate, + 1.0,// epiworld_double avg_incubation_days, + 1.0/2.0,// epiworld_double recovery_rate, + contact_matrix + ); + + // Copy the original virus + Virus<> v1 = model.get_virus(0); + model.rm_virus(0); + + model.add_virus_fun(v1, dist_virus<>(0)); + + // Creating three groups + Entity<> e1("Entity 1"); + Entity<> e2("Entity 2"); + Entity<> e3("Entity 3"); + + model.add_entity_fun(e1, dist_factory<>(0, 3000)); + model.add_entity_fun(e2, dist_factory<>(3000, 6000)); + model.add_entity_fun(e3, dist_factory<>(6000, 10000)); + + // Running and checking the results + model.run(50, 123); + model.print(); + + // Getting all agents + int n_right = 0; + int n_wrong = 0; + + for (const auto & a : model.get_agents()) + { + if (a.get_state() != epimodels::ModelSEIRMixing::SUSCEPTIBLE) + { + if (a.get_entity(0).get_id() == 0) + { + n_right++; + continue; + } + + n_wrong++; + + } + + } + + #ifdef CATCH_CONFIG_MAIN + REQUIRE_FALSE((n_wrong != 0 | n_right != 3000)); + #endif + + // Reruning the model where individuals from group 0 transmit all to group 1 + contact_matrix[0] = 0.0; + contact_matrix[6] = 1.0; + contact_matrix[4] = 0.5; + contact_matrix[1] = 0.5; + model.set_contact_matrix(contact_matrix); + + // Running and checking the results + model.run(50, 123); + model.print(); + + // Getting all agents + n_right = 0; + n_wrong = 0; + + for (const auto & a : model.get_agents()) + { + + if (a.get_id() == 0) + { + n_right++; + } + else if (a.get_state() != epimodels::ModelSEIRMixing::SUSCEPTIBLE) + { + if (a.get_entity(0).get_id() == 1) + { + n_right++; + continue; + } + + n_wrong++; + + } + + } + + #ifdef CATCH_CONFIG_MAIN + REQUIRE_FALSE((n_wrong != 0 | n_right != 3001)); + #endif + + // Rerunning with plain mixing + std::fill(contact_matrix.begin(), contact_matrix.end(), 1.0/3.0); + model.set_contact_matrix(contact_matrix); + + // Running and checking the results + model.run(50, 123); + model.print(); + + std::vector< int > totals; + model.get_db().get_today_total(nullptr, &totals); + + std::vector< int > expected_totals = { + 0, 0, 0, + static_cast(model.size()) + }; + + #ifdef CATCH_CONFIG_MAIN + REQUIRE_THAT(totals, Catch::Equals(expected_totals)); + #endif + + + + #ifndef CATCH_CONFIG_MAIN + return 0; + #endif + +} diff --git a/tests/06-mixing.cpp b/tests/06-mixing.cpp new file mode 100644 index 00000000..18abce10 --- /dev/null +++ b/tests/06-mixing.cpp @@ -0,0 +1,138 @@ +#ifndef CATCH_CONFIG_MAIN +#define EPI_DEBUG +#endif + +#include "tests.hpp" + +using namespace epiworld; + + + +EPIWORLD_TEST_CASE("SIRMixing", "[SIR-mixing]") { + + std::vector< double > contact_matrix = { + 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0 + }; + + epimodels::ModelSIRMixing<> model( + "Flu", // std::string vname, + 10000, // epiworld_fast_uint n, + 0.01,// epiworld_double prevalence, + 40.0,// epiworld_double contact_rate, + 1.0,// epiworld_double transmission_rate, + 1.0/2.0,// epiworld_double recovery_rate, + contact_matrix + ); + + // Copy the original virus + Virus<> v1 = model.get_virus(0); + model.rm_virus(0); + + model.add_virus_fun(v1, dist_virus<>(0)); + + // Creating three groups + Entity<> e1("Entity 1"); + Entity<> e2("Entity 2"); + Entity<> e3("Entity 3"); + + model.add_entity_fun(e1, dist_factory<>(0, 3000)); + model.add_entity_fun(e2, dist_factory<>(3000, 6000)); + model.add_entity_fun(e3, dist_factory<>(6000, 10000)); + + // Running and checking the results + model.run(50, 123); + model.print(); + + // Getting all agents + int n_right = 0; + int n_wrong = 0; + + for (const auto & a : model.get_agents()) + { + if (a.get_state() != epimodels::ModelSIRMixing::SUSCEPTIBLE) + { + if (a.get_entity(0).get_id() == 0) + { + n_right++; + continue; + } + + n_wrong++; + + } + + } + + #ifdef CATCH_CONFIG_MAIN + REQUIRE_FALSE((n_wrong != 0 | n_right != 3000)); + #endif + + // Reruning the model where individuals from group 0 transmit all to group 1 + contact_matrix[0] = 0.0; + contact_matrix[6] = 1.0; + contact_matrix[4] = 0.5; + contact_matrix[1] = 0.5; + model.set_contact_matrix(contact_matrix); + + // Running and checking the results + model.run(50, 123); + model.print(); + + // Getting all agents + n_right = 0; + n_wrong = 0; + + for (const auto & a : model.get_agents()) + { + + if (a.get_id() == 0) + { + n_right++; + } + else if (a.get_state() != epimodels::ModelSIRMixing::SUSCEPTIBLE) + { + if (a.get_entity(0).get_id() == 1) + { + n_right++; + continue; + } + + n_wrong++; + + } + + } + + #ifdef CATCH_CONFIG_MAIN + REQUIRE_FALSE((n_wrong != 0 | n_right != 3001)); + #endif + + // Rerunning with plain mixing + std::fill(contact_matrix.begin(), contact_matrix.end(), 1.0/3.0); + model.set_contact_matrix(contact_matrix); + + // Running and checking the results + model.run(50, 123); + model.print(); + + std::vector< int > totals; + model.get_db().get_today_total(nullptr, &totals); + + std::vector< int > expected_totals = { + 0, 0, + static_cast(model.size()) + }; + + #ifdef CATCH_CONFIG_MAIN + REQUIRE_THAT(totals, Catch::Equals(expected_totals)); + #endif + + + + #ifndef CATCH_CONFIG_MAIN + return 0; + #endif + +} diff --git a/tests/Makefile b/tests/Makefile index b878ffbb..1d4bb2ca 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -19,6 +19,9 @@ main.a: main.cpp clean 01c.o: 01c-sir.cpp g++ -std=c++14 -Wall -Wextra -O2 -g $(OPENMP) -pedantic 01c-sir.cpp -o 01c.o +05-mixing.o: 05-mixing.cpp + g++ -std=c++14 -Wall -Wextra -O2 -g $(OPENMP) -pedantic 05-mixing.cpp -o 05-mixing.o + # Check coverage using the main.o target coverage: main.o g++ -g -O1 -fprofile-arcs -ftest-coverage main.cpp -o main.o && \ diff --git a/tests/main.cpp b/tests/main.cpp index 510a6bfa..3757a556 100644 --- a/tests/main.cpp +++ b/tests/main.cpp @@ -14,4 +14,6 @@ #include "01-sirconnected.cpp" #include "02-reproducible-sir.cpp" #include "02-reproducible-sirconn.cpp" -#include "04-initial-dist.cpp" \ No newline at end of file +#include "04-initial-dist.cpp" +#include "05-mixing.cpp" +#include "06-mixing.cpp" diff --git a/tests/tests.hpp b/tests/tests.hpp index dde5d153..c7210143 100644 --- a/tests/tests.hpp +++ b/tests/tests.hpp @@ -43,6 +43,33 @@ std::string file_reader(std::string fname) return res; } +template +inline epiworld::EntityToAgentFun dist_factory(int from, int to) { + return [from, to](epiworld::Entity & e, epiworld::Model * m) -> void { + + auto & agents = m->get_agents(); + for (int i = from; i < to; ++i) + { + e.add_agent(&agents[i], m); + } + + return; + + }; +} + +template +inline epiworld::VirusToAgentFun dist_virus(int i) +{ + return [i](epiworld::Virus & v, epiworld::Model * m) -> void { + + m->get_agents()[i].set_virus(v, m); + return; + + }; + +} + #ifndef CATCH_CONFIG_MAIN