From cf9c8ebc7cfdc00d7dbac58eb09523be9eee868b Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Sat, 27 Apr 2024 16:48:06 -0600 Subject: [PATCH] Mixing model (#14) * Adding entities back * Draft of entities model * Adding function to sample from contact matrix more efficiently * Sampler * Adding mixing matrix! Still need to set dist_ent without replacement * Moving group sampler to its own * Adding function to distribute entities based on data * Renaming model * Adding idea about how to sample more efficiently * Working on algorithm to speed up [skip ci] * Working on mixing doc [skip ci] * Solved issue about sampling * Good progress, need testing [skip ci] * Example with mixing (new algo) running * Test transmitting from one group to another almost done! * Adding missing test file * Adding tests for the SEIR mixing model * Updating version of epiworld * Adding the SIRMixing model + tests --- .gitignore | 4 +- .vscode/settings.json | 8 +- epiworld.hpp | 1575 +++++++++++++++-- examples/11-entities/Makefile | 2 + examples/11-entities/main.cpp | 55 + include/epiworld/agent-bones.hpp | 4 +- include/epiworld/agent-events-meat.hpp | 17 +- include/epiworld/agent-meat.hpp | 36 +- include/epiworld/agentssample-bones.hpp | 35 +- include/epiworld/entity-meat.hpp | 14 +- include/epiworld/epiworld.hpp | 7 +- include/epiworld/groupsampler-bones.hpp | 59 + include/epiworld/groupsampler-meat.hpp | 84 + include/epiworld/model-bones.hpp | 16 +- include/epiworld/model-meat.hpp | 222 ++- include/epiworld/models/models.hpp | 2 + include/epiworld/models/seirmixing.hpp | 520 ++++++ include/epiworld/models/sirmixing.hpp | 493 ++++++ paper/mixing.md | 159 ++ paper/mixing.qmd | 138 ++ .../figure-commonmark/Simulation-1.png | Bin 0 -> 15312 bytes .../figure-commonmark/Simulation-2.png | Bin 0 -> 20910 bytes .../figure-commonmark/sim-with-mixing-2-1.png | Bin 0 -> 14857 bytes .../figure-commonmark/sim-with-mixing-2-2.png | Bin 0 -> 23912 bytes tests/05-mixing.cpp | 137 ++ tests/06-mixing.cpp | 138 ++ tests/Makefile | 3 + tests/main.cpp | 4 +- tests/tests.hpp | 27 + 29 files changed, 3523 insertions(+), 236 deletions(-) create mode 100644 examples/11-entities/Makefile create mode 100644 examples/11-entities/main.cpp create mode 100644 include/epiworld/groupsampler-bones.hpp create mode 100644 include/epiworld/groupsampler-meat.hpp create mode 100644 include/epiworld/models/seirmixing.hpp create mode 100644 include/epiworld/models/sirmixing.hpp create mode 100644 paper/mixing.md create mode 100644 paper/mixing.qmd create mode 100644 paper/mixing_files/figure-commonmark/Simulation-1.png create mode 100644 paper/mixing_files/figure-commonmark/Simulation-2.png create mode 100644 paper/mixing_files/figure-commonmark/sim-with-mixing-2-1.png create mode 100644 paper/mixing_files/figure-commonmark/sim-with-mixing-2-2.png create mode 100644 tests/05-mixing.cpp create mode 100644 tests/06-mixing.cpp diff --git a/.gitignore b/.gitignore index 37b160139..fbadc30d5 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 639f0a55d..4976af687 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 e18df010c..927837e48 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 000000000..0c4d9f21d --- /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 000000000..6db469bde --- /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 1e83c069f..cd703a32c 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 8da043fb9..275896e61 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 a5006b5d1..1d907ed80 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 c8462c585..e0d1285a0 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 4c5a30d42..25fdbad0f 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 5d7d7bdcd..5ed8afe71 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 000000000..88ec690ac --- /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 000000000..6e9609659 --- /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 3b51e2a83..7959c608e 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 63bdc3e23..37e93379d 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 72cfc8607..a61ae91d4 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 000000000..91c86a93d --- /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 000000000..74c2169aa --- /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 000000000..2ac7d3735 --- /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 000000000..adf719eb1 --- /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 0000000000000000000000000000000000000000..0a55ade6865e0adebcda025f6b70af7c2ef178b5 GIT binary patch literal 15312 zcmeIZXHb;gwl(@7B8mtiS(3KoEFd|ljhIN1vl1mIo1DZz6huHoqJok^a!!&ZDj<@Z zEEy!{JKMeQ+54Vdr|ysY{kZkj*Hv%5U3Rm2t@X?`=a^%R6{xHze~#!P5rQD+?%YN{ zK#&tj2!e+_eFFYvbNRM0{O^qYZEYt6A*sjxhu3VEW`ZD>kUQw>YOZlhBW`*crbowX zr1toJmy*(wv8QN0*_chfNlCY83lTOOT*(yvq^S6KE+pvvq&RJ$@n;KLg-K!_VcHl$ zowph@Z{KsJuv7$nWZy^G<8QBoJ#z^ODjSYl820a9|IRc!>@RBfdZBV$x?ilh(X8MM zg78yo_TeFj@1)8qg1kFVK#KUD`9E*}x7&x;au|fY*rfGMd;9(Ud0#wy_l@5l_(_rJ zA3uI*Wb2q3F7>AO^eD+>^OM?=5EzbC{rSc;;2$8CHyl5T%YEDcDTL%no$_H1lzH|Gm2d=4g= z^!4;2%}zX|5RSFk-&)ehRM$Os`tgpHhxpj*J0w^feiynUvC(CI(zv3}FJzZqguc*t*$G7q`Tel-yzx;`q&rxNN9)gH+ z!JKn_4i}=hG~pSL&5hL@7uyV3+Ir2#yNEk4rm@l@NN4MLyYGAkm8D9F(%YjRE34!6 zuyzTT<>B7mUd4Fvco`pQUS3{lY46FVR~dpQ5bS-}d_o*AAq@YMiO@-wL2J~6l{74B53k#c9K2O%#*(v?x_iL{rxAhs95lr+~NqYl5 zskwK&UJof0%BVL@r6*M>QOt>zfb`unXKm~4l9;MXA#yg49 zC{$pe zd_>ZA4*U6WbR(}}wImzc?7{B3MrL1szlf;luZ|7{j8Ls4gdjQV{ngdgfY8wDKiNB9 zzI?f6^)rO#>I94k&j^x|`q0mv9oHpILgZGhX`XDn0J<*twh8Vv1rv*S^QLf?^wi0d z&;9-V1$6278F+&Uo;-~7>&N|rzkjCUBga?~yJACUZq(qIl1P&o7Z+D>VBq(k!^6V~ z-g|Qkacp?@+T2N`#*_!@a%lA8Ktbxqk2f&qBU)NoZdV<9yDb&XZ+`hgcSNzf%c8E{ zvlwA5Vz;%v?g0_bKt=Ys$u!SbIxY$8e~yldjO>H>i!vcg=F68$HfSp=Gx$=;!iR_7 z+ZCSZ|G&T29LkU>egg4#=rzsL$h>y>@`qwWUdT2TE-Re{bIOUdnb@;LL~kM@iVD8s z;o&U}m6&VrJQN5m=XM^e_JZ9?vsr+JA0Mvx$Q+s5w*a8fG=9a zYyO1TMdGWk>91u{`4&v(Ye{+SPDjI;7rEfk+n5M3u>slkJ=eKpN{*_|C5LVE?3?~pvtZqI*{qLM$8DoqJ@7N4pR z8PsPDiVGMw5$~fL3k!wHwhfBhloK7avvh)~`Mg)EwhO02LWl!{9>fx2nzabR#Pe-O z${QReQ{>K)g_n#zU?AMQ7~Pn$dAL2sL?d*Y(Mt_I6=SZ}e4XJVcHxygQ`C?xIoUwr z9SlXUm8A8pD(9sXCv&p!m*>>wuQLcVU{;ZN(xv=ggT;>8+*tMl0#aE`T!wk(N0`I~ zQM3X1x|=S#@q~C17B8(17kfr;7!ZUNU5Af6yiq9^HaK9-SGnp}yaO|6>;7Oli6A*l z1YuU%USnPp6DU+ZCUc$}@uickn|$9ClFqwx%}Pxk4;e-}|BgT*`Co(mPj77p>Yn)= zttZkM{`qAbmQHdK!QK*Kje7f5*yqS|-%>+EgJC%`QH z+sH%3PoF54mls|VryF=}4U~+!xApgHTCg6Fs#1F1xpPNRQPD;>Tf<6D?nQzdQrpf= zl}CN>`}^HwEa!doRfzRmx3z84>ccHv?~0R1FZum^e5Cd{yDNg+5AY*|xA!q5jIxK} zF)_Pi$44c_+1lFLFhlnlypWfEgxE+TA|et-Nv$E%JY1cgo)(wS=no3Ic+E0}r@W#f z_0uN;0)ljflSGWZszjuuhGmYXHoB9OlS^Ls@e~vRxw!#ab~ROLX}2*Kl~<96&X633 zD%}c-J;%Qi>~`4o&ODWdx$u`6H-~=ejAhfW2nyJk|0GWGpCPwVjzJ4tk^R9pjh?O{$S@k>5{h0pX&-j2o-QXiO{#Uiv z{;8T(zhCm5{e~-tiv_&;acdGpOrtNEoJ)!!r6STyirq9=yLjwJ~=#DyqlPDN|GC zP>7j%s+BEK?fn^K{>n}2gN8uz{f)VsbrLbG3MNMf^V3r%u^pd<&LO?3PlUpggL)aXckN!~f(i@v>V`($!yiGiHF zz`QGgLBjdT;;*%WE@{>H)4rpuigE8^V$#J49x}e%IcwkqkZEPG$ZGmq6qL~gfMRc0 z?y0dk3KI|{dhc04k_@01e*yu7oqy+)90r*T0IOviB_*Y0z3dr$%F2M!^72sljqarD zlg(i(zw(mzf~(`&S|hGDhXp$a4pw{p$cdkqtySlry~QNuz9Aa<24!oDLbUJ zQ&xVK3#MB5nLSi!(bLT2Lx28!l6a_UjJSF^hMX)$Sm6S>dMaLypP0xVI|oO&*FLkP z9{M@HMU_c_)2cIsKOrgO!QMh{&C%Y%TOM6_+(PW^w=wKUXhR&GEYmHH9feq74B&>b zu^QMtF94B4B`aRrBeN?jVvJ6+zb#ROUo#qy2;}eqT6+RA! zO!7M^De01!m>3_Q5*lVEEiIjK6L02GfFCId$!L{_OMi}jZaHSQGd^3Vs63&|X`&&J zROaAIaW*s+92^{FtpNzVdCzBs^_h;L3KzQ$b2140&G|m7&Nwkr>7Asn*1i@rAF<70 zOiYxNQQ_fNAJ+PIcXtD_4ZQ8`?QLWWbVhpt|XkqvMs_WgAW-QK@pzr zTzpCj3R3o1-(o5XiU#%0fLE^=#2lx(Kiq77MMr(%!pB$Jm6eqZ4GpEGr7+v(#zu#` zB&8TQ182xtU0q$#^74ke8lGBOzusT7x)nm>AeG3*%xqG9c>R;+&eo8f8@{rgKZ$~v zyV}vA9a;Fq1mQ)I(ITtp`TNQtVj{lMfU-~*XmKyMgPpaMl$3)5PqH(oPoE|wB|Ude zo#9+btyu)YS(DX=WDqc%2JUw-m(rU55e}{jFo9C%F&2b%I!#Oj$?4D{@KJ|W$zyAA z{LvX-F;UA-Tnb1|PPQ6)9~l`L5i$9t+n3zY&yPviS__SSK2L$r*Q+?$0{j;=X?Z` z39C1 z>Z2+m)yPB7_Of)tEREe6I=A1K$%F*fgn0VgkKtFw8xLGJ4+giH-8RB zjL*u-`uGuL6Eoxi%}VZg6SeQ$+LxH}9ykYQrKKn4IZj zUV@sM+A4F4Ko&}xy#b7RhJ)kHYvn3Ns?S0sAupiZn$OI0UcQ`hvo?U;m94;Y_Ya(G zi?-RWM9udE$U{Af%7_6-Lgw4USK`G313b5AU6M~CdXxk!&rztUwPrcAfFOq}QmbB; zK>~yqec0a${%^(noe3e*E$8N@qcaGt1G=HH@e@?m1j)i`UL5(PPcOyhf!!0p4A6;PU^1*4&5HRUiXgMK zVmRp0qUKWvTUU9wHB@9Hg_&1VQNbl0M7dFoU}brk@%(w^Aq^1`5lzj}JOt@I`xE4`Q_+Cy6I z&Yi${ZjsDFp(|Hl*YF2*Zr^UXHX2sf@d1lh8w00)cqlve$s#+-aqg;Gr}fDk$*Zx8q7(4FxSU-kefA>mzh_4M*;|N; ziOIl^0d`D|lS@!gP*qj6sA#4mR>b<3vJDF1CnsVIqL&i){A05A=K>iA$F6>fE#zy> z_%lfRGwluKl#IKhK86n-Ftf=0%;MnTQD-=f+{eY+{{!j$$FuzJUno`2CaVH23fBo0 zO+!j6dIU5`fS<@m9^K+Bt5H4oEef)M5+&Ep3kdd<%IdAYPr<~r(46PBP#dNxs?Kt9 zaymLW6`4X7OyDL&N&@`+@#rOrH;3%#{9kJ1Qc$>ZhKfc3=D{iRC9uGd9JM<+^oXkdWzjSIX?RGpviIR5oTJapS5<(q)S$;rvB zY)7CJK^w1T3aFMAT6Ah&Jy9S?ppsm^`4S#uifXDNortX(1B}9E7U1Ut1P^TN_IG)v zDyY}~q)dzw&WqT(XC<~H8RFDP71RlPiklRkw{G1EW0Wd14Q~8XQgdhn=R`?*cdG6g zVc#iytbN=Y{EdIq!lR@86|3s}hjo7NnarP-D^DaJ3jFh99M~{C`V$S*&d`eK%*4c> z?V97m^wfl)fB-FEg$7*_`T$kOzvU!?O({Xh;)kwXRa8WTvCsF84jUkgxHNNyt2|1+ zS|jwqs*a_4H@*T$-|oI=;IVYA@9F)okTPdxw&pW)&CJb9tooPp72pR}poJ%45Km#I zLHs8p>3O5$<7G^4&%fnK*DX$|-uMA^I!XNX9h%<;8&tuc{#E8roH#L9Vz&b5eR>)L z&0+?ilBiD{rdGwzrvY?8TW#cW(?@z+(n_^BNf|#k$~>P(!{-~I0Kx9 zqXPy54Yh;A8B$VB2Exudpuyzj<+~tOnnGx!s~E}1u3Gk{!DHx(6I%~tOHh7uPF5cq zIedO=v-aj}z1VC3|+)UCpheSZtA|xa{^wiZcec5w!a{w_z+s|Q}AT)7=oA=k}XBZh7zkmM@ zfnRQ|`tKm7xxc>zBEHCGq#To3QL)cMKw#hIy%sqf)r^OT)SV~bLxul+0q#<|d-q#c z7j-y50icWE#4H>$!#SrHG`nx33v4b#!L)A{V*AwjZjM!dj^>E_Vd@l|aK!{k?V+Ff z`T4Q!gzsJe=|IlN*hzWGw@hq4ayUWCBf(?2?7d588V4jYF0QulGCXtREWetJfPerg zk1aJdwKS3wwRrTHySxB3#Kg*qvqu0RHc5R~u7OieN+%3W)eSiHB$KaI9RVHv3(!a^ zw&G#F0b1&l&`l>Hkd=XRM&HOz&{c1SoQ6$q1}F=NITX%4^f1tBo;`EsyIJ4N^fU`A zYg1$6=JJRHFK;Mdj7_}MsriavYNZ6pG92-kjiWCmoYi7m+#rdR6;LXZ!ftCi?5uW1 zu(?Ne1|HXY5)%||rawVs^`Qn3=&?<)TU z1+?_0e+|bg@%7U*gOhphf$ArDQnK)1YW3!;87nnot9_bui(#ZPn0dWhrx-{Dt-6!; zv|AGrHmk6L0?z*f#h;MshK2GCa^|dmEFC&vz;guv_;9HOsoTV zz+VC1m1(?7L80f2@TOGWuIQNW{SyT&cGWS6D!JBy@D0qG6& z^Aj}hBxVp-)3y4UHIQ%G*2^j&4O|rTdJDf?V7AD}%3Qy1H&S!#7@B1G*}<^yt{MaHZPsy7 zX^AQoi8HjdUH(gj>pXbyKu0IToUvET(C|Bvv`o$7?rf?;?B}Y*>FH}uGxs$$J4FpU z)o^gIbTU-pzu0aW85yZ{S((PhXwRo@ zz&F~1AWCT zXX4c>Ro2Y`Vh4b#}ze5X4+TeFy{KZ$jb3?cD+qXAk2e)h( zoy0wyzs@EG1IO1=Y%}J4;Lw#IrEdNJXhTE_Xj7&1+)C6)dEj3{FVu}k8&xvGiuzIE zB2{x;lz?=-rL7H=gPyl^z$HpaU*LaD#gk!YYdrodSiG58;$357@O1_WVUlw;^cM$3 z*;|f6n2@x@9i9tjj|y5LbCL$fN$x^<2I8itSjcX)QaaHeUs*l&=Yt>wp>xK00mNPHLA&YOG|-;g+?|C6#gZJ)x0BtqIXEPQ+ije{!G)Ua z&ixcD<(C1Imh&-!^5~SpUaIp+qNK?0P`_Nd+{y=s25X~L8OaE6ZbYMa(D#*i?K|&V znwoz4%8wwY30BTC`5deT2UFimRqR`c&gng&+~u|LT^>-y;2&Vi`d(i5f|AGwG?AiW zXImRNt_DLpzQ41RMzYyvEgUPkm&UEHQ?)$tQHTUfANQnk?6r1*8NO~PMGBGb5 zOq@Hb;QD62p<#t{cCoO-M8oQ;eNo^Ovq--&=pf-SLJw|LFtEs1{@-X9d)ufejm#p$ z8t=jPjYe*uMU+`8X761aWmCA2t*Mc#!SE>o`j<&NA6+la3V@Rzo)L1m5 zKh1zOl0#w|WiL%^ZVhx@ogE#NWYFI&rFL!rvjA8I`nCA0#!aTV@@RAj(93`bA5`wD4W>W@P zuE!@MbLRAGU*CC;L8O$VQ+)IhE69KEt_qliJXy*RzY7|P?gns_VL?IPT3Tp8$A0`c z&Doh5veVV84miXKjF)9%5=m|A)lobz;91!2zr&7?lQA$_WcClE(6+qQQ$mXPuMnFqUR6qR*C-Vo0kDX1fXG^w{9`YcuN8^ zCnIBs2BiPBzW(OTo9gnJR7^}VcB7K}md?(_U-?Vj($g#KYP?)<#m&u4fttG%wC^Gw zXcPbg8Won_?VRi(73Uo;jtB< zcdy?c{3A~}GAC2~Wm$m}8XX-)T{wTf8~X6RK<}-3Qkkr__QK!V$fi0TS1cZ_ z;4K_ZSuZf#tPm=m_GgmTD*SzXv~o;wD^oK>a@6j&g+;pE!-t!bZ0M9P%F9)ng=h0h zm!j$pcP6)f{8%b0L!p{O&jp92Vi~&LiBz+5aTUmz6_%E2K@WcUW=4q$KC7}e$J3|z zruR1I1IPei4@TvqCKeW~@@Pzb--TRl-dD&>OFQr4(%PDC;Bu{D%2Y2G`h+7-vhbQC z%=#M^9GKbqGZ`jx)L(O?Zu&*nEMPMFk)7vvP;o|N4ld#NWfm3z*Dpu3Nr_(j20I%) z)CS(#XdnV~{Yl2Ww&UmTvzq)`^S*ia-aTU@BVL0_b`${ajLBNE&iJ;u!E3WonrQUa zirvPFPbSHeOgJSktJdTk9P+Hkc!3{6{h0`MoU@a7rz-a%OOIoaIQRMUiSc@MLIJk= z&}}+L;PS=V-UbTaECB6=Ui9{4(?^AF4t8aRMm!j*#&t57%X~<@8?w&Qb_#ho@aLd& zCNDb^&MF%u#Ku;XH+9qlCY}dEj6d8WjbSB};{2sSp4SEWCTm807FcCd)0mi;tt+xx zT3U(HUcId4XI`?qE;Ui>Ta}H=>W}Rh0ns}qWX;Z_8x+9IT)H@5EW`Kv>&Y1Rg@(Y8 zH*aJZMU|`eTG=3e&H+mQDUO-nfSO^rD>Gm%yfpCmqNb*%j7&{=*Stw@==OVZD|Exy z*zvbVIC2vz;l3p_?CuuZg47&bTX+iAB#CxO#Xf#Rz=5{_GUzghBjeA{;4Mx_Gy^4fvd1$(_Nn($d?H9&HY- zym(PO(=q*no@t0J@wQszW1kWRapvi1fhUhQk2Ta;!o&O4&1&`)FAD9rMFdE$-`|Yc zH*RL;x);hoPL`+BXnD9>|D$I*@il*G2y^X_4?Dw zt@j7(i73>aJC1uBU9orXP4oGL<4ivnE~I=}omHTu6nNyfvlkTLd0bUmdQjPF7V$l~ zu|9dv>i(VZdttd$WO@Bs58imLKBbprM?KKvZfunCc`0aa(w8xxZmF&9xV6|7>*zSO z*)Jy%6_vmX3BZn%_idSzX(JI?IFnE0S!(Z5lcKw@z&_VLJLtEE5~TV!ZTR?(6T9pu zo@Z%ak#GqDyT)C`&aT+}e88i-T|OB#OHolWiIE)KdL_2;UNRl;uLV~;p_K!82KBh; z*%81!7m$3M7^$Np*c9yzi`yWCo)MCc7FtAjZ9XMrJWWJ&AskFoE-u@(wX#rg?@KyP zweYd6uB9K1b8toc_;F0D+9j>)^{g$3XtyUZ!{4AIad&-ox6bEpI(v=#aciR%`i~FPttzBlmEXLXY6#3++-r}SpDr->F{pOzhYHx|B4r*b zWaEw-h(%aqYcc&OBe#5^0iB+Pi#1dIQu{h@j`-4we&_WNFl#`T*2hF#;Js#^$5{wg zD8Tk1At7xBq-W2@LJN!;-9dqQBQv);!^XI?BERzE$BC)d zh%Fhvo}M1vV(U539j8mASXe#*6|Q+$%*Mik+rG^q1Hm^5V_RUg0^-+jbI3MN>l6un zn%CvD)=QIBGc~RqiK~FSf#+Y6&rC^SW@BR$71cW0POQ$z;JkLN?k!Ho0g!pPJFBq@D#01;+ z>eVYSC28Kky}m*w*T^*(NxJ?FgfF{aU#??rA)$4yguquit%B0MA3$~`$8+bo86YfR zMnl5F`o$=H%WMCB!*6?IW4PMO!)E*QC@8l#Nu|`F-EzukqX(9>YJ0TD);dKQsB~bS zDvNZdIE@MhLK&spfmG#ezk~vX5AgI>>~Q*L4ye;D;a4`tFmiHo*49~v z18na3I@!QyCcMejS;~>w!Qz{9xaQH;PBaFCD8{B3x6>gq7ANT@0H!rC3xb9YaY-xT zoYiQgooDo7KBO1~ODqn5sHJ{6oDr!RD0N64B!<-<9UWN;Li7IVUUc=&L~w@cuX4=7 z7B~RY702(#83AesMgt!8&k}uVDyqGmonClDL9w&9$2qX6fFm0J&8sEW{pisn=&8#( z{gprjYm!JJ{TF!qvw9{;`e7x$I6p8pu#@TvYIsS4O%lqf#10I&Amjc|XGCsvMMXv9 zM{mWCjYbYGE^NH}62Zt(P+cb~JaG_$^V<-#)U+ z9`worTM7cBHl<%*U%m-q@D|t@tjDSk`!o$P7)-o~T|%+N?w{%9K-uN3#lhBaXq>n=nFK+&2z#3u;G-2177<1>P!A(^XcF^wP*vu{k^EeM=zYMcqCi5M7)EXiD_kededL#;KPjTPo5_57HB%jUz^*_ zrh#J%;i#{!{wdoQ$Qsz(z2-E4Ok9^PUE<|cv?2c83O~z2c=AOv-m?%dHii|713tQ| zMrdfqz=@1sQw5+Dq3ntu0OkM6xr=hfg&oC-? z;UN1l{LuaTA%AGDy@CcF$Xqw#v`_1G*Zc6Uk2};}G9J6rEHR6>bL&Q_SX>+U}lJ^ zt*r&Qobc4CQv?LCufSEs8}~*gBpiVrShF=~HJFh>qb%3=cLIHXC*amBL138K%;>?* z&i?pYB*#kSdgqr%^=N3s4AIbhnjF*Efb$`j_k@X=nPQ^!VRp&rpRKK&-SUL1L}$;| z93Sn2&0}Y0=lRQ*im}4kAOg=4gNjQXF6K1zWWLt|To16-kVXcwvAc-W1V=RY?ZG3j>nV_}CnVp?{K@5Nt zcf8ij)bCc~3}$EG{s980e!dSZV~5L_mk9U9UM$Mbmz}sLIK!-lrY%*xF+8sSK_C6qSqE25e87WzA2P9+G8*s9UVov(yO)j zFh!W0Hqg-tV?+P}$mV*}BMze=h1(w;>;mX;Z~8J+>R=3M9g5UqE_)!}iwzl3cF9ZZ zk|Z>b79}=aXLO{t}CQ2uXgaLtFvq zngCaKcTVtu00{yb#P-@+Ud_iitP8N=t!eVHFY*;^tn2gIHu26db$&VcOy2As-`RXW;31Fn$lR|9cTTZ%9WFpbJ!2 zu3lXO^)F*rz8-WEz*6w6GiOV_0*C@SCeiax;9(S`14t5tgoF#U^~*5*x#2RW4$a60 z-y;e}T+C1~5{B?^#d}m{zOS}336tJ?3{EYeIg6m_21o`u4`fT7AqqM=t2`|#pclb^ z6>aD}(AS57UlT`XxT%NOUIiW*hZXAU`yqiqx`o>U5bm>yJ_cRH?_s2}ii$8e2ap8p zD%W0un;e|Rn_-U%1~|F6Muvyq^1+A;8Q%O!g3svbHzY5mNs**6>~~6!a_nTbbaf?4 zx;cVd3G#I}=wva1=5He-f0g8D<)N6FanvH@*n_DErH#NR1CXz7xx4Ryo(sC+G6;`= zC92n&=%Cij62AvA(Ph53wZ8re98^HXfLiYw!I=zdfr#*#W@V9~>5jQhLL$y|zL$}Z zBGk*AX@yrfR$&vaJRr4oqJ;l{xpq_50Zbz1)NmyvD>&Nl_%lcx*YB^gAVhP5e)JX8 zx2mMMGw`u*D?ydUx%0u}6Zw7iQtgK(T9wt#r2+Fqbp|}&yLxQnUUn-vvDZT*DQN~gBf-bCJcZQu3%`{*{#II(Z_cD1n8(~ z;3l{K1srKGn=Ji`El6o+sd*{FbMz~~3{aGkVgf}T5&+v6 zAi>|i8Jqy0vfqNpz$m$cyHW#P7|^wweVR*5VsPsNP$hh*f2{YYc(f%YT?4o64Mqyg7M-Z+n%@J`&(LBIfgAZhUyGa0!O=|ZUu#H*bqA{>daRn zJ?{|g8a#RiO#WjI%?$8tV`Y4fU@PFbyKl}B^S4$w7 z9&hIQWM{j(x_Tb0H3Nw$;W%|0rU6E^!z&d}<{s$j#l^*W-jrDS@hSLF5S*ZSVCI5` z&&(GxYGicu_FPZ%(KqNdAf=~OW%q-o^x`ZvU>)3TrF+r5HXsOHtEQ?o!A;eZuS684 z3ic`N+#Jp}Om^iKRBs8KFJ!0~9<5dGUO3dQ(x!|Kvx+8zKkMsq4Q*`;`TcMQ0L^^q zQteStaaNW0p~U9qW>HNVg8;NMIm6$tgh~Fuz5zfMbX!B??G_$v02dDr*a~PNiCwEc zU+ekLUx$t2Vt1YGhOq26@ZjO#m;AE2&*81Nzy_9714D?BktJpg5Hq{(1*~x)5S^5I*9Ev;hK{wPTy}t<8 zk)UtG-K5IVDVhgs4XNa+EQs;)uPzpWO$NLM@MxjTEiILdB?y=XH98v0NCIsNIdMy{ zOEMxdasw`sW8dFJ1FWeBRIdzKcM4S#dNJLKHYSI5bHK$R{|sr)h9-PnM9L zg5R7SNPGd`kZmN??cw0i+aCWt>9Wo-fP;GjCkYW!a!J~oH~#Qy@+sm{IV_A8iB{%G z>8E`kC{yRtQlfX~SYuP}!JW9MF>#QNn4G}U&s|e|P#jEAx*+n;=@QRPPrd$V{*cgr z!PaZnq&qNIU;q5e`-OV;dMu%rA0QCOiKx&X9QdwuuXY5!*+CKIAP~xc|Nq1PZ`owG?GdU{#g~{0I2lDnUQ2(Th={ZWd7DntUyhB9{r&rQZOs7A&CP9MVnRbh zWBqLH$&)9mUE$Uvza(Ax$uKd87jI5_S#%rSl(SG51_!OXBdI2fb@32{+&(EN43CdH zF1PwIGc&KB(4Xy2I@{UN$7Yk0lV6MP?Cj__R!BEIZ?@l*b#`_}MGa2obH2X03Hbha zXw}+vbajc#$n^L3@BfXU&^M+eC(n~M&_`^xm@X3y#_dn!Djst+Hr`xr3mA-NH@?5U zA|WBkFrWaZg@DsYad5_S6|>94=I7_z+S(Y{FfcI{3~8_tt@ZV>UcQu{*c#7Qm64J8 z*9)33Y-;H^0Z%t+Y3XnhenuuHM+XPEU%!5R4+tRTb@Yp2G!k8kKkPAHx-xe2zO;cLm2tVSPoa(MRZM)3yaSy8hge`wa+&ip)4F! z+XJJ4DqfDt`W0q5Rr1KZZ(-agCm;KHoY4z4gjDhB>Mr-#=NdJy9m~o}Y{)`FKAS*y z##!~AE;P{bZyqeXoiuBC#(`q7PuUq%>`h`*`=&UpFID>)a@pE&rg&Ev@ob~}#ld|2 z-R&(o1x4?1w`P@v>;6ooorf?vIZs3)gT*|($654Hc6fzZw)^dRP__& z7a!t#-d6=i*{pxerDLIHCe5pr(>yCIx}Htg$(4E+Yk0`7WhZ>SK6{WyLqtg3_xjXm zV$wPRzstMGj`M6TG}1nXi~n;L$?5EgaX1NX(tIk+H94SYXVM-k=#ld|9Byc-YX)cnJ+p_83D0xXp zZ|~41Z7$N}{dt&4W-FaQFZL|=ghri4~6A|G4)@S=GGXh)b=9EWDGFZ11e;`l%y6$8FAf);o&bJ@l;yJVV88 zI^8I1sx{JFeSYtLdR_4r{v|P5AVQYF$X+nZTRTdGCs@dE2+uOyFGddg-CdSuI709h z99-Co&V@oGXZ@{ii)EQnji^iY0YKJ zNu-&t+osf96H{`Tk#m2*Oh2hW`$oT7Gf7KJiETc|bql1u6IT|A{}|SkG4G#>gE~w6 zA(1m~H@ae=IUXkYG}GMCS%!f?Mx8*rfK<0xQ|T&K($xpjG6bnzOWJuBWENQ+$YO_r z!Cr>rsI(Y!jp_LAbo6xWS3s>#WS-qNSF19Po&71!uVzi}#ynY}Cx}!elt@`Lw`UC* zzupK|!%)|8Xldys7al$nQ1&!^qphf@sIr(QcS0f0z`(_&hHa6Vm?9Gr`ht3t~A^T7Usf7YTc#*EXf4!f4yR2RL;;U9{ujgWHg`?jNu z$$>$<3Y-z_LB!~!A8z|iRh0XdPQw_Ewo*08A5DbL`Sdn;lH(?*Ke%x(3_r}q&@4*c zu+>%xYfbs7ZKDzpM0;F(0J&jFsAanL@(5?5DSsq8rUfo%>6>Esj#7-$#jXpaFO?_QssBru4t@RK0 z;x&f+b>p~I^$OwP*Gye*-I9pVFlL_3{uWe?@nGvO=F^3p!@jOAw`y+Gn@53Btf871 zg{m$&Xlh-+Z5ANx&qm`+HaAUXDkfCu*)AyF5WS6U0t<)#&*e9!_07;5qS8 z?X1V1|FeRdHi*)KWM({Ms5E|WAca-61Y5y#@ z;R})kqV4s5@Z5`LjuD0YHGD??20vB+T5<6 zprJV5Yq*8vl-o@>rO)-8O|*1B|67Xmwh`&0Gh&t|mXpHO!n(=BDE4iWuj` z;x13gNl4K04h#&C3XFVU;$SY?f_w~0D~`YVTaj5}o%#}w2H`nsrkLx_YfawfOS77a zcP%q-4eL6?1MeEYyPmdci+dH>+J4AJNdNweGClN+!8m$%fumX19M+is@q>M^Diw1p zHiD4L-TQj))f*0+m3gJ22}Nj0;W^{dQ17R!!(`8O4R$5fNPVWjAq9!V$ulPvF|S;n zn+w#l9pv#mkp)s;WmPSXkIngMoD?9$({|Q`SOx z*@!tvu4NzMdMjLU?iTf0up?f7lzt@^N7nYFT~hpw)MjCggHMR8>qQsxnodpacCV%*%TcF0`)(aCT_toc6m22^za}}KuPnCT zljf!Y=6sQ-+l(uQ-Z!S~ou745<#g7EBW1tXR6eN#yMBJYAw!VB?P}B7#P52Qd%BbL7EPgr9Dklp4wBN8(vE%V{-)(IrrR zqT?Bk$$XuoH#;E#!PTgo;N3;i%*>4WWDzls-I})t331=VqiOl2Yg7dVM_a>+XU~{Y_EBrn4NTQxX%+7r4#Rv3 zs^Jso6q)|=o{b!c#U11-mQ;Ds*vY?CDHN7KGMj)*7F}*6vI&|n^J`y=yxiN{bK0LV zG&k2gW?{mv?v{yX6IFOTjvF%(5u7&Z-dH-r&CJxeVBTEh>(mVG?Q7OJ zbcOj|?#9OPaW7pg3T_2_sqaOKulZ|rzBdi}17u|6=|lT#m}x0I7UrE49i3Vkg09Cm zue%aGcssiV0yLBpcGB4#%pxtJP(mbh|4fRi{7@#>iy0BgwwUbZw%;-?H@I+C3rf4= z<4smOZTS1G$Qa?wjaEB}u2I=dr+7C%>A1Z~eA0*GDQP&$ z6DEt6i*Mqvkgw0$Izy+Lh*vusv*&Z!A_H-5ofkbf4-W+`rb_N^E`pQ&pj(3OM{wG8 zY3|N1i2RdY@689xB-$)LHJf03?R{1?Q$nGeRP+6}OfU)SmuE@D#EH%Vla~*J>Dt;? z{XUJ3RpDJ@x$;Wi)5$glWUfn#s+^8@E#{z-Q9lb*LQQ8Jq&u;3_q1g_`E7sqotK(cj@Esu^=HUOP#sSq2+A% zZy!Vg>H56#Q?P51vPT=Zqoa$iSlyL5s;o2#OKWQo`=LbX1zOjL{e44?d|KMm2QYMX zIv4Vcj}!tVIy!GzzvYa)KtzsTn+xAiBMRCkpb;465kGr|^rC&O+vNyODi$}4#x!oOFV9tfi0Fgdb($+b>AU#} z>$V3$%DZ5elF#2LT{2>mV^x>qw@fUo{IF>jTwQgm&QTfScgS zbuz;N+lEh5VQP5wVid=*(!ye>ipkF2_UeQpeA>u1GIG&wqMT-?qJ`A^4f#_eBXG~S zxcYa}l&Oy#jyEPz@R1EZLAnlq!wrF!px~jNMargo#n#XOdNOaiVWv=HlpNl{Uvtxz zdE(J`Tz4+?TH6Oqv<$|KI_!oSfw17tEVGIB@5rh5#_4%Mk58AX!7J43ZfWg%NAPet z&E7udOL@Cn1F!6ePA{rKwTT|HFe23xO{u|K7`FZwz1i%AUqmT7n>U5|SKEyup;^Ti5q)-?OB^#WFH}DUO^m!3Or)Hu^G+taTe_2u(FQ zM|gO7xt(6jR<*P~b#mgieVWK65=iu0sq~U247&BetSjc^G{@CjUG(F-U%CntY;SwP zt&i#UmR+TgK568S&Y8xhC6i}!YlUf!4smWSVpR45k~vSM-ZqW-I7*?0GED9RkF+7T<`^GC4<~LjqBQ@dh{FZt-lJ~YG=Bm$b2TX)d)shX*Fn@@i|}A)g_3DB?t+b zOrpDA1pSQ5v0!|BU-ob<+Uv?CkFGfgx%_i; zY4rl{pFB32aeU4;`-BKjm^9MV%Ofd{5jo8G2s31e7&SoR!o>C+&o^h7ZW<|LJbO(~ zZ)2DDrw{T#vgk{bK~KR)1phq4TKByNwB)UeTZXQnfuUp<6XTKmJmze5+aw;pyZ)73 zO??3ReN$1^nsrLF1g9;a1(iGJYXA1SJ09x|!Q%*hZ41%LR2KmWlGB3rH9&id9=`mx zMY<$*Oz*KVzl~dQyT~{6K}dri-|742)oq#m@)AvAVjefIFnP=VY){*B3-z*5Y1{po z*7~vWF`LV`O}AK@Re!pXOG~+((HFhLH~Kr!GDMCy;lgG;%Bpsp1za@@Qt!7RU0smu-iZZ4MhL)UxBkW9u=$jF99 zvAMYZ`v=sa;3rKCAoLd06z+S=8X!$!rlHVbPX2-C%J3d@8am= zGo=?7Un4O4`fM(k%_gJ|9i(Lx74}RZE1vR;Fd^QAA(=C4a{kovuUMKil{r%axkI&E{jc)oh6j;b${m=(eh!53q zI5@~xGk1oI>+5$585$-DhQ~0mdcA-By2Op6=?<4cH7Gh-!ecIjCZp%`mEG+*QA6W! z9mmS*@#&z+%quQyVO?EbQBg#8BQRB%idfw*iL&cgTuS)`PgDJbeo9{Z@n1P33Gj({>lvXC#YM*(NeJp)KtFvzr& zy9*618GpYyQv?erHK{OJEQGK&vNynAv*7X(Mf63^<9Gm0*2Kib&8==q2r6MPc>>Cc z%qU?zl^i-7x1){?HUn0~UR0pZ;NYMxC~u*mGL3F$4ULUY7`5xh{`^T($Q0w{2PB1K}rP;8_aMwbNp0>Zmr$=&UpUK#bI$? zULBxfx3;zxy>2YW^ORLoRF;>Qdt+H1?r#q9m^3Jzq{*=Y3Yj{M`5G#%(`QP^0Ol2b zO32yp)}Yw4EN^NUY7i^^>NboR(4SMo&EDQ{nPQ=4GZov#&G)=6N1x$Rz3*$|TXt;8@?g4$x%|LRoFG}FM;H2KF@v+wdA4{%P$!5C0sU zk(L&2C`&395fvAj903)Gui>t#N%Oqf%)~R68rxkE%-dS~?)$t}g65iq(&r5{IbA|X z=o^GQLW0sSV0d^qA|hgAVW#>NEVc#t9bIWr@@Wq&-d8D3y2$3+E{cbTG3 zH9lu%hI@(mt;LrptLsI|4UCSF5rQtRWZEs?65BC7&GYu{TNai>Qe$c=s$_r5QV916rHxj&(EKnoGfTyBrV-{ ze|OVM*x%c`O0$C5omO9BUF=G16OqC??4d|hOpYkE@LU`=OfS?DPRw@>ZwzksV79ur zK*93y*tbIwg{+H-YO}iUXcfBg$-A@wF#(^GQOZ2=BL^CYbY*?O=H}++f(O>!-3~g!@icYZ13C(pt75kg zF@)sT_-qPCe9i~2IW5$hz3w_giBzktw26s{;lNnH0&8k&va+%QUxD1A-$j`rm*M3( z*ARGhIzH!A+4uLu=?1skE1AR~;^3UfsKA>g5Q2gxLusz38;#Bm14%qx;iS#ZhYQ{w zI1A5#(uX76ezRn_Qofa$k-_V@qco9IU+?bX;v%|12B_RU&@C*qw^P1(^WoahE=#M? zZG3F3!iNG=!8Hi_r^(ah`jz9Z{?dmd_5!E?f4o|YT-;kd5Uq;YvW+gsE2o=-ik}{j zjLgiz6n;Ldxti&6W4UBr<;SH^>UE3B%q$@T)EsEIc)1biVoTXEy9FOFlf(8>J~l%k>>Me#iw z$cZ#H3dz^l$3s<3O%14;A3vhNf%qXJN<6v#O+b#-1p-L3@alJzZMJMlS!*D697q1a z1a5F0$$O-OK1|P#U>5x>{yExZg&^zHY{K~7Gk02Wgj~$L_N&k8obj==+@ni#TA*tw zK*R#2&rHsvtL4mE0$-0@c z%@>Fll=JaNA(I<un0i+?A37g3^I?2mj@J;E$p3HyxC$7 zC=YRv)Bcat@FH!BaKRb(`3zMS5O>{t|D01+*|9fmGE;j!l-{;G-sbyxFC{FR!8Q2x zFijU|GCVa_yHvjtbf6FxXHN|G@$dA2#NbhQDFZ`>4ZaJoL+yWTzb~ z|H>*Oq584RvwDY??a`7-96kaC<$s1z5E6p7_^^|K|JFFb(c}HWg5p-ma-(-D4*QHDBBrGd>;5~ zE2M}yMn-BTVpwiW3xZJccmYxKv1;5_oK!(3*D|5@;cOH~(T@;F$K!lNm?h*1OqzZ# z(ze~kIysGq_(tmO(KIkQq$LATM!OG9M!4Y@Uo#`#dfxjW>!xpE=pmt}a;cgtOS9f3 zkCh8r-HkcyU1E+`=WOj%*58P}dOE^}<9(&9Zx)oR*bycQ(Q3gS_8zWb6&4xJ!?`VH z;P<3fFGCOy7mh}CAt(3)X5)^AkPu%`XYVBhh6m2Z#Ox-+eUFJ5@ml{s|G?bv=J~&S z0sj50=wy%^c3B97H0W_)Hy!P1%9&$ih)Ru1SsNU3_WnmB(N)^#>6Vu#jMK zHKKBmL+9-}jb_9%ASA$1V)f!=?lOqY@tkOZ2MrT?e41W7Jo4u9o5_kOD0vAVc?NX@ ze*&%w5!_5;s_(#Rc~{xT6iM06`JcR1Ok*8EV3c7X!VbCIZT4D8&%$5#vAhon6HBcs z4%i_Xv7VKm`wW4R8B{FU+uH*<7vv696qNP>VI@(XqSs5acHiYm79ooJh3vdy z*1E1ExLAQ;QnpgNAuT{+9Z2BV9Ta>(_z@MwZoeT-Oq}d?wgc4G$sC0*>pjs0uH=}s z4Ar$YurTW}-*3HPl{zF3cU9`FtI-KOP~u>y091O8n#+Ep42wh7F;|Ha39Mxi~v(1R4K$wd+w}hA74Vg|zVT@;-b1yjZPN|L*$C*47rz z>;B3Z4v6x}ypA(9wqL-9w*Tv7__thK4_BK2j$Lg6LHzRaa&ArwR8nkgY*A5B5Kfqw z$^sx`egcwuG4q`r>Z8I>%cgqdy}q~i8XrF{HufT%{q*DnOl2=2Qv);{x`0U-)D^C0 z7DO%c2aTPa=Myr@!szJe^z@a8kxspn<v?gOl4>Ngd^Wt@XIG=z* z6>GQmx$>Fi{5hJ?=x_Zl@{AJ&WE8KxsZvtz(TuX6EMhr>Qv-DibTXUCy+O ziIpGSrN2ljDfz0~Ym!>ITS)z+dx6%0E)H4$P=zAE!NI}C#>T@#@HHhRrLwXTfTXE( z=$ibN93H0;dM#dNW-$hw8s8m?6BqAAbCEFsK%RN1ct@3ivFKW5s^}y#M0+XAKMbShDI59&(8Pz6(!H=$1yqZFr$e;!4(%3b=(=tEl_ZF zzBoEM3hEP~ylA@Ye5Y1BV+%T zOC`a^9tG<{N=oVxxG z$>V_nye3engrww?wCrr-FJHcxn$kJFdt6H5Dkx1qLjO4)B|9iRe|$|)3iwYC9w<0a z33O1H)?(khp+~W`)!%$WjtT0_+E@(JTTZb6bd!k}6BFD{(t95uO6dl>a z^`Xh*1)0^3hL=}W%6daCrupufxcJ(}$^w61+>EP+jU3WbO03dTOa)QU5qA?Cu1_N? zQcxUkDGWPzh>B-s2S8Nud4l$3l?;!-KwP6NKi6@L+Rx^E*yd_Uh+obg3@=uUCo^Lv z&yxqGr#o?TQ15fj=eNJ{B4aY$of|xZ`P~yNPs9F(G0}c~`k*YQ#BucCNbZ`M>AR*S zv;0>95aw=YLj?p@^KPqnY3X!&_1Mb9k5?iR5s6Fu;;l65MI$RAXQc(`710->5S+q3 zAM5eFaj-O=zrat+&)1Cj2Hk3&Z6KJhTRzZs*bl0Rg zG;*AroaZpudyvS5goH#pfu7Uajoi`#w1RTuFIRSP3E-i`#XaV2@Fs>!-{%-c1_qa) zgS);y-x-&Zle0crY5@zfvaBrcuLlHX`vPrwpUDY{9BXh<+b%OVce+@&8Ms+GySg5Z z+AzXN1+At^(u0G0Nxg4d5cxqtX1qx=H*0HcwOVL!y*aluG&I!Ko(HWnSnT}n=X;6) zk5wxC)ESW!53kzo%uG{L6J%n+hgj=V0Wy1$o@@!y3y+2u8V9zfu7PMVym-`U_6_EafRiq;dn(cQVdShrSdRf-^o3pKXB zizwxXTM@#jV}dPBa%c{W^_uTX-`laVR(kTqN>J z!92!sCYisv+Ak98>-E{P27i0^m!YrJc@}^-TEXR4pt9$!y^=Tgt}ywwKWn>!u+h&> zI$Rp9I1oKcs4q@IPSA?;UV$nqo@)3wB|J9Nb`M{tS!}hl#!9vxsv9a10~RcO zyjbYcVYteuTZZvqE4!;JpxW7upox!n^Uq)du}j@JWK)@W`~Lo^u`Mb6R#=!*3YF*t zuvn@4omB|G`_eihBH~yq_uQCky4SDfT!5LSn|GH4rp(OJlgI+2{bS;cS`;d62wWrY6uGI0j*U48qvhjJQq&3Y%;IImOKRoMXAqYDgFc+ z$cw{8-3pUoe>9>_(C^tS3!{Z*$jv=kP|}2a3AwBot=?w@#l!}%8^6ST_T@{3H3mF< zDLkgaItBJ&i@AwOfm$i~+qXwO(F~!49K5zG-=_kh?bujP-Q5{=!;A+3b3!dsVxK(~ z56wA=cq9Zs%!9%ulgR1n=H_%T=L%qjJ`H)E&Ghuq(VvwLeD0DhX*K+rYH{CZCKJ+Q z1`wt0zwIbGABBYUK;_0D6?lb^4we*nu8(h1+Z+HiKm2ol|CJ*>gl8*-erUgC+Nvx8 zd6)4~%wOf`6UYsKgt|CA<+E9OEZS4HQ2L?K?7{w-klPm#5hazPTCcINN=z%vvvZ>U z?QOWIP4d;itmFtXFiHkWbhh!5SB?210nieHg5Yenv94>m8rIBIb@wbFO0|C{aPQoC z`1tS^P#(gWbsE$fEeHq*Bqb$Nd6|a2Z6#Oj*F)BL<_~eI_~0gu;#WDUtXn5DHXf^p zlj_kYFr-J10b99ty%R_&%;sZ#eFXr9uEz%F-1>Exa_y%1$`-O(;3p-$Bnl?=wpef& znMtYxah&n-H`c4$snRck9&C%{*&k&zi6q@(YCMSq1g4rvE6wz7&fDXvcmn{F2#F50 z!0|9MGczzSu(M-(Y;n!a&6${(N<}c0lhBFZ8lq&NO7}J06dg%>33%{+zO=Emf1G3n z#KXj3X~2|05Rai2SA}YWqf_1}M#Qjo#jj+9-ht;$dl$*h&14q$HTz6R|3k6o~UjD+oa&IvSN58v5$VzV4pa|09LjiJLr<3-H@{w7k4grz%XZ zPxWEdoqjY493qWo#Wa1IZV{8yfbLmENy+g60L0EKk&LDpVqr0yi^uYe z1^));`H8UnON1xr=$MimueeTSUC8OZ=kO`>xJQ5MwIe~YRmB6)x29OB^E-#q7cJw&QU5FF!Hgtavw7NHS-AQv*BrIzcgFa_NJxhrVj^SH5)f6H~nu>zaHn?D2S zb#?m2z?a$qUNPb2wXwzILjaUX)0y~x;yGZ*YwLl@U!i%tAeUdco<75-{___pAptN% z7=k6+ff9_2E*W-_5#Km1-$X185E~c*KiF&jAWRBCB14p40hNTF1_)RK=)046#4|Ik zgBWMK6tqa_=o!=$8F@gfQIdW7bg16hgchAt%z|-eye&;AkE71XuUYU>kstx(r$^H3 zeuN^F&k3SINXXIZ``ly8E+o^eB#5RRD z1&&~I;A1Mw01y@(&d<|GL;q0FY|C0!2KR=!=qlj}Ge_Hfh&4R`x_LnWTCT4eMh0E& zSidAh@$?)k3iU@*X!c6J?_XVcj6G$n$DxSt&p8g}*&Vg`WG$N{Ydg|1ZqNCNTW54ldo|r=H&&5yN zje(BC=3oF|w(h{gB3O~@TyiaM#$uR-PQLNMG z23S_mrF!0+0|&B`jSaWgZ4GEygQ%z&@Od9Pw@06oY zbTo)RXwMq}54zghfqCTj(Iiki@>L3-4uP&e7>@-ANkN36NN}%T%NG*HxaJR}@XPDH zfdhJ`AL{Zd1!#IQ*gzAO7#BC!;971uTQ$NCI$UW9iOaq6K>%T>^xtL(V5l*kf^en=k>aE0n7~WqX!fJ zK%%AwrlrXuGD-yCKRo`{)C5opw1*lSQP6Ux0&HpkxZbuMis^@l@j*escC|AU*zptC zjl*h=b3RJQbo4~mI3TO5gSag>ZkxL;F4w3O&R4-91ug^KsKDU8S2#Ga%-RWYacTf~ zeG&ovMi>-|go1)ag|Cdm%34w`@(C-3rE8=AEr)2iQH_WC9$=tJ{{bkT`uLv{yEtnrHK*1hi~7$xw^VOaG!z->0O?h z8XX-S8TqMbSVAZ1JsraQd^H{(%gVkcH?Mv>lzE3B^Y~G2X$PbRj1`CHuv#yuG~v8FtbudP2dZJyO_=`D$q!Sy?$rAe6)0=F@zUj@dulrobmE+RM@} zU0zGlU;{>&;bQf#oE+K+i5|LF-O<*L6B@e>Gdspwm(LBR-ddat;Io<^;7Cb5@<5)H z|Lc7t`T^uu&TkHbPnR1hmv{BPz{4RUXAWX$rZsGe=KU`HJ+%zzu%Loh4XLS%kyKJe zy3Nc;KCoXI85#NcyTA&nt)&I&6SHFx5NR_|(5EAt5*vD1A7vYG*yRVzxn%& z5UX%Hqe$ zIxlWJH91{x2Z_G3eTH_9@#5jtqe2=$9->4&=To=|ivaEh$-16cSP=i=daeG<}DS}pu_x}cj0Ffc&3utD`*OA2`Iu(RVn z#owsBT3#_VW-i8GKjuDC;s+E`&}Wce|E%$$agB&Z#KEDZkKrv$`?}h7+=h`y;dCAt zwDZN#Pa*TITLmp0H1uh;rzegVr=kmz*=5@!d2bu^d3_|xsAy!uaAxWG}f!RlR>}Mox{`f zS61#s4y%#KLxNQ>uHKe&v9|#KuF50}+x-EoYI-gg9aqStZ)fI@ zt~F!x_hoyB(95O0hRxngD}tWK9v+JTHN{|Ick<{$+>GABa`|Xrp*{#jrLc7qwXRrW zyisvCK3{Z|Wm~1Wz9AKPBo08Y(^X#l3Ge8pYWa2x1h3CTE)Iann@7L^$!Whe zhN+}t++p_@Yz5J3-X}g@%{CBffT(&z+2-Nd#+8&{bOUQuKl6DSmw*Q&2m()A zN@OnE%`l7<0}YLqHqB5f=~A(b-rQrSG5dd)*4?UD{i$mco+;u$^NEmx<|kniyXGt~hxD%~wH$&O)* zB~_-2BoyMps-O@w-JU5y?SA-TKj?jgz%bN~C=(jK5)n0XSF(ONimE1-QP#!ybQ0f` znJ!k$Jc4ZBtZe6aHd1DT@8%I@=C-n*kQ$!rXlzGHrLB#rh>fX?u7w;S5gEnL*TAcC zvpcWar^xm*Y$5L~r{jk#yxE%QD_t_{FTSXFHfzegJL2BMjUELLS*5Ggbc2b&ngqJj6^-+RtS9Wk3ZIq%{)6wc%R8k|S;vK2*4=|Qi@YPIXoO4c0CkhGWApN6WKogl zIe^t;t5sWc;&-RGj+IKl5|I|Jxw8JUu~}B-lYTum<{h6_%V@U4I*QfY*95|(WIbb5 zJ4#h$DJrJiF^BPxUs^bKD1G7(%IX@QY&i;^E7?SAHwT8dT5mfoVmAM9%4J^nL#&D% z=5nh$wm*x!_trTf%f4qcnaJ`-G7W^Tp`)e<1Y+}7J2cn-y5NodLCf?Y&wsx~me(yw#XMjlLD9j%G3uB6&N1$NN-YQkX1zp`X%<~G9x zm1mupukqGkJU!zlWk#OXs~3vHuS7*J=Q$Hr*LBSOewTh+2~J|3pZTqQTj~G%w@EG- z1n?&bOtkoNd%Xs*TY(}3_$e7QQWpaUKMk49su!-oSZD}juT}665#?#kJo))8Zbc3Q zwBEZLS?E9I;IJA1$1y!^IicF``yBFIQ#c+F)=Y}8XEYa(b4j% zstLD?w(ea0S8F(z)=i{d>>YoF`wf*mlV@hA#3%26ow(K)MCSveN>LH_YO%ti%hr?# zmznoC$=67%jXwIV>b$ceqQ!TwKL;lzF(isKqo+QMGB5PWLX?CIvX;O$#27|(ySr=0 z^mf$%I(quLx;kd9nxT#k(RK@Enak~LUpK=IRC;u3C?6R?Bc=V_ujigB)jwrx2(3@| zf9^$MeZkrOC*7s}VdNS~5e~=Vayb5Scn14@+bLIOIQ`zwU z#-j0)p|n~T1<^;6ef86({fNCA?9b1pz%NeqaA!~E=84X0jkWWQ^v>@rQ8w3d&FBR^ z5#G^WeOMoB*MjejG!_e2CcLWrr1Y~ck6v0@U>!pQczUV+!T~EoM=)M22QMI?sbvyw zFI0y|MrO(kd%e9=daNGC%r`!=R^sT!(RRkucisyuG_P%?23zy)BE-7=4JyLx;%uVi zPa5JgY8v2O9?k@2JfFUfd~>?yFr?mx_GwU@2sAZkmzGAnb5&l%*Q|j@ng>iBhM{%! z^#K6^ssMHs{k@AU-f5#di^TA-0Ou;>@*YOFP+2V@cJppgW9c4IEHub+DtB>4C3$YVaHqSS zl!lz#Cqzs1zwFlyj{`gAq7*ll=G~2Xi`0zE+*DU_SW&BE>M+o1Yi9|ID|Hthg8Ki? zkx@A3g;=baI$IufMk?Z7k6jyzXap~ZqF^L?jW1POPNRC4{*cYlQQRIUivCo5b7Ja* zR#-hEVsD>>D!+JGQZO?Oj~g(-8OAehA)Mded~7E{POeVg(YA|JTw$A%A3v2g;7{#ODFy}WD;eJ1DZ8oL zFTD4_Ea7}>pejF*crhqhb7~>pjE$OjvwwR6ImyKHPJZdm>F}#NRRtwnhH8*MPy5Bp z%k(aKI5l-d$t>Ze2o$qe(ym}+%Ma0_514c z;@a0{BUvfz6ak!*T2?uT(lgBJEyIi@M4?L{1j{Tqlf^_wj~8jz_pDRhYaA8pAp1T; z=H>M%EKJL$aoDGS@1=FEKoh%ArV>`B7g7XiFBQpUiRr@QG2BJo$||Pf3B;W?rz-#^Jgnx=_JJ+43_;HcJlNb2OsIY&ysX`@Ow}-mMB;3MODM~q&MzM zV3W{q!E=Re4TC9l0-dt^g_WRSGuRYw=W+50F*T8an;VJU2q`{ZG+<%Dr}lB9^X)vn`A02gKOcI)kMPCPcCEhFr z6L#p$g%**ZE?t3WMaBNHZuMNk&y*FL)z?7$1ZN+5%4mc+m7^K`>E|}-$ zJ`MsL6w8c5>Q6awLPGg7Da|lvI}Ok4f${|r5uA&IKv6obLrB8Uq}7w?j$f4>;coYz zYSqg}Sdq(RG@L#lPuy5!|uSAOzcR2RZ<2+@~E41)Tvm+X@k-x*Vi{`q3Ju@jn)=7 zQ_%x}WsKIMqIHpzlLIUTMtm;drL(hDCrrR!0bpbb77P?s#<5=bnGzL`tUSIV;^uhXCFB$!C zWy>t%mD|h}J6Pci-h9%8*;iNCWCGsx?OK)byRq>nYdXkUmCA$m8poBp#Mk5aHq6D{ z2-OPglkwwyO~jTZx@@h+L}vI6_QenijC>Iil?OpH9bTY+c0X4*M}rK5+G0K#wCIKY z;uRy4UqQhN^3~OY`Qg29M&5;Bn75jf-V1Pgn4q!yo|jby7cCi{KdUyH3~ z^5$i)6BXXDdH*12(R$k78hJ?{UmcI%_oDvi=-;>Za2(BAVuv1hwUzTYzh3AeUiGIO zlDan&`9v?)2R|3MxVkhS+_dLzaD6Cc889A9$?Qwc_qZKZHt*S)%7gJJi(uFpJ=At5 z1kYPO+?L*_eA9g2ELeLA8}if%5`=l|VQy~d3f^%l*)MyH^K442ek+gSz-#6?AFYHv zKYidFA9aIO7XSSAL{#Bk=ltx@qiLbxC{J!3gS6%~*e9{Bd9pcsx^TUpn&G(7btoYC z_U-LpLB8Tdv6+b7k+KZ{aLwaTavejGVx>Y}3tEycgtls$MQ|yL_Q?Y_&MT(w`I`m2;mn9?o3;;7Nx`h;f3bh-ZELtnS1VRFo?%{k+i z8G_ko{kcTOCY<*o#oISJkcQhD!~W;0$vO`UYu=soPWN1nuP!|q6B0U0*N4$3W*R1i^$%2aORwT?6;Y( zIqy05zVChV+!VyhKfjx;At$(d?O2Ztm zJCwUNRXk|!>J1tnXv&Y9@5$X53`Rx#?|3VJ<+0XqG%D+ful?%cYa!rYmt_4v``R~S zOPsE7!0~q*8|}+gS2f8vXI)I<%B}x9lt3VuT(W8k2n!3tVs#4)vhvc>(t3M)`!-&o zfpBnexSap-w(v>&yYY!-Ap1o_L)tnzHNolEDT#@@l$D`3%nO(h0MEz0y>Wo2iBrFo zHuO0C-Q3&>q*-z5Q}{YkusrTl%UA5VL-_3o1wqDfA2=Y0hQYWJlB2Hl36gyM!;1*= zrEAMBtVYeeTEq zjDjHFrYEZ)M=^grp+6tx5$n$Uns@7#Rg@Tw$KxydfN%?AOG+D{v`)>)=%0D}b!gXc zLs)hOnHB*JXh4-ed-kl0ii|<^G`!g{Se7_vH*G2>`$8>vb8QJ;z{SO7NAI+>94{$n z*V)ERWp)5?-e{YU^MJqi)Z#E9+pEmsR|t9zXZ@`RFA7IPK*xGQ%*o<8DJdz?l1H#2 zjzDGoQeKW3C3f8J>hfgxM_3*xIkY@P23TN+2AI(A3T+@{tZ<~<`DS1$DT~jHz~-<% zaYJvczx3wBAg)s&O_mU^7^_U|XKC0`T++y(BA zJCjy?GE?vBs=l#t6_59_-sF5@V0l1C`;Yr25(%C1#$-Aj$AtI9eAqA&$qMo^yctJ? zUP1up6ocZ%%0`P-a`witN3YBpxj_wOlS2=(?t+rrO#af((-mJdUb9w6rPR0M*$z zqknj~D(Kc^FMKdbe}7(SsbhT380VO~`wgg6KaCa7lzSAITiMvCsHnh&7m!L_wY1F4 z%(WQ9>?hH5Es-lNJpjKE|2wbC&!@_m-w$KJu{b-KNaN3C` zlFQ{pExWlhO)+pDBhT}eXiVmOZ!Bd-8zyb!Q^uR5eFGx4`e&aN!r^;x@NhU&2xNls zbr8pCRBBCaZC7U}Wkb#eSk+j{;J%)oHZn8>YE)Ej5|PLv>L});A%XeVu&S$xut+eT z1y7M8)_MkZrNMWvk_3dP{d4Ep5(p-~e!N*Hh+8{-Aa`|tISj0A!s>C~xv3HlYZ#3%akwDQoENrU6nDEGpJ_;_1) z4Gj%QP__Owec_rIG#YJNi9(?q933lt=yjO8lP3ZL0|C{_H}PwP%u)laMieAp|C%c^ z47OAr4%E%3nD0%sVVNI1u+-Dj>%j0IkB{F`Z<9!h$2ffj_^XBa`RHYg@B8JgV6m*s z%oye6ZFVOE0;VS>z5=EaZVOdaRbV^+rm|BgoNwNki@>onbD`pX5>m&zB;Zal7z~iN zC;j|nj2W07_$AWJb3mvHM;9p+G2$5mzzR0zxTIIjSDxoqS68>Rw7_)W)xobs%g)xb zH;76+y_Gj6V5uCu7y_NSAmDLiG8x=I5J0dOeBq~eS#iz3rn{$ydH?>Uwbun!D6i9} zZ*c}`_|()?7lDcoDfop%q6>KnjS_Zt(hfQn6%%Xc&avzr>m{l%O@B9>)dfKuYIbj3kwT{c5$JQER&K1akp1ERSzHj z!`%G2P$&$sj1)Ki5EIjDcZ?K#VWMrfemh*q+`XTPaBxa^fpdi(deKk96n|#vzuaN| zmyeE%oz|s?82(rYjA4g_*`M4hOe`z{*hdi2GZ$Y?Rs>?{=D4YN?7Y1?3Fyph!bXW* z9DGx!?NtJcsvHyCNYp_>X2tac_xSi_P&&CBPOX2$`SU-EwUnA*J58o*p-@G1`XBus z7W>rHK&6BoBT$}!oUynqRocu0P9`Tuz0_Ou^5v{t&PcFB>U4BQlvWSx+C)Q60;zFi zWP~VE9{=Oo*B24z8Gvg5c7kO6Mkf2J&=Ux0SRA6x?s1N~0;yvN1kU>U6Dtyug_)Ta z#uUdJj9ol(j74h_fT2w``(W778Jt5=%|aZoYQrXDsjav!Zm&KNy^zw!Uy c7i-JHI`_!Pt~j#`5Io57W8SWHn2?lz11NvZSpWb4 literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..abae27850266be2e2a831e075006088b6d5285e6 GIT binary patch literal 14857 zcmeHubyQXR+vi3^6bTiO6vP4qBn4?iO1eQ>r9-;&5+Wet73l_Pkh*k>grrF4pS(SJV|GQ!(@xm5C@Eg$oVKiH2=p)E&`#eUXi9F7{yEMd0wY!2Hixt`?*5_ShJoZet<@5qk7N5*x0m zU%Va@xlMyZ1(&C=jvx~6ys;4>VIn>Rxsg19fqcA%MSytu(m5f>J(2(a*8f{0kn!y; z9>rJP+OKvvLR9!-DGH0Y$ZZbtj*aPNYLESwsXoZ3#0iYS&v|*4xox(#wj(7L+1!4&DTaTRupijIzsJvIW#*f~KyE6CvI%1T#hX=xi9n^&)1?X8Xpo*mn!5@OiE z%+n=`hfw(wP{;GTI;g7;MzX$sq*Ik)cSnd5)`r0o9v<0vKJLel`;0jp*2b#~3l$ed zyp*4`;Cp= z47YGTVxzBALsEZH{B};5<>?K^%e7nV89w&$Nk$4|39ie}XmX?;WvGIIo z-1O8GKD9vXr}HB+l9YmCd1;CK#*IS9HJw!U3v3g2T@cQl2_J$F6kKm6zxfq^jT7`d zHA--I*qnJzk$z#Of!7J)42z25c3PjP=zsh6?b42m<{XB8qX0ZVySyCp`Ezf&T!O{` z6^|WU@tRINEidm;N{B#~nT3nX;pS``Oo?j(J;o7ip#03mp7gGGeqLHyqy6=~Gu^S&1Ax@|PBdF=7xP|1*?f7x)FOiXf0 z%36(qmY&|kr~6VcV*~3l+rLK~$13dQ1P=5Y5^(aBvOKVBUpBlgc->s`W=v!8Hj+bs z=loz|V&dnoPD~WKoYZIPIgCGwrCb*t zJu#s&Uprc>4LPViBP4tf*&G%Y?RS;lN&RoxbSYWcuNfKgJQATd2?;}Yrh>lN%uei~ zJm;c}=$z1Jb98X%?CiuRAn;sB^3P&Ucx=`$#a#W}GOXg_B^;>~HRX4WLBukcL%x1~ zi*`JzN-cT$`9Jko)ZbE3QNafV4*~_xUf8XgdmPuE4Lu^u5hup-Tpq0<2LR59W8QKR6Uen4%_g1^^m7+Xc2lHP_`d|22RyOorQBq1uwDiq4oASEk zHA?=9B?gBz)qA8lQp8x6zf#0Q8vF@rk5(&RzI^%YEw1xc2UnQv#Y6?GBCe)wV1N$iV#kBcejFqf+NaYEWJl2_KEfB@j=i{rmPx4(z+7`wR$xa?fmp~ zk@d9L?hlMl2SYI(Iog%|;m*j5%+!e?+H_e8PV)>877-JZOs8XOtSzAn{S`yFUWb{< z6R@AYH8eE8Cu28gCGzz2?CGH0+}hf4ad&fblaP>bbK|2PqQAq$^n1nk$~PN}b^GP- z<+gzjyvK_AI}DLtIby6X0dAG8-J7SUbIm$BM116QY;3=VhPZin4p&2_FflPp-;9}t z#dIL@qV+{`_{NH5L%4dFS#(Zjze9N)S4ThG;$6~+iH@##)fl+@djw@9iKsSq7_Nvt zt)EZu^jPgocl_|-L+wF_r4y5iAWT;A?hh4J)w%ik0ekJmuMh25?J^Fm z_4T#aQW77+CDjXXBS?-U6PD$-^GCwk1Ic+|vgtj>y*vJw5#*r}b|cx%n`u3ZpJg!+ zFa81_?*u;Q!e7@BB%gih;rfM}uCDZUMvy&m zCakUwL*+ZHtgL1YmnNrQQ`PP#CtJfS&tV`t40~t)78jZCoNIpbM)Cm$GRmw%#K#qb zfz(wfTSt#}jhl<*4UDi6AjqyPthGH$6GL_csS_cl?$%s3WS1X>3PnKqMBe=e z>@W_86?MNUM5re<(7O~SOwUX#rCu+i#&s4J`+!l`&1%<}=8cUIw=;g|_9f|o-Xx>u zO|mMgr(pqD{H+bwW?-RRNRV+io*5TXcV;XcM0!CPg)dp?M}FaloAJ? zX)`p7e+7#uMNwafiHM}*`O2%RMhj<_`m%B6Fc9^tSg0@BWmc0ZKG=;F_9{V8KOW&D zb*6Ot54w7k-cxBm#=WZa2nU`S#3m4OhM|QTk(!z+iH$J>KRTb@r#|0&s1VB5$&y;_ zaKUT;Yp}6OE z`Yfb#?NRJd4acipsd*h3A3SJkZk{h5&{5)fptlEDk8`IUo0?Ix#C*KUS@>Ogf`<#t zHoIZ_yTHI3!oF|ca$1h_adL)xW0QXfS%AJQD=S;;aUvap#tR7vs9I0Yw&7-KPd;8= z%_L)omDVr@l$z3s&R63?Cv1YcwBcg2i1_%gJ&V5Aay5#{Vd(Sn@(e6YO-+w>e}_pK znk98$U9fgqZkbE8$El zU+<}tQ&}>`?@#r5*Q6o>uMoCR(>x;z##A{+1I11oQy7?+pttp_*O(1_g(fO+x<3sC zRPi3^ZL)Fn*b)$gK(C#F=KDCDfQ)VCdjKK7%eMSI;@i|XR9+V}icMA2)g%A@`|pGV zX5XLZd;AuUojec8Z31{b_x`+8ERR;h-byWS)GnzFJdyu6i4>9Ww(B#kii z9{=e8S4ke_6lndAJ-}tPe4gmXr4|w5wTW`h$3N}#4Gr^4-6GU;G1QrkMrwE(3#Ecd?6K5(B;GhkJ#3 zSlth|hl-5mXJ&Az1^8~?uIuShx~w2CpE$6*w=zUXK~ZTv-Q?;VPwjck`rtv@*W&W> za>s-KY!3P?x_$D2v=8szWmsF7K5~1RXJfm!h7vUSb+65vir+;#icPd$L|#i@{|i7_ zJG*UI7-l&*xVV6mmZ5*I&N@?h9&flG?e-OySfXlRk>lau401y`97qdS6%&&jDWQkj z&B1^W+x{LzrV#$bJTf^I^2<@6%K8H7Cs{1$m?XL z9dihUVXb;kN87BZs3<5XDEB6LoqNH*H+jokGbLVLUbEuiihds8B0O1BgjW;ki~6az zktKB7P_5uBh5``FZM&p#w3wbyQCX=;!{Bu(AKoD!l#|0;v;V3+k~JLbN*U+*qpkF>FG5sbk9xx2t7F4^CYSw&&ZxIOYK?w6dKC#@Zo3t3w6=^ zJ~~u?{q;?8T}NkXIS&vcJ1-r2#F08?18l6W;!F%=Kebl{W}6&AvgM`fZHvgDqe~Fa z`~7(&b-qpPm8_Th0XYKqVY~|D>mQSz*tj@iTr5OUQ4~HCl{WRwkAaxjc55>$+Hw zB8a-*!%iQd-uU6#!GzTM;(=)!^ zNay2aT?TS&Q();{?wcFe(MXJ3r(J;%_{G1X@zxWmGeDsLIvbEZprySkD2TO%d>l=c zqLquMWH)?o@K0W(>^-r3G`ae8KNV=1JoR2%Wi1@t3KF7G0Ta3h|h{yP= z+Z?{fU~pR(pnXsh++Pd?VbTQ zcfZh3=i>%~z^1=$LFb@*$i)RUFDW`i;80>w+E@$P;tSYov@8sh7oB^R6 z?X9W7Nyu<7ZOleYRj#$4l~*#vT{dmzb)w5ERE*`mdwkN-kg2q_xk-fQD;5zNy0zUS z7ei&n5zRI3^z8D+x2w#QPplH*R|Se;e))!GsK7*zewU(Dbqh`Vvf#q;Z}Ds|bf-wB zUYL1?=ADA48zEP3GDZq|*39;hV`KZ*p8U2<5expcXt!3i^`USBF_B+K%$NGqbmM?+ zjab&YPdDMU3DeT1%fxV2166tc{5g&xrv59~H_ioK(m=e!nRS75?>y43=v9C00d;bF z`*rgN`Wv_(-#^3|hr!pWu5er%zs2i7$?K4#lb2V+Nl%{%g-+4{^Cu|>Ih3xA?iBG_ z0B6kb2OuVZ7~#v%uA8Y`r5EkIaABXoGLC0Ak`)rFgP$W8`=qO&OW^VRZ&TaxE*tp>CvM{e0*+DeU#FqH}>}Swzr9dXXfXZ zcMkXW_dkC82!H3}U1b*|y~dPa|LWuvD=RDKjO0I67Q{cGK9hfhFCVn9umEh+?46)vawHKL5T+4R9F{k0Ei+p-?+PmLKW)wWvL8&)u!OFTc|xdvBdHjx-Txi(aNCG z;Jo!Y=tA!uc6QtKiF*57JUl!QNphLCkV2>#F`YrxF!U=!Ma3GqWwtsiy_rh(Jhw

PHSCT}tF$bCe75D40j2a0D-Gp|lrt4mupTC}4I(|4!eP;a7$$GwOu*UtUGeJ;y z$O*|YKHOe}q3-D)0l*Doj5(seef#$Q zJ#J6)8|dD2bQ#WDbFm~>B_8}~?*FyYVOD9K(lW$d;<%P|Un%{oR{5{KK9zE589xmi z0UDZ@Lxl#qo+riK#{V$d+1LQ!toxNJb+ogbuE&W4(*Qy^*FcQjI%3XdT2L(#v&kGx zS!F+E{#W&;e5`Z8FbjZT006@i=9^y60Vv9*EAe36#ss)B-*|x~>JLYJo{>C}Yp=tI z@G34pEe3I6CJD{ToC{zmmk09la&mI)EmT!iU%pIwY>QdOOsRZz&VtR8NC$)9e_PSrgEwL4LgN#l6NP-?9@U5jY&zDg`~CL|~5 z|DpoEicLVwK+o-Y>Sk`9WAL;soC%aaQ7l3Xuj<5!q*s}lv|C%W?q}7`#vD9vhJ@=e z>&~ASOL3aHUf^8htp0+ V>n8ZwBt`>){wnmc>pjn$m=D#>y4M5A;}IimSpcmJ_O z=?uv$CQtm7tp$M@8Qn!@{5?i|9)9?&ubJG zQCn*(y5E%B8iDkk?+xS>4KKO5`zrcv6#Ua3fbzsudl1}h%z;;)hDF|qFOY1MKzTb9 zY6BJ)7P^GyaGo7>d49&vc#m!^-KJfuBSwd7qpoZ8ahDys@PVF5>$v_d=Cm00_4jAE zCE~87?N?iSASfu9+EcfhId#_sL^B!x(9j-`pUyTyJZENR_74yH)t46+htlP!MZGaT zDg2vu1q%1sQg&W{MHckxA3uI@c4A6^aw7;m2!`DaWWSCWu3FdM(6^wKnJ3@~r9AwP z+P8)-2s_j&D?ac$etrQ&4_fa)TU$1CRy0p|Vx^a-10AsH5QWMqaNkI)VDrVuCx z-UgR8S{?3g<0U;cH6{G)83<^UBqWX?J%2twzJLEdIH<0sHUM1>bj*iw*FXOCP#~FP zB@{FI{l98q{ClIef%=&GXgvNB#y&VweKrtr=haK#ml6K4%NsyPWqQ~hi7!#n3 zXqyq-;+bitpYVE|F?Gug$|aZ)K+@Bsqo&r@))p7(eL*}QSl-##*r1`IfxmA0*ex&r zQxB>kAtGAXIfTLc-qiGFv}^<9c2?Gn)Cip{v)_$DB+R-s=DF=&{Lmwdgd#X6!N2L%K;`$jxAL-(-VAARtaixOh&(f@n9k=J-RIY3B@(}pt5tEKWdck3eO_DpE9IKkHRpt@*^geWM0x$ZUueCtBs z!z^nYq0UYRIy!QnK273C>4q!3*Y>l%{v~)Ij;Q>Ly9|ts(&;_)aTj1k;OMBEBxuEO zTIe{>hE?>x8LJ2zv)Rop{6!~XzEe9n?q>U)Ku~_2?Y}oJ(Dw#u{Szi9ZfBk0mAaPN zDEySy8pJVHHZA`t!*9nECzkS6Pz2$)NK$gaybeDw6um^t@`?2*9OPllG4;RsU=lYsV{$1J?E3?y>b4p-PImBLH9{NlJ{b-`|f;r z@Ky!T7W#wqs^$4%)s-G671mQC8>O{T&1+S6ODxqmuIq#0FmK2%MEKn=*GmkPPdFK{ zu77(VbeX0z*7;;-;Pu+5QfGbrBp~ElZ2IpClU6>%z~sx}LnkS-xqqLJFV)j?Yu$(7 z&~pn9%hR)0as2`O6tA?^zs0{NKAmgECU{8RSOcx@FWtHA#r1V7E-tP%Z!&g6=@?F# zdJ!b?1`*%Ni|0QY`?5uAJ?19+bD@Q+KuOPl#>6}DHEw;6OAy=~A`*t^OGvM`w^}7s z4MAqy))OwZ?k5Ox1SQ)CZHg@Edr*?3gR|@9 zF5mVWI;0oze{JWAcVl`nx$RH-_pQ0P!?w<1_=*|q1uaigLPqA~G&B9X@0B@g<kIiEg$69KT4P`B^4Gvb9SDa45rXVc_w&cuc(6sabjO!mE>GzG1BLLM1-d= zN9WXsZj=B?4Af6pJ;r1d8jGFN_t z1A~ZfCM?FaPcPDp+f4Ryfz51DmjK1JYl&J9ko@a7REk@>UwP0km_^tqAa!YI=9*6^ z9`a_Sq$4T%Z13ze;a5~I8ZQh6{8D7~sA3`SIJ-|fKR?e9_q_UFdFEbXtOl*0Rr{(A z=3+p8>G|fo?txoegCJyi=qv#8lkSmL*~e1r6jIK1J}Nem?kC1JbWVPNXDoqMgui_L zJiz-Uq8{b%|I@T1AT|JKimp@Z~csxHtu|# zrlvwWPMic;#Lv%T=Z8om{i@Sef}^{5h?@S-d^! zwMT3MbJ?KKV>BNyG@Jw5Xjy}ToIGn_*&-@kK}A`4M#Xk(P9~M`u?fiZeP`UbmwbK& z%nW>)<4wE`6%8y@8A-_{ut9RqyXE}6(X>RcQE!HVYwdC2()+=wSJ%lhjgud&wEWGz zRx4;UE4h_l<<_4667Y}heogE zXC^046&BXBvZ_b{29%WjFk3d9CI(YQMz*-JLV4>}H8_#v*%!{eJf1hBPXH3-^UU#* zOBtjZk87=rs1<+Fu54|UoV`FV4Z4@_c`yCX%@k?dkVo`sI^53or$b~8WTfY}F&`!Y z>M#L{2#dS2vhtd6b|18B)kP76Jm}(<`^sR!+W zm?-rvT2Xss7!!@o_UGoM$q?%3%uTjOy?F70s7mI!+how@=H|=VlU!2nLQ@5eXd7$! zJf-vx+7%pW8uh$Q&9+mqeb4z&Nn$5QL7W!tAxttT4xn}h73+TWjMWOJYsf(p_|HwR~Z~?Uj?D!}F#U}SqlbVK3ef{HA$68x^`(nmHAj%)X zyO8e9fCVh*IvG?MQlgT3eC)Pyy!EP!YO5?m!pW5(`tn{EoyZQ?_7v)Db#EaBc32VP z{>hrBAV`&1pN7~UckQ*dwk~0J%+AibZ~magcXTL-2oLu--e~Syd)M}PFflAoYY(c;sIptZ%l{%=)Vl#ziU4k^voY$yUMAXt;+47(BI$Rpwkm{ zDV$m1yDa3Mn~PF0oB=0?0nJqG1YpDN$8tlhD$QBz&(%~`{%mn~nA%|8)y8HMY$g_6?0GfLN_Gpk#e&%gy5a~^Q>UoUGF&QjPqH?WlNI*& zm9?eHzRZkLpIWUR%y%@*wvjOFQ~{i88xS}<-r5<>Mu%5Ga4Eta8yyAtAhy@e{T13I ze&mLB*$vvl`-*o#*)ER%&So>yGFIuRrn}j4H-vA0!Y8d61_xR)9zK3a&dKpsmlW*Q zHqbap*_)2bl7&G8!t=Y<*f_N_&^XsIH@Sxj+P1S^jdga~NN)LY@MC}i_f;#=%LDZ% zwBl8J@%x)nTkE~_V2|wVT=HQdH67A2H}9MK^-ILTp>NK{)@&%&-`_bhj;X}_IyQc9 zO?j*2OODeSPv(a=+9HB#PdfOHcQyZxerRLU-XIn6-S6sSE-nj=zPZi>63vo#KG+=% z)Sg6KZ;p@WO3YaaY7fL>dG0^bVj*DGy-C8{U&!6Yo4+;R*q?icwP3g73JoQt2OA?J zQ9!^Q{xLp24?cb_t}!;EZvc3}?kmELd2HseI?BDOs^(L;i-UK<{ba zquREJ#u3XP--~T>;xjXs+eAI(Vy~s9s&s|N9EDnpov^7)+b1@O0Lh^xoNu{4xumgE zn;^jSRV$NRr}oiTos%D-S819))Z`z_?4Ili zItR41-Z1J^#)A~c^XDDZ!)?r0Ng}eC)eFBYDfL|?`A(v%e2XM3*6k7Vn%2vtS=uVk zv(rqabdC7-k&((ns=>HQhk%6f_NmPDc*VZVc>IUuLrs5TJGP86R^oqhJF6StNVwjN z?yq^CX$}m7mIFnw(((zqf`O*h#sc0?ax&y%?)ub!3=2~^K{4LFdl%&8&CD?O{G-Cqw>niv*eCqvw_apPP^8wu3G!$ z@Y(Hjcnhe4mo8oEFX9GCF5EHgFjun_P03*dTt~Ga%m8(fKZe(K?md`M2DTB*x(-mO z`sdz@O#@rcvXo8f24gsDiP?onE$+87a=p9$h4@kN?J9l?? z9V)k+>-w_Q9+L;}cII{ko*xi3k#2yR3lcQ?i+5%ozm1BD0(UZ+m!!J0qr+mTFm?8* z+IsqJyLO+0?T(2Gfb;=ft7J?$6zudY+J%QCAuYiN`07Eoq`!bk2tuG8 z4pLI*hMnUF-Z=8{W5swq$<3RF`udTwZ~o*$9H}o}6~aG1S?WY;#lb}{7giJog?!^v z2eAFXv5n(({16nx>$Yztd47OH1^Q8TwwZ)P%XqbGSm)(Bp2#ΞP#)$36rCPk=I4 zS_-^m12Gyz^qj-K`xt}FbG&*Kz$ZyaB0$s0`f8>`OG_&VYOH~Q0jI2iqPjXMF|o>t zVsVzDhK9?=lxUw0)ePt|5Y#Eu%9gUQv9(Qb-$mOTii@{pg0kk0KtTdzVOCxd>p~jb zgx1)Ar>gQuFachjn40>f5ysWT2MC^@gyoH?ash!(2UiSun%mgaZY-@frB}kx@E#x= zc+O$?ifCJ#_~uQ4on9qL?94F`@4@JisUL~VnsWsq!NcPeQjPjWOyenqg@xeCi^Mr{ zLTTRI-d3KZrY^NajTgMe#zKiVn7sjwtoH?QLyA8cJP6#J8iwDdekmR@+_Pj}EF8Pr?1^6{Y6T>B$Wd&kHjG8ph2 zMln&d50WJx$7Fk5p8)&eRA)lU*L>Ryto6Qe06{Q5_u}^!ULAwYIgXxDxyXO$dCc^0G3} z3IgdLKD2=~<=}vxC*`Bol(~SEVe05R=j7zO`1b!n;@5KuAOZv&tQ;NpHa4Ojo2jCy zTVIT24H+rv!qU<(T$05|$&kPZ@L3W@^$AdfR#&MgDXUET9)b$##0OCpt4U!11Ut(E z1L`Vr_eQ?tGQ`yDsTUc_z|duRe9FhLm$fREv3O9-_r7d;wA@xBcbC9(zhSb$|J=WG zcJ7Gd9r$M{&T-tNKW7N?10j^$n}^GJygWP}yF*5tAGoZiz5_ko4)@PpRZ z*VSv>$_>zAtlnZXHG}EujV224>Db1sikjTAX+K%~>JQ6>piobbEo5%O@55P)C168?uF{$TcDnFn?KxI{{4F_>6t&oKO!O`uF^oR!tWa7Tf>=m_x8BRgMmkK2F*`Qtb^AF>gNH3 zeu#v@7Bq+M<$2hX@Hu#Xd0BIL69$&iLVE~E1VjPB+(MV(0!x%YoUEFs`LdzU?@R~p zYztHNOoLVGPi_Ehx@U)K!D2;_4}@GTjc|~?kHM+H@=Xt2;A;>bIlQhs*K_COu)4V_D1#s~M zUdOM$%2x(qO>&tJ=0mTTn9xftO0>>&_GT~6efM1Z7Ik+_v2nNXW$-rT+CG7?2O&== zzEJ`PtS9%4UmY%YZ3cYAixu?rfQ+wMe~$A;BZ=$w0_)3fzW1fVMMOo3iHS35bntdW z6rT@6hmO$oIDm!zTtb3Wz-_;)slojWKwxr>6%^9aaeWRG+-uy|6EFppHY{ z8`NhO z##AE&gCS)FR(GdAwFlaCgdR1*sF|&r3vm%ah<8C{2Qvg3pebk&u;Iud)!@DHfr>vm zJREEah+{z_lmsdmKsBJZCLGb{H9&g^imCW}QdV%E#i`vnvdIB|LH}$#OX_y&<*x7J zERO0`R9HIZvp<{)Vs_bG0g5~d8wDng6!%hF7P#vNM@NvA0Y21F>9|(nUIO|PsDiMl z8~gem)6p3N+J?|oYtu4h$H!_F0P6@+{J3PY-0A>eQG@my;6aOUJV}N5i)wwaRNKDS z8*}szoI%PHzkuNXiTgFUnY1gE7jerx$H&K^A+MrpYkxgmO}L1AqADO57l0@n+Cm2> z1Z;JlxYg{u8aAfA-Cdm$^Qc8@2;HN=UjZUn3L2G4wkiZ3p&uT;yM7l?Fpt9uDGAA5 zktRJea{|9>Nk~YDzkf3*$8n!O>uY#`EdX`w=-_~xm)CKr*90a2c37=(3Z(U6#kd^I zMnGwTPr%vW-V6^%FTkZ}@nEUgKHn&bkh2RJcQ;zNj-c!J^zR}I<*s?55_rhdTqhT<$OY zA=6>M^bh1G0R4Z}sn)pM4L2S@J_?9BBpBK{I+%booZ^X!h&Z{rt_>DUqm%k9FNGnw z$i~hNrUaxEA(K)LXMm^%bc4s1<8Bb}gdYa!%FW$fXehQXUJOAd1LE9bA1F;+-Y%iO zfy6L6HU?-9lq@cb;Te#4C@9P{h{?$h!A68!ud1^0_R5egnEd@GG?=D_F(eqvZ&D|V zsPIJUNB$zDmAy;@IuS!pepS^mL^NP6H&q&rfr4u05*ragiG zV8l=UE@a`ikd>7c{8W$zRKIFDZ!ai7J%rw02Gem82Eo~ni+@E#42(4#A0NAv*nngJ z0tUK#q$L#^*tak5GPL7=%AH0-YFSUn|JR{lUQsA(8lNUJwunLhYver zxs~|qK#Bz{d$O#1He5IpWfb4NyptBP1HBuaK!VsUtZa+XvI3|PJ8sNg`miDNc5*qa z`cwX*bK(W$@hLJA4e6OQOSV_XDnVq2tI6A8(AL&oCL<4Saw$=im%r_6Wow&VP?euA zkG|)FNlH3IrftnVg-PR|d`6CH{m;R8i;SA8wf)vUCJWD{ZU_#=th2tWYB{r;G7&9}MQ6Qsm0yV)eqZ(B1+ueX&U*JdQ=HyU9OanwfcupZ$>@|0? z&gH-1rX09O2<{xx-984bxqwdqc41JEl#q~6Far4=uw&iVKY;8cxL$_|{}z6S)&@bi z29wJL1S@#`bG|Eqn%hQSTDslR^GK8CI}o2c zz6(9+A@-$Hv$Kz(DMPujUm4s0(gQnc3${~&%TlIlt_J#63{d7FMG6KcA9M^?5Vi(t zPgPP53~Kox3@a7Rc+fz7`Y}Gf1{6Pq6tTX+TOdXrfZE3RqVW63$Puieme$s2*4JWi zpVsc`n5K@-X|v!dD#baNNLqTuLk<{bBV%bIk!~fA7fTney|82iM0XcD=jpE{%V>+QxNvLYsx3G(BI7HR7O%!>EJxT-W4u-B@w|UCkA^mP)vhP1}su zVP#>Vg_g&w9sSZDLPKr#e?VydWlZwFbkcsFVPBT+GF5c{eaB435|y#E`4 C3dFzw literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..ad53d2ce2d3c31fc6a57bc89f3858f479e4d3189 GIT binary patch literal 23912 zcmeFZbyVC<7&S%&NFYEK4IUu4J3$jzBxrDgTW}Aq3GReo!C~>m-Q5BNcPF^JYlpX| z{l50JJ^j0W(njLTHv|Oq_Q!uux~+2!5fI)Ye1d*da!%TvH~xUF+K#wXE_I;w(G(>R_S9QV z9@-J!wIIXrse{0c8A^GMSSW}EkxO*&cOE}dlL;9DfzaZW4m}AB4BT8XXvG3A$#D50kVs7Z|Nq7Rx5q#MhxbEg zJip5^sn;!;JqkHFT)3j5;{0Gyr_uE+*vq_=?lLJUsivl;pnw5EOiTRCHh^`WDq)7I8T#bq!V zhmVW9#ML?M)YaE}MaPkole;oAGBTzwsFMw}A|oR^?M(DfOmOk=Tpq9W{rj(^r>ARF zTaAy6ot~cyx?h$I;Rgi<(qZ}iZbe2yLVEuEZ#?@(p{<j2N}a^zRm%B2^+s7he^^4k@E6d zhyF37WG)-)59^EQS)*a_tvXK}WMqKKqTlbZIoE2gH47X)$dMTxVo^z})jRGjG`cl5 zHWn5ZqWxDaEiBwUJmivjaj>yD%*JTa*)bu*$Ui0CU64lm3w;FZ zgO9i0NEC^KOB!!uq*wk^S&+=BO4IHAZ>%%>u#10?Sn8Z{K0B)}+CPQ}CLbB=$T^>csL%ail>P+&Q@dnAv;LEpMtYSyq9GJWNXdE-h{Btb5Km$^L|NC z{QbaT9^TzELnk+8&q<_ATC{Cb#*xFq$p6(oncWB;e~rs~4v(0;O8&KqNBu5MN9fI) zmcOqm%?(DWB(hqP`3{f0506fV!ARBn)UjmbZ=nl9TFZM z9v0R=Fo2|MXt=t$X<}@AQL5eGa$Sf4BA1 z)}YbE+WPPG23wY@dz36yPwCa7wG!v+4c?hg0}JdLJY@Y@8J$rYTT)F?cp~keIhi%G zt8K!xwRG~QQ04GK2Y=$l8RDArsp0fvRFylhhw~@0-Gyvy|MEtj#`Po_%)Vb%=vU#p zH9YPhGc}sIip%Wp4r~9-IyBxnU;l<3Da<^p_`McsU!1|v+rZAi;I=@>;d$3`C+RWP zY%SKwhQW@Cbkb&>fp|%I`AN4+{*4i6R zVYPd0b7wHi@l?k`7y8qjH?^AgJAcYKJA1;;GCcn7Cj^gnPY8H8Ir;isPHP?v{fT;s z!JJb->Rg+^<6wJ2!8x9Cpp*PrlYyBjmQ*cKO>xgfF5|?rvW;DVl3C~BlYM&P$D|Ci zx#hv4=H8!i%G+Ntt;i*X%9dO6TbIedNUglklGpz0ms5;HnN?=w3h{>(>{7=UBV|!y zs`k@JjYt0y{&Xc^o<~_{CK~r(Bo&f1-pM?ZxWFK!_EWildVp|(vG7*cv1(45f{$z| zhucdn&rU?dMXbxs3(gilH=};Pk!n|~xxc&Xez8A)A-JaaIq=)z!{5=-lkM?*%%Igq zx311k%h{?uEsczf3`GS6K{F1pfiuz&dU{GqebZer<}IxBwT@w8&rU0r?xy&iT=GHx~*7)z@9KusL^=W*CT~l z$3+_T#FgZNSI*+$=9+r1o2P+j_x*xGKF5f)WPY@Us95&bjr77DlYzvfNjL4gtsApu zlx%}E-rR5N+GP+ofR_bpaqXB6`X!s zt<~q0IXQMy8T~V%%uM@$^;sEjCx1P08#{({jP3K}1Opqb#B*tFwC$aY%9$8*G2%tF zB}+3DRn@uLXVwQ^_wqTJX_?QxbyNy5v2R;`JQENQ5HggLk&&UG$Z)b`q@j_JmA!Ya z4GhGNRtnx2lyBr!QbJhzq}19$5+T<;lhCKA0Of2lIUDkdl@IXrmlE$npC0@0j%&d6 z*Iih|RQ(q{joC{^R!g0xK+EpBeq_Od*{10#)@rRfE|>QvCh2D@l}Oj)9DHF6j1S>5 zHyrPomC&16v_C@GUa57^%{}1GAy=jf{vE(^VgXk_*!Bum^7L$=z@y^nt$>Zn(LuIverVfJKsyKC)lEi)#Iag^ai#nc4_B=~W;>Zv&Ln zn*>GK@?MeuS$o}Rr4#hL)2_0*X@@+>EPZuP=fDs1>=!qoI#bM=6lXcgSRY30IxY?# zg;sant*#Z5uvB89YaQ{M3s(+SNMxkVxAq$~wGZORszqroCl~uSFKIjIuDw*0luj>? z4AK=}LuhepiRz!h8sS2@d65%;-#O05U7(;2572&D8JS+Kyh%y%IJ7_0ju}EN)U;EN z7n5n8RxR`-HZghQTyNY1l~{DeI8!~d*1Co#B2 zK1+MY?v=jL!=I-T4mFmIQuLE!M=%t7=cT;i9b;vrq!K1tD2q##^-O}FaI||r*0^6h z$LbBEPPq3|guuWLaZHWDsp#NfOsg?dqRK6`H-6Ws*>Vt{wXF41RLEY!`udt)w{k$8N`f&)5VE7~Opa4d zrH7q*`1QYc7sH-G?RdXKH`N-8Jj?382ul`Oxks+WBGSs-!NU&d;HE!StafD(t;G0Yw9Q?VM1eXaq zm_m9PFrJSW8p`sG3KQ$u!#wFpcZM2i)QtunljbE2VsTDRwd>h-B!#$Vv=S7GLnT1G z=pN~C8=4T+)&4H1R*&;6NC)q4B4b$K=$WLP zpcA9e=2v^>dbJhV)*Vqxd<=GCqPMs`PDQQ@p27kM@p7nLRkb=l>BIiQhudz44~G-3 zjR=gLa+m~G?nm2lYPTpAo2y4~*39#w`DwXQ|KYJ-RmW~P38iv1(NN_lys&d?ugJ2K+C1RsY_62`pHj8TX%!o z#*i=BxMqSnYSXp&4S#Vx+1?+wp?>2_rt+|?Jnah~Y_7TEpgOIMDUXyfCk*Th7uztB zBzcH8S}}rIA78V0n$o6-8c*gLf4<+Xg?gXnh>lxw^gXC$Dy(jZ=sYlcAz;2#dRTFv z>g@d59Zu}+E-g!j;@O|T`B1{i#AMa(kH+9k4p&f6(6xAlhW2O8L~A3yyPP6q1KW(# zsdF?3d42wFdAPwd+^I^mism8uyvxf-%D!lAj9zFeXFVn{PbET9o0*T9E_#w%QYA8` zszf%StZfd~bQC5KkqFkWGoQD*)s&uyi?TqZ5A)^4evee&#E zFeky6oWUuRy{WFJ-!qj+oXGN6EHp*hJMrC7?=V?=TDEqpA!GW9tfNG~lyvtL-KEy6FgZRk3QUaVCwiiF~=#0|#p{ zL2j3wTG@5{=B_cOW*QkP5n)?S&M!y0lQ-WSc_!c)ju#616|Frf{t?^@Cz>NWwh#Bt zD&DTH+o)Yx7n&1DA38FBiBs#3q^wTvm>4ci?XT_pz{pxn-Ru^B7a%8Zj>dlKV6%?6Nu6`c=DDC=}!IK_FC2)3Vj)U z(D64JjxjotX+gaYjeJ6)T%&o}x#Xi#Dvgtnt?wLSJYpM2aW0o`6nz@_AyP! zhf4<8+z>SyoPU^&_jaw?-5%H9Gfb{|`B&URVIhL?W6CZg4jN#E7H9>ZYJN3C*Lf7p9 zhf!v+4sB1p6*I$svuV5J6`2z7YtR7s8ZKTyxMlF%yUE!N5*n3Jt;UU$P9$6hC6bAx zv#XutiVCsod_=^A?3`ZdzuKunca~0{9%d$Icl--U_PL22%yanFU3sM?d?4wi7z-=j@> z_ZMz(Zp*xc_@J|IsRshat}D9o4$W&#PWx}S8O+8$M9~h7U~apm_j)7BK~E!l45#Yx zREz&~>Q|W1noVosUVY%-mc-z&MMe-28IJCeBC)iaUyZs|(#}}l2^T7(FD{q2-$;

syUf$7sngD{)8R(2aq!FE;aiB>v;+qaEMlt@)&8h=)C+gCP`6qsv+jK3 zM+h2YeqLrihs(--dicpr2v<)%N8eP#Vuh(9Hr~nxsbri5YaI(`<38+@4><*;6Gom2 z_s1wEQs<*@ulxmbSs$uBUo}&6VozFlLMoo` zL4;%DScwRCnDSqhqP?}Rv!l>)HC9rY>I*O(-4TVSh;bp1*UNEI_Gco zlSykcV!`Ce2IJxVSG9U#>xAYjiYX!E&|l>t#a7Zg<`C>Dn~l&FR}b#};_T&FmE2V& z92=y^F72VyHH@=K*z>@4Ma=tR`$;(UdLlfi(LUWBIhDtWuSiL;e%%a7r^*ps>8c{& z%N~!S9vTro+is)Efq#E|sB}y$2(JH}W=(zjG=4lHK`fzB)4fC;9Nx7}87m~$;*b1{*Z7GZA;BE^#%b-J~e}amkP&yT)W6GU$kV)n2aL1x~kP-%bu-LNqD`h zD%P?}buZPflfZ7ey$!=jEH28U^xy4Hh|sv_dQHrW8k-2Gd;Ii`<*a{(`0Bq=Pu6-8 z?+wJPyE$sXM3TZ^><$gBMLr4J+X!!?^wv4@QhQmG$pmIYi6^M&QT16*&-?~tm9bG% z0ypofETWm1-xxS@K|}(CH8l2ir(?vo$36!|i;8laERAF!D=FLF@`QfFyl{Cy@@YcQ zHKDOLgO8{kLZJs2H1Mg~R4z|+TDRSeoPTS|#zXB}ygO@fv-DCCO1qe^3F0$2`|VYV z-c*!y?cUVnu8ErA&ya1FpR=(o?S>9#=#hIwj+htNUM0JpO3-?k)tSY5Tvx8+`H)4X zNnyQe42kRB3G(fF{G^C2tLy0pF(k2Y-=LRF+P>xN4WiX5yW5zn7p~QfZhC3jCdF|` zi}#z0IyDVYuZYJ9@^{!tqoaRm&LL0aOcWn4K)|(o-QUZ+UR_bAF;#st?AOzd&rn(OEihWYqomYhd-?~_Y zNUa8b{XtiJCB>j?L}jlHYwrH9o35JKqRbc!@e1ctgBSw?oJ24-UdLG!~KxN(##Ic92W`~3MkKKtk9CS6F7 zAOe^vc-@GCZvX`+&e?YB41o^eZc3A>=;)J@vZn~~ubw<98IdhEPI-A^v^lM~*#IT= z1<$WMvo(5n5s*RxSM8p!H<>3OOpp}&RLp5LS0_XG`}bQ8Y@@D``D0HT0SbFHnAeZT zkElBH_`hmRxkqj!=OSMdwKb3X0`GI}F5pQuS+pL;NJr8d9c9gI4_QwX>i*7EQu3OY zdAPa7n;Hf;pC{NdEM(NS?PhVH``gc1q@{6hshBfT@8#(OHCMhHIa_l`cNn3Bl!z(FXH+)RSH1PPOSUp0_cUb8pEvE2svWUk4<=?-lI-EicgUDgY zsO)p$0Ltr!<#dn>VMsS;0pC$r-7#&gL7|JDlO1c2P&D5@k&SRdW4T zHBTs0_f>c0+-{j)F!b~jSmzS)ICya$tY(?37kyMalRa4q2!DE7bM-p?YlE~iwqL%9 z*c_Fg^B(N)4+{W)7cq#5zKN}>iEoxbQ+Yn5(O+BZ7ji_UjK~8NOh? ze*M;7njB7az8|^MwfK;?z{BwSSpZHY1%=b#!R+dJ3YM`2Zy2F5hgDMK(90mY@50#b zpZ_b(o_7fDS8H)$gxm{F9ydejB2VnM$9iUFT)RT?!P>h1pjVGoXP(ruHTqfc83G}p zZ6en^W3}(z-h5tE*3ecUU7#=#Oio(bVy0rNuXuiLPDx1#B-p2B`5_@v`ueaB!-2WS zk(C_c>&w#6(>CVh?!q)P!!9?VuB1VSi)oRJTQ*M04ay}!GWOX88I zivgMeHyxcB=zpWz1uicIM48}O9bN(Y&j`|p?+9=3$vAOFvp-UD=S#_81X*-)uLx64V*NKD|u)?a*!Y3~3GJp-Rs zpJDT&zdTd*W#+Xg@em*pFf%k*xmX*bB5dY4tC77aHx>;F8c*f_xHq%JO9FSr&RtsU z{GXWD?5wP0A+LKa4NAvezM$tH?Hw(+H(6gT9!#8SbNjn_2 zEiEK*%o?+#@4YkM@w+}ZHwTz5Ia-1Yq`-VqQhg;QB~F%&_FKPjgnA^E){+I$R_3Uq zn9p8E5QVRFAMI6JlHMOysvO43upyLGKzAlAoM5Wh`ba{ax0iW}IaHLCS65fl)6CeOL0kq@rt!+IE92H0(*s0*SV^%-p)=`y1>fqVJoQ zAln6E+-Fr)j90I2cha5R+-}m+(pr6xQ9*`@9pI@C6$?XpjGz+H3l-(}LCUAaywsyEG-w@M+&}}jJ*g7QY(GR`L)M`&k&V3Oaely@`!yDXR8*~{#tG2 zRx#PzZB0l1Z(e}YL81Jyf}oc&MW>6iiX7MeLp_1m5D@9GCo{wWsd=_-Rm&waj0Qr)L#(SoRb1>@#_1`Lcj z2{@lJFzD)$*^7~jc0b%=Mz+;DDBj=hms2J)!-c{>O1!XIuqZWdv+r4Uc!Do1P66k5 z?1K$SlJlI}T9tDrn~6y@8XzDnY81y}7(RZ3C{uwYs=(tK|1#cn!a7-Ju*7kkSr&&5 zQ$FhvD*KK0GCsJf?@dijQrNw(+GXFj@RX|g42^m9I1b_RIYxbB8x=><7E|97xovmf z@H#?OEVgqaF&|HA0!>jlXmIzG^ArJI^_OJTJvF?FR(OKd$HBIF3gQcRC?BIe8%)v% z02CC{ju&{`)qC`K6gWRfxgsH|4vQk>QzkCWEu;LC8uSw z2{_xBn_nb_=8gdWL}UwQ%AcExPdJo9mW)l<^ux*Y_q+N3epz$Jzf2CvwOe`Q7t|hN zgBlHQiWWY9QQ@&NZ(Rzc6c4v>Yx=gvo(=}W3Je6M!r1t(G3{}l*SSwQJ3?+=yfjgr ziOF9%dUn5XF{Ro^v?3wL)WJs}h!-CpOKMNBbj%~8nkxbAeZ;N8a3h9K3`OS?x^0Rv z+wd*#RWb#F$7+9>lLm|k_q3*Jlm3j8IGor2q7!HGL=i&!qzFb$P8JmQ%JG`?)j0Vv&(X(C{Wqs-w8ZLj#*|bY*CRUMU+nR0<0mwU1s7dBmVv4U;bf)*d#D42*sL z3ET%2+3jKV+urNs3bJ5C;BV5fz{VHvbA!uC1MO%AAXkTN8=#g}Y0Fd|mS>MMqL&bB@t(CeRy6CqfzZg%VbfudJ97Ck@moHa2$P zH|xsf(F%g3q-3^SYEf~qPL3=wAt50aRuq6nJ%((+pfKFbzLQ^Zxrdj5>~S+2eoOJu z^rG%BfA#vp_SKii38#BJsy(hx!Kv!FH}mB8`)})jzoU{(0`Z(SMiTf$gq-&00AmaE z{+Gm2fr3h(URK?~@GBYmb(hQHD8~7m7C7z$tnv$}J?JesU~pl;Pfr$U?MxPXeK`as z3iT%UtKX~!+E&LWC*>6tnvJfUs42OaF1S${%zdn+lR3A#?S3Aqg`+fd(dK4HECR#S zIe;fFFDo+}ho2qq{QC9Fs6QcQaVGOq)YzC3nH^-EXA;f>@%lCWX>2Sf!lY}APTSV!l#A@nr$%|zZhT~qbl&F--e;kv*9 z{F`=W7OfH;j82pceYp5bc{UK9i;9YTkkF2gj$Xca;e5I&pCJ~epr{DKy3N5tV}^JH zN^)f2J(|N=m`G7Ys@8?9{io`MZ3S9KWE|~Vt@G|N?DQyzG7bn_R#sL{2Mbzt_FFTR z=21~mPb}tYAm72n6&D|XG3Mvz_Z^UcKFHoNt{mj|Yc-b79E%9Te^&GgmHs2Fm2I(j zJWz08amZab;=s^jW@ct;nwOH2l9iPe+$&zc6!SlzWU2Pgzy)Sj2)ES-88i}-h7mO7 zkN+pE=*_n z-RDHVOPFp*tsEK-gfsnLdzyp;$YB@1`8haC+?EfY3FrtW@jgA*NqzI6yhnKF+YcJ1 z?eMcg@b&fOuwIxO8XB6PpYQ&!fW!m5?Cz%K;85DMhn~ugD16yO?f#Ld5l6av-zaBM zNu zhc?PRP!qyAR?#VnIB4f9DGEig4HN|V+uOIg14-|#!`xh5(LsJVQSb%8RbVcve+z`B zR8s=eFL_-$i_}&#lk4NUf7g_1iPq8jpT81Q>ZI*>B`7GkFgbalrBUy&gNucw80iwY zsS`g>Vb@{f(vejrm)9{rKs-iiUU%V8S|=<`76_;EqjCN};5O#ykC`oki}R+i$wEi% ztN|psxzOcAZYJVg%=`W$^v6?R1bg{m^8clq3g(0rBJO)ibjVEWw3<1Jyu=`zIBFtjr?^U&ovc`7Z1|-7#knjmIcnupPdNI!kgXIW@Rfb6v8-5wNS7| zebGR70LL}lY*YPbxL&97n1ZuFxsX$*LmT{*_L2X*LnT-qY5U#5?VmYvi^&<5ko}aB z{68iQ=G8jOdNp~ZN*38bP<|DIbAisWlz;L6e=kvswHZ-aZf1H+4-x197lux2k`MVAAStk zg=J(|&(+wTtjln6az1T zgNTTDwlyj)CACuw81da$Car5)A4ue(g~j2HASZ$Xt&)h}!yb^qUeW6xSvD=B+_TWP z-Nwema=?>aj6VL-u~u+%vG|}#t{@FIhj53V(jXUHt&zadacELJ03E&|%gM}-=NaUb z@fKbHK-%-Y+p#2E;>9;`6~@^iF>11H;=#8dr2!eo<+jevl~UL-MFi=Jx)HHlJ$5=d zWcUFK^lynX(St2ce7AF6S92C|J|ua1Zw2*lnvwDH(2YP?KD zF`GkDO2nfcrd9b=U5&c6_2Z@q&|7T5e)*o9yx!B!9@Xt9z!DCb6nLZ3C`ylWd+0M% z8o1fK@+u_|0imd<#{Ok|yz|jhN5{`=%*X(Pvpd7+ax$OOtHNBqjluoUiEVD3Yc0jZ zZ0rr|v1p&%vs1y<0ZCiM`1fy#KBCz~N1?`Vetgwp^0yxzmxerOTB-6mV(SO-k4gs!~2Yl z$GsWx2^fGi)Yb3zYLw_KCySEdaJuiAnVG?8qM1nql2$q#_XSe5klJvoT(Nq=&dDc2 zt@@DR?(Ely#YIIazkl;sOtDZ?=Nzv@QOn&ghcN%`@ApFy%XAmhAU(Uxkjj_`HXZ~N z6c%&bQ_f(G>31)P!vY0CM7zAXI6!{>yZ~5?cc)5$MR{mq!t=f0*WPHr1#p|RgnY&H z;80;K6bYRSSATEW^-4{hoDEZ}{*TuvKKZd4g*^xP^>6HDqK2xfb6_Jb)@~#qAW#QH zs6lr)2r3dTVfZVaS2sI86DdLzT6Nf^lWg8!ePnhLV5W)y&kc8Hz6L=DI5^k(;&Nn@ zdBJf85T>0A@i4WXkv}RPG)HdSj@|6a&nKv!(QK@xNb7O2M@aubA^EPxKqBC3ucQ>M z$yBV<+yt0GUPvTcS+%&t3k3z6TOKyH68HD^4t@hrC~3TKl8;x2pI;}aU07Hc*t@L_ z4d2-l!Fqw;!JY#XQ`Z`UuilQ(-d@+5Sq&V^9oVHyI)Dj~ieu7RXeNa;&L;C+9!nV- zXE<6y#oC4;8dYhImXGH_PL5i@=0j{wPR@3~{Pptf)d|e050pQDF2_uRumKUjIfoxb zMMecX?WoGUg3earD6`H0+G-| zsMQ_{5)t835fKjSB7F;%t)Az6ggyV-|Mw3AiZ71g3%NZCte5&YCG%>_TF*SBg5w<( zIV%Jakx{=|d@c^qDS?9TJ9&{)qoT0z@Zw@(5L^EI`2#jxT!>rHJ^g$puv9pd9)^^4 zi-}FuaprdwZOdi+sC$Y;0YgcnVZZ$O^Ct*VgaD(}*WUm``Q>FzZ7n_l!S&^(ufIP^ za(J%7cA6s-&uDSx>Y0RNxQWikxbOtL=Y7DMBua)N$iq_u7FbnPRcmW&9J3BO#(xD3 z4Umla+1c6E)fajx?C-8G>L!9qxciTf+af6fLpxFgOcsdllK0FRylmNl1M6d8Wb?<(MFdQpo-CzhMBzd4O4OuCJTS%>fS8*Vm_G!-)L# zhUG(al)z!f(8|s-dE3_jTz{vT&44xTH zQ>JJ8mmBMS^>f@yOaGzWxU}LkU+)(#3dQ*d5}ycCAx7l$>DGQXR1yU zutOS&yL)>+-rh(UBn%QPFv-;;gUg$h0b(AWrGFx0g&{JDpK}HtxXR~);BWpPu;4N;uKwi_ z;6&w8IayB*jv7oP)D#JP=WJ4|t>~?_4AdQgw6$OXS*E6={h>o9WC9jDE-${J}ULxq)PqAmMGDF%Bpe4x3vZIYG3S->6pan zO2ji+IZ~FiIaz^$mO;6aKPs8y1 z`?CJu)QDV~2Kh^{PNwTvGc#bUsO(2}J1P{W@uFCX3<)Yv}j?EG+Z zetGAV8N&w*%S$<_L$&2LPl6==fx*liN)c$!m9KHZ3L^k(1|NT_W;oraStxxW%O4GJ%ql9E z*XFNCUPVd(6+H6;Ks>=a+sdmO2LIY7=jZDFJ(~&zWp1v-K)kkJ`~=XBYCPodcOYkxkO!N+4yV$rwe2Yq(&hqj zY3u@z6FReYL+N5T@!=7K=27YcdsdXy9M|(6KLDMK8_(4A@aW@C)B2IfMtdJ>0Q$da z_PIL0Rv!q~GcT|7cFaiF3Yjap-{CT%=jD^qQWh>G1R0r7(uiFP1`l4xlinG%TFm?V zl9CqVRJ|P)6LX+|gRueKC7Ww_i4jre=+!xhq2vC&wO3|bI2 zx7*sD>gm-~P%K<`Jpm|1K1=DP(IF)zB|bj>`ue&sGTC)>a4?j2rrlp}iUk|{gCdqb ztBnqhpWkCxmiz>6Iw|DYcLw9^@9$4cOaRs^I4Fp=nUjIxZ(Eyit_3tz!j{$Bo5hpg zXrN=d{6o8cjpyChWOH+o+vC)AKU5nB94=~zvI02y=~W7T3wqQgCu^k&dh9PW&gRNz zJaT2k!JAmYK_A=w!M5SyX=y>gqDs3kO#F4O#cKFsUb6c$=LnM4ne@}yus$5cVc2< z+Hjy?lX%3>q|oA~r&e<3LLcN%Sba!Y^*eT_%ST5?0gf~?G<^Qy9VTEMEdi_r+#WDf;_$v5?lv?H5_kHR zl|iJ8PV5b*Luhb@>~OZLQxmg5z*?Yv3DALDJlogf)gD4F3xJWFIdfsKGSkuQ{zUEo zKR=*lZf2CLNaaz7sOdKp6wLo{Pc|#(Q1u|2VX|N$SD*p|4Q~5|;4#a}$rUT-PiQ3i zqZ1Pl61v@-TLPxW1M}#pb$}k7SpMj= zdimMS?e*!=$<3mtPh7exGno%&q+Q-~L?N4{58z@U#`I`=y8)P&gP9pX(#6qV2GWiA zz8o59^FbaOLik?5aJn0WV;LcHll^ia=|?fM-M6EYL4?Eo=0xY z`t6&=jPF`sWa>k~N4{qQ_&Gbk3-ibs^6)rMmFRigoD;JfN$<^6o}HaNG5q?qv!esB z*e$cOH{cf0ic~VyrQqBXfxM=HP{i|dlHvs|U3M>wAlD`MJS4&IF z_wV0q9E|;>aox8T8bXnq}cxsrOBv2@BFOPtwTSijyQI#k! zcXL4(R<%_!@bD~%P_-9nIbR++*oJgoTE$Bi>RekB&4`@ErOFo) zAOOS07>#r?ovso&90hm-z%ouNDk6FFrUpnC)m+oFv&e{uXe9h0t)qWpf~c;!DAsMCI)IwVsbLtV-f^7 z(t!b4Kt7N{bx0Lx+1g0nFem;fVF8x^X*P3pUHQnUf5+3mHzN3)($b&c;T>+HLr=G8 zC0m=pl}%7+K<4!zEQ;u%E|j|C=5W_53nTSJc>A!oH0P-5TFt7W$R1G2I#Y%7=$#=8 zRz8Hr^quWou3$t*r`Z}*SV!Q{`UvJ%)_?hF8o>|<<#9fwKU3do&N)RIc$`}}QPB$J zHM687i{T3-bf_6{r?DENH>mf#gA8vRyy=hk*|jh<(m%JjwH^G!vM;x~_7VhDIMKjT z(0xeDht<`$5lM@H3%10;dJBcmAuauReV^MHH6LkN*d2hDVuxH1u zBM7QiqjXpBmv859s`akZ3Q}lT>kT%RydSt0e~0(!lT8U-TW<_N48**;`LuATn*3WXXwgy6m#yUz3vd{>!usc+B&Ty_z ze`sn~EtY)+cFdtnsFEm9@&uKPsO;=O9u)cFnY7$s93$F&XXS34!`>{j9R9@aDrXcV zRj&3$pKdQV9Nyqg4%PryK}Q$}pB5bgc@ z;dux0vz6}1q+L7daU9xiJ-t>lYVual(RrL-qpus}_$}tg10C|-FkA2D&^c;hfoDmHeKd0Ul5~m&{zTrB#usl3|`h9iv_x%ZpIBtz5N%(iB z8b`FwPGBRhG(}&sWhq-<`a$;3Gtu-D^k1Oi?mSv2`lSV#mhy{5Fh06&R<3@Hr+AlByl!DK8{%t21SmZnW&d-Wiyn>KnzN6Oy%3Aa)x}7Me+lZ*ByR!lBzbz@!-)}622^di5THM;D)^xh< zGrX=0)b^`H9j@oL`!<{l0RSQ(pyVI0!sDZ+z;nxuOI_Ia=YG_>O#eA(Hjow)5Eyk# zl^Y`zAjr-JY~{D=)0+YQ4?>OX6Gb>2%zznh7dcfDMMO53sd&u%G`Pp)Jt~~-{F8~h z!bAelfqTvELf^3b#`zdplL-PTqVb_pB{m*#P9gmIrK^*Hz~0np-Td$B8_4-TMQy*? z=8qUnsU|sDk8vA_qHged1ySZrYE_K+8RKW&^zm$)h2^sOe4o82uYKp)Wsi4B5W|;X zY1_2G7j(hFCgbG9W`#5A^0PB4eAv*|@nChF)+*SY`3UFt8V4y$=ORI^?DZ{8kn2uyWo}TdPZ(+ zyjFpZFjVyV@X(H)ImdFAhUF)hC2ZGXrZOZcp@^~JCUPmi*8VhiDMSu$F?p@G(UB?R z1SYDaB?FC4l~c_5RqMG|g0W|MpuIg7L`2*(lkm5HYf>lX&KCh;2AONPU23=EWy1TY z1~Hb#zkSwJTON0u#c?iImtBhX=N2am`{(DYYn&rA9hx~!X}Se@w5H~YBG3BcMLt&F znatUtpl0++_6|}L5)x^UAWS`fg25d13sX{xLQ!#hB1}@!t~Rlu^L-6>ma*b! zJ0Ta>NQ(VNWP6tV9q(XnR8Ro%7xiD+RK$YPMB&FNuMt|jP@iDHrN~Y1ee^a9t0j~i zgmM-yzr9zv;H^hK=fQ0HX!M6RgHa=4vUz;mwC5LVEcZBJzIYg{`5xx|l-Dh`p)n7c zx!>G7@cba>e<|j?f|~jkH68^4MVbPlR8cz8YiJgbVnPcEQlkU}q)Cw?q4}5Iq=^U; z=_K^tK}A4{AiWc$cR~}8c9-Wo-I+Uc@7$M@2WHrrowfH~Wv{iq-}lW}ghghRO8nVS z$gd!S)ZvGC{}2d-uCA`IuuI%mLRG<5;n>(%hvOyajb>c~#sVXry5yGLfhXpiJwCHE zf^}w$&&1x9ZL_*Kr`6Yot`48UpgEb%CIx}S^>u}4X5u0=-Yh%#E{l24| zy>qv@xoK!hOZ2+S&O??S6RwcN#l?Z(RqwWZrCBVDJGyOWch?jkC2wy&h#z4+&dJn= zy1a>v)i6nZ7?s)RBjEY8zy0F{mBE7=Xgo4rL*`?Oq~xpYXsM_;mv@(GLPKiZvrdkP zRlB%UmD}Pme||eO=B2QOg{(STE;afFQmS+%Yfn1xOU6mDFIfuP+6rvRpNF_9CQu!G zRPX@64W=(s6X*K}Ur|v}SSZ4PG%(1_$x&BVr{h#sghI!LhYd%8d$zM9UA9fejnK87 zSr&M{_aLDtlC8+gQs==~brDVFIxiX7gLJE4TC)!<9%^#M@{{v&Y^Oi(b$ee$v#!q^ zQRH|yb#$cX@?Vmf9Nlb4mGZe`Wz~CfWND>_LTP!wfYO*}=C(aPh3a zTU>l@D5n5y+V6O9?=Ba5DBct^xt&O&q%8DV>@%zh#NNtSxYYdfAfN@~6sWzFd^qk; z>VlSu7-u=<6jh4I>G+*RvQI(os_icK^GQhbwYTGEIgNRSrUAVJP*_~K8<=N<1?tjw z-)W(o!M%q9AXH~V)!rSK4?Iczz04BE|Do}vHo7M@B&vsX%t zGJ4~jMcS^X=T`l4H+yM$rT}7TxL?Bq0@080uPQE<^E<5{zp7p1voRqL7ztX(i*p=` zpTZ=Pd^`$uM@&NxD`sGBMlYFWXYePDO`*CPL=pTbeTUo@;sR_|S`%lBu*?XS)jsnt zBM`T~Z$?_`KRX-hVrc2I;Zwt1l=o(EQr2fDR63 zX6zJuDpvKdzGmSJ&Mv;g`J-n#=|dm+__vx!<>%;DM52h4hhtjyMz$oHmyit_VeIyx$zr0Hm|O>(kN6e?3miYP73O@CeMVR($nD~WNO6PBSqTM(twKq7|l6hMT&U5RQi>qCCeGmn-L`h{TcWfnXVROq%YFtQDcXq3kMrSz#wMmt=mD_ETXri3T+p=ya!ib-*Yxq+5Mzw5dj} z=&kiQ9BLO?fb?r_rFfQs?)w#YqvfS>6tC5eVOS8XdpXDa zh{$$`IPj#^RA*21TdBC`t<2ZLcL6htmkeaEQ!6lk6J#zDWU$?8S8L({W>t86Xr)6& zfw$T<>v%T@4bBq}F?p!IHmzoazyjCH^)I&}oVX;-zB(bmn!zo83v-ab!mRVlxAK1O$fMgoz zUhB`lM_{hGqQP&m15Rr7f%$m-XC;O+ikViZzzX`s} z*z1fjZ|@1H@*}zvO$%&ZZ)|q5?Zd9JOa4k+Gf&!Dok`}cw@}#`NpMu3l00mXELtK} zRW;slOJm(Ut&J(Evpy@yMdyoG9W3?~)~@6?tj`kdY(Q)t6pc(CtnR7t&a-MAf0?LxB{ z%vRtkXOE>-!_Q724e5a(DuC_9Z~p}2w+le!m0#`J82sx4f}H@tpfGi~LypBNVxG08 z$oFGQir-#WIsLTHYhuECGp*)bVq3(G7IN}>U8k$UqJKue4bB<1i}I>YEq(2DUNf$< zxNs{*3e6`TBQ9+fahOG!vw^^V(l}eTw{}bmF=1$0;=dKw9xf|;`tJ5|%KZ4PTU1M-&pwW^3jP#U+_o&D1<@3k0*ZWJGKXhR^G_^$(hD%2b zkZm_dKbDQ&32%|5_Pkx=T9av2a;hy7Y`Hv~XTluqDC)AnJGR1Dzx4{?DOcvggCQan z;j4H1x8k}u{luoG4qLTp8vcj}1i+Ak{Y~#75CYjbB10sbZt$b?Y z)oHB9kdf0iFra9S*%zI=;Tqn;w09@GlL#C5HGJY{DVr<}9XBn=9sEI*P$HVGYv-j0 zYtL|foI|y2Jon%=Hn%}a%`;(tVY88X!3d@W)?otwr`AcdIX3fl3`xVq1+#KX+WBf? zW{&LSbE9KBr$+G{o3~~$Gm_&qt}w&N(%B*h#-3CAy6$=}Gr9}OPG_AY_(Bzt1p=wm zwfLO(&4c8yLn2jsZl$X_kM!}?dF?!7B@W3-Nev7QDUV}_)%ViUq7BL`oFtT#%0;iz zAO<2LzITL4NaA&H)^*-Bu|#lj{={)Dugk2Nr;Xv~I6Ghfp%LEXL_twWQeBRRFh!R+CLFJ3V4 zy^_m1YO)4Z)ll8?8}K8YJDGmTwdm4lBIu3T)MjEHlqeXN%+95+Hx`8jKlZC0qu=qkAyQPV;((dVs6 z7hZKLc*VIhvYJO0!|7|)D9ZmOKG~`3j`Hb)ork9%1ea6J83x;#RW3{(bqVb~jFpm2 z9A$b(L_RKlIa1gY;cj9RkCowd?VYPTfpc7DHtmS6K~>l;x~#u6hb9KCn4&*zJ?M?y z{D2ORk52)1_TY63N1|_cA_ZicCm6o_!(hWw^4`6<4Q&hYA0O1X89vHCxOeZST_b0G zDjKB*76tE^Kg!+9Gej5+?BHG*z}0KI+$O4S+zys~DwrDxq#AUcws_dVUWo5^z-dD0#oQUtu2=G zq76XRMjcUl_io;F+rEu7&?Po9GOChM#3$%!P_?w^IozqK)qf!-FIR8lZQCX>qKwFvtI zQ^((@kMx8#t>$cvk_MV~l46?Gn?v|lb^E3&_BrOL^Y0{2nLTq_nAcFv61b4l^r*c- zv}7|w{JFDVKq@bTrVv4rc7ngwq~E6+Ut+i!>@YZ}XT?*ABZDY%@^BDPJm<}wOOu>u zTTgZ@B>~3dU6*3NAWy~tflPZfthA-00bBq!l}^1%uLC|qMWv|)=om(&oENI1;{k ztE@2upul{5d{wjT%2{)Ba}yJ$Je&P=I$Q}_;bCF05}{87L4JNQDXF;~c6xf6G~*qw@&%MsWLl$-LjC5P^kW4PSU^T_c=F`u z79GtQAua8utnBdEnDJ=j=UK1@==t_dm5fh7S95-T9zfz}z72*W<4#QuR#sW@!Kl!a zX(LfCuD^T8L?RlTg2~Lz)=rTV=iyQN`vWy~b%1&7{qY0-_dDP{^5iBtR(*Qs@O;4C zQM1ppTZdRS9(cmjy95{){S&AlHD&EG6<9YZa_I0g3f3`r*7x=4vsvp}B_t$-s&dl? zXRJM(76%#>;CF!cDt?kmk>cNF6o>u~MJa-Esc%l7gFqfr{a>%{=n)LIySod(P$@~t zWEHszkpFVNf+1UZmo02lP!fseAOTR@QrY?}F zZ~NkfSYjs)9i3A}Ze-|_$LtH5WU*RVK%D`Ro1U_G5sIqv zMMaz}Ec|Fj28NsTR5JD>FMfRmn_Y)No%0+D?y#E4#^=Y!$1phjFNp(U@qh|wYhxp9 zR%_5m)vp&;U$3a6qXSl%$~aDD<{=k7)~jENp)gjmOP4O`8ycbk9yD5J4FUz*=Us*i zO!&}m&f%%Zu?k0XaXMfsW+{<$JOiwF03n+51MHsf-k@QaUE@x-0??SWY_^cVgj-Qw z4sbKZ1Ult|&r?(8r?S}@jH^Qh$zo6MzkFE%fbP6w8w@uKD%h*f|HFq5z@GUxku4XN zB3FWL^4)F{=VA(g@2;-_FEP+nkSTSq=Dx=OwvHzqgp{tW*@MT>B8~Pt*#k@x5{VSZ z={2MaZhp!V)Uzlg2n4sw&GgwBiU10ns&c|rR=VfC`qw@7#Zi2lcPTQfTiDzz%6)8X z%zW+I)hRfE{2Vd0qoZScdmE%4{nb%zY-}7IH#Rqk0|Tn6K>k*>#9%V2ODig*o5g^r zsAad$)cotkv!2ujz%l_F3AdJp#z>RjS-w#jc)}di)ZuMI3mzPBYEW(RBEO#P=Y2-R`T3iH{RfbL*?7F&2|Om`)-|0&pV>hRi4l9KZ ztxe<)7mZqEL=K9SJV{)_!opHgyB|HKoLyW#DEs^Yjtx})V?8so$mQG*GQ6*!;%Sp$ zGTy(}K=knY_be@m8!T^E?tZS~qo$??>8OW6ujCJ-(m+ja_f#$~FNdjtKZ+bZ-Q9(= zM_6fK7U{5*04-RpWoQT`%xrCKje!!f|8^e&2KaX=={O?7!&&8BGSfrO;c3aS-ZNXC zQ&UsGSFCFQpzqziJ%EzigO7s)9XmAhF(W&>lH1Zxzf(n~(_5gs@T3FPPCyV8%0`<5 z+SvOjf-PvT@4`laHyDuA7!2kJIP^=Sz^`K|iQTXN0SU1jimnSV@csb%OTfu)3qS}{ zk;D){NUtm|f_@}RM;8so^|sd5)5Aqv$HqRB1qLJeYIjL3oEa1(4~6pb@d4A&vp>7b z+h+>zqx4Oie9ap7f%m(!s|(!EEp87Fk2Kv8P}TXgJDKOiFJEP3oSU6RAP~1iJ7{i5 zX~5w$Wco}U;}>uh;I~17Ey&LWNNnOA6W#U3}bqiWt-He?!4N|=N|ati?`vA;IDlZf`k z=jC;FcY}yqq`NwVR1RKfBvr6R1qL>B^@}uBZV+QECL@FT*O!0aL}U7Bw3CGeH$4>? zeg*~y6DDaPig%8`(?Z5d-In!TU0p#27JwRr>aIaNDIWhr>CaaSGVp-p4x)p3n8SUr zSg=_KUjFI-x*o3SbJJU3Nd2?FNd!w}b_I9(USM_W>+357gakcR{NxSDv4VX%rGJ}A z*)zDdcK4sX)%k19&CQ#Wbyd#?1?Sx_Ky)tJ{ktxu|F<{vf27L5)ZvV>QSZsc+l}>7 Qf7_m>nzm};{U@*g3u%lBxc~qF literal 0 HcmV?d00001 diff --git a/tests/05-mixing.cpp b/tests/05-mixing.cpp new file mode 100644 index 000000000..67396a25c --- /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 000000000..18abce106 --- /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 b878ffbb6..1d4bb2cac 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 510a6bfaa..3757a5566 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 dde5d1536..c72101439 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