diff --git a/examples/10-likelihood-free-mcmc/main.cpp b/examples/10-likelihood-free-mcmc/main.cpp index af2ba1819..a5d55c4a5 100644 --- a/examples/10-likelihood-free-mcmc/main.cpp +++ b/examples/10-likelihood-free-mcmc/main.cpp @@ -1,7 +1,11 @@ // #include "../../include/epiworld/epiworld.hpp" #include "../../include/epiworld/epiworld.hpp" + +#define EPI_DEBUG #include "../../include/epiworld/math/lfmcmc.hpp" + + using namespace epiworld; epimodels::ModelSIR<> model; @@ -53,7 +57,7 @@ int main() model, // Model "covid", // Name of the virus .1, // Initial prevalence - .9, // Infectiousness (par[1]) + .1, // Infectiousness (par[1]) .3 // Immune Recovery (par[0]) ); @@ -69,7 +73,6 @@ int main() lfmcmc.set_kernel_fun(kernel_fun_gaussian>); // Simulating some data - model.set_backup(); model.verbose_off(); model.run(50, 122); model.print(); @@ -82,7 +85,7 @@ int main() std::vector< epiworld_double > par0 = {.5, .5}; - lfmcmc.run(par0, 1000, 1); + lfmcmc.run(par0, 2000, 1); lfmcmc.set_par_names({"Immune recovery", "Infectiousness"}); lfmcmc.set_stats_names(model.get_states()); diff --git a/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp b/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp index ea57b9a40..688f4fd03 100755 --- a/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp +++ b/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp @@ -83,6 +83,13 @@ inline LFMCMCProposalFun make_proposal_norm_reflective( } + #ifdef EPI_DEBUG + for (auto & p : params_now) + if (p < lb || p > ub) + throw std::range_error("The parameter is out of bounds."); + #endif + + return; }; diff --git a/include/epiworld/model-meat.hpp b/include/epiworld/model-meat.hpp index 03075212c..5ad8172bc 100644 --- a/include/epiworld/model-meat.hpp +++ b/include/epiworld/model-meat.hpp @@ -2039,21 +2039,19 @@ inline void Model::reset() { } #endif - } - else - { - for (auto & p : population) - p.reset(); + } - #ifdef EPI_DEBUG - for (auto & a: population) - { - if (a.get_state() != 0u) - throw std::logic_error("Model::reset population doesn't match." - "Some agents are not in the baseline state."); - } - #endif + for (auto & p : population) + p.reset(); + + #ifdef EPI_DEBUG + for (auto & a: population) + { + if (a.get_state() != 0u) + throw std::logic_error("Model::reset population doesn't match." + "Some agents are not in the baseline state."); } + #endif if (entities_backup.size() != 0) { diff --git a/include/epiworld/models/seirconnected_logit.hpp b/include/epiworld/models/seirconnected_logit.hpp deleted file mode 100644 index ba871ac35..000000000 --- a/include/epiworld/models/seirconnected_logit.hpp +++ /dev/null @@ -1,333 +0,0 @@ -#ifndef EPIWORLD_MODELS_SEIRCONNECTEDLOGIT_HPP -#define EPIWORLD_MODELS_SEIRCONNECTEDLOGIT_HPP - -template -class ModelSEIRCONNLogit : public epiworld::Model -{ -private: - static const int SUSCEPTIBLE = 0; - static const int EXPOSED = 1; - static const int INFECTED = 2; - static const int RECOVERE = 3; - -public: - - ModelSEIRCONNLogit() { - tracked_agents_infected.reserve(1e4); - tracked_agents_infected_next.reserve(1e4); - }; - - ModelSEIRCONNLogit( - ModelSEIRCONNLogit & 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, - double * covars, - std::vector< double > logit_params - ); - - ModelSEIRCONNLogit( - 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 - double * covars, - std::vector< double > logit_params - ); - - // Tracking who is infected and who is not - std::vector< epiworld::Agent<>* > tracked_agents_infected = {}; - std::vector< epiworld::Agent<>* > tracked_agents_infected_next = {}; - - bool tracked_started = false; - int tracked_ninfected = 0; - int tracked_ninfected_next = 0; - -}; - -/** - * @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 Reproductive number (beta) - * @param transmission_rate Probability of transmission - * @param recovery_rate Probability of recovery - */ -template -inline ModelSEIRCONNLogit::ModelSEIRCONNLogit( - ModelSEIRCONNLogit & 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, - double * covars, - std::vector< double > logit_params - // epiworld_double prob_reinfection - ) -{ - - auto * _tracked_agents_infected = &model.tracked_agents_infected; - auto * _tracked_agents_infected_next = &model.tracked_agents_infected_next; - auto * _tracked_started = &model.tracked_started; - auto * _tracked_ninfected = &model.tracked_ninfected; - auto * _tracked_ninfected_next = &model.tracked_ninfected_next; - - std::function *)> tracked_agents_check_init = - [ - _tracked_started, - _tracked_agents_infected, - _tracked_ninfected - ](epiworld::Model * m) - { - - if (*_tracked_started) - return; - - /* Checking first if it hasn't */ - if (!*_tracked_started) - { - - /* Listing who is infected */ - for (auto & p : m->get_agents()) - { - if (p.get_state() == ModelSEIRCONNLogit::INFECTED) - { - - _tracked_agents_infected->push_back(&p); - *_tracked_ninfected += 1; - - } - } - - for (auto & p: *_tracked_agents_infected) - { - if (p->get_virus() == nullptr) - throw std::logic_error("Cannot be infected and have no viruses."); - } - - *_tracked_started = true; - - } - - }; - - epiworld::UpdateFun update_susceptible = - [ - tracked_agents_check_init, - _tracked_ninfected, - _tracked_agents_infected - ](epiworld::Agent * p, epiworld::Model * m) -> void - { - - tracked_agents_check_init(m); - - // No infected individual? - if (*_tracked_ninfected == 0) - return; - - // Computing probability of contagion - // P(infected) = 1 - (1 - beta/Pop * ptransmit) ^ ninfected - epiworld_double prob_infect = 1.0 - std::pow( - 1.0 - (m->par("Contact rate")) * (m->par("Transmission rate")) / m->size(), - *_tracked_ninfected - ); - - if (m->runif() < prob_infect) - { - - // Now selecting who is transmitting the disease - epiworld_fast_uint which = static_cast( - std::floor(*_tracked_ninfected * m->runif()) - ); - - // Infecting the individual - #ifdef EPI_DEBUG - if (_tracked_agents_infected->operator[](which)->get_virus() == nullptr) - { - - printf_epiworld("[epi-debug] date: %i\n", m->today()); - printf_epiworld("[epi-debug] sim#: %i\n", m->get_n_replicates()); - - throw std::logic_error( - "[epi-debug] The agent " + std::to_string(which) + " has no "+ - "virus to share. The agent's state is: " + - std::to_string(_tracked_agents_infected->operator[](which)->get_state()) - ); - } - #endif - p->set_virus( - _tracked_agents_infected->operator[](which)->get_virus(), - ModelSEIRCONNLogit::EXPOSED - ); - - return; - - } - - return; - - }; - - epiworld::UpdateFun update_infected = - [ - tracked_agents_check_init, - _tracked_agents_infected_next, - _tracked_ninfected_next - - ](epiworld::Agent * p, epiworld::Model * m) -> void - { - - tracked_agents_check_init(m); - auto state = p->get_state(); - - if (state == ModelSEIRCONNLogit::EXPOSED) - { - - // Does the agent become infected? - if (m->runif() < 1.0/(m->par("Avg. Incubation days"))) - { - // Adding the individual to the queue - _tracked_agents_infected_next->push_back(p); - *_tracked_ninfected_next += 1; - - p->change_state(m, ModelSEIRCONNLogit::INFECTED); - - return; - - } - - - } else if (state == ModelSEIRCONNLogit::INFECTED) - { - - if (m->runif() < (m->par("Recovery rate"))) - { - - *_tracked_ninfected_next -= 1; - p->rm_virus(m); - return; - - } - - _tracked_agents_infected_next->push_back(p); - - } - - return; - - }; - - epiworld::GlobalFun global_accounting = - [ - _tracked_started, - _tracked_agents_infected, - _tracked_agents_infected_next, - _tracked_ninfected, - _tracked_ninfected_next - ](epiworld::Model* m) -> void - { - - // On the last day, also reset tracked agents and - // set the initialized value to false - if (static_cast(m->today()) == (m->get_ndays() - 1)) - { - - *_tracked_started = false; - _tracked_agents_infected->clear(); - _tracked_agents_infected_next->clear(); - *_tracked_ninfected = 0; - *_tracked_ninfected_next = 0; - - return; - } - - std::swap(*_tracked_agents_infected, *_tracked_agents_infected_next); - _tracked_agents_infected_next->clear(); - - *_tracked_ninfected += *_tracked_ninfected_next; - *_tracked_ninfected_next = 0; - - }; - - // Setting up parameters - model.add_param(contact_rate, "Contact rate"); - model.add_param(transmission_rate, "Transmission rate"); - model.add_param(recovery_rate, "Recovery rate"); - 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"); - - // Adding agent's parameters - model.set_agents_data(covars, logit_params.size()); - for (size_t i = 0u; i < logit_params.size(); ++i) - model.add_param( - std::string("x" + std::to_string(i)), - logit_params[i] - ); - - // Preparing the virus ------------------------------------------- - epiworld::Virus virus(vname); - virus.set_state(1,3,3); - model.add_virus(virus, prevalence); - - // Adding updating function - model.add_globalevent(global_accounting, "Accounting", -1); - - model.queuing_off(); // No queuing need - - // Adding the empty population - model.agents_empty_graph(n); - - model.set_name("Susceptible-Exposed-Infected-Removed (SEIR) (connected)"); - - return; - -} - -template -inline ModelSEIRCONNLogit::ModelSEIRCONNLogit( - 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, - double * covars, - std::vector< double > logit_params - ) -{ - - ModelSEIRCONNLogit( - *this, - vname, - n, - prevalence, - contact_rate, - transmission_rate, - avg_incubation_days, - covars, - logit_params - ); - - return; - -} - -#endif diff --git a/include/epiworld/models/sirlogit.hpp b/include/epiworld/models/sirlogit.hpp index 81ed8616b..9af13166f 100644 --- a/include/epiworld/models/sirlogit.hpp +++ b/include/epiworld/models/sirlogit.hpp @@ -3,6 +3,34 @@ #ifndef EPIWORLD_MODELS_SIRLOGIT_HPP #define EPIWORLD_MODELS_SIRLOGIT_HPP + +/** + * @brief Template for a Susceptible-Infected-Removed (SIR) model + * + * @details + * In this model, infection and recoveru probabilities are computed + * using a logit model. Particularly, the probability of infection + * is computed as: + * + * \f[ + * \frac{1}{1 + \exp\left(-\left(\beta_0 E_i + \sum_{i=1}^{n} \beta_i x_i\right)\right)} + * \f] + * + * where \f$\beta_0\f$ is the exposure coefficient and \f$E_i\f$ is the exposure + * number, \f$\beta_i\f$ are the + * coefficients for the features \f$x_i\f$ of the agents, and \f$n\f$ is the + * number of features. The probability of recovery is computed as: + * + * \f[ + * \frac{1}{1 + \exp\left(-\left(\sum_{i=1}^{n} \beta_i x_i\right)\right)} + * \f] + * + * where \f$\beta_i\f$ are the coefficients for the features \f$x_i\f$ of the agents, + * and \f$n\f$ is the number of features. + * + * @param TSeq Type of the sequence (e.g. std::vector, std::deque) + +*/ template class ModelSIRLogit : public epiworld::Model {