From 51ab550dbe3a019077e021c440dab16dbec2f6eb Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Wed, 13 Aug 2025 16:26:28 +0200 Subject: [PATCH 01/19] CHG: functional smm example with agegroups, rest is broken --- cpp/examples/smm.cpp | 51 +++++++++++++++--------- cpp/memilio/epidemiology/adoption_rate.h | 8 ++-- cpp/models/smm/model.h | 33 ++++++++------- cpp/models/smm/parameters.h | 16 ++++---- cpp/models/smm/simulation.h | 34 +++++++++------- 5 files changed, 84 insertions(+), 58 deletions(-) diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index 1e7f5d285e..5ee794e6d5 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -18,6 +18,7 @@ * limitations under the License. */ +#include "memilio/epidemiology/age_group.h" #include "smm/simulation.h" #include "smm/model.h" #include "smm/parameters.h" @@ -36,41 +37,53 @@ enum class InfectionState }; +enum class AgeGroup +{ + All, + Count +}; + int main() { //Example how to run the stochastic metapopulation models with four regions const size_t num_regions = 4; - using Model = mio::smm::Model; + using Model = mio::smm::Model; double numE = 12, numC = 4, numI = 12, numR = 0, numD = 0; Model model; //Population are distributed uniformly to the four regions for (size_t r = 0; r < num_regions; ++r) { - model.populations[{mio::regions::Region(r), InfectionState::S}] = + model.populations[{mio::regions::Region(r), InfectionState::S, AgeGroup::All}] = (1000 - numE - numC - numI - numR - numD) / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::E}] = numE / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::C}] = numC / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::I}] = numI / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::R}] = numR / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::D}] = numD / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::E, AgeGroup::All}] = numE / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::C, AgeGroup::All}] = numC / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::I, AgeGroup::All}] = numI / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::R, AgeGroup::All}] = numR / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::D, AgeGroup::All}] = numD / num_regions; } //Set infection state adoption and spatial transition rates - std::vector> adoption_rates; - std::vector> transition_rates; + std::vector> adoption_rates; + std::vector> transition_rates; for (size_t r = 0; r < num_regions; ++r) { adoption_rates.push_back({InfectionState::S, InfectionState::E, + AgeGroup::All, mio::regions::Region(r), 0.1, - {{InfectionState::C, 1}, {InfectionState::I, 0.5}}}); - adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}}); - adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}}); - adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}}); - adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}}); - adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}}); + {{InfectionState::C, AgeGroup::All, 1}, {InfectionState::I, AgeGroup::All, 0.5}}}); + adoption_rates.push_back( + {InfectionState::E, InfectionState::C, AgeGroup::All, mio::regions::Region(r), 1.0 / 5., {}}); + adoption_rates.push_back( + {InfectionState::C, InfectionState::R, AgeGroup::All, mio::regions::Region(r), 0.2 / 3., {}}); + adoption_rates.push_back( + {InfectionState::C, InfectionState::I, AgeGroup::All, mio::regions::Region(r), 0.8 / 3., {}}); + adoption_rates.push_back( + {InfectionState::I, InfectionState::R, AgeGroup::All, mio::regions::Region(r), 0.99 / 5., {}}); + adoption_rates.push_back( + {InfectionState::I, InfectionState::D, AgeGroup::All, mio::regions::Region(r), 0.01 / 5., {}}); } //Agents in infection state D do not transition @@ -79,15 +92,15 @@ int main() for (size_t j = 0; j < num_regions; ++j) if (i != j) { transition_rates.push_back( - {InfectionState(s), mio::regions::Region(i), mio::regions::Region(j), 0.01}); + {InfectionState(s), AgeGroup::All, mio::regions::Region(i), mio::regions::Region(j), 0.01}); transition_rates.push_back( - {InfectionState(s), mio::regions::Region(j), mio::regions::Region(i), 0.01}); + {InfectionState(s), AgeGroup::All, mio::regions::Region(j), mio::regions::Region(i), 0.01}); } } } - model.parameters.get>() = adoption_rates; - model.parameters.get>() = transition_rates; + model.parameters.get>() = adoption_rates; + model.parameters.get>() = transition_rates; double dt = 0.1; double tmax = 30.; diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h index a56dfcdc54..43fc44aedf 100644 --- a/cpp/memilio/epidemiology/adoption_rate.h +++ b/cpp/memilio/epidemiology/adoption_rate.h @@ -32,9 +32,10 @@ namespace mio * The population having "status" is multiplied with "factor." * @tparam Status An infection state enum. */ -template +template struct Influence { Status status; + AgeGroup age_group; ScalarType factor; }; @@ -45,13 +46,14 @@ struct Influence { * the sum over all "influences", which scale their "status" with the respective "factor". * @tparam Status An infection state enum. */ -template +template struct AdoptionRate { Status from; // i Status to; // j + AgeGroup age_group; mio::regions::Region region; // k ScalarType factor; // gammahat_{ij}^k - std::vector> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) + std::vector> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) }; } // namespace mio diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index 4066667455..4c41732fc8 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -37,17 +37,18 @@ namespace smm * @tparam regions Number of regions. * @tparam Status An infection state enum. */ -template -class Model - : public mio::CompartmentalModel, - ParametersBase> +template +class Model : public mio::CompartmentalModel, + ParametersBase> { - using Base = mio::CompartmentalModel, - ParametersBase>; + using Base = mio::CompartmentalModel, + ParametersBase>; public: Model() - : Base(typename Base::Populations({static_cast(regions), Status::Count}, 0.), + : Base(typename Base::Populations({static_cast(regions), Status::Count, AgeGroup::Count}), typename Base::ParameterSet()) { } @@ -58,10 +59,11 @@ class Model * @param[in] x The current state of the model. * @return Current value of the adoption rate. */ - ScalarType evaluate(const AdoptionRate& rate, const Eigen::VectorXd& x) const + ScalarType evaluate(const AdoptionRate& rate, const Eigen::VectorXd& x) const { - const auto& pop = this->populations; - const auto source = pop.get_flat_index({rate.region, rate.from}); + const auto& pop = this->populations; + const auto source = + pop.get_flat_index({rate.region, rate.from, rate.age_group}); // Why is here rate.from used? KV // determine order and calculate rate if (rate.influences.size() == 0) { // first order adoption return rate.factor * x[source]; @@ -69,13 +71,16 @@ class Model else { // second order adoption ScalarType N = 0; for (size_t s = 0; s < static_cast(Status::Count); ++s) { - N += x[pop.get_flat_index({rate.region, Status(s)})]; + for (size_t age = 0; age < static_cast(AgeGroup::Count); ++age) { + N += x[pop.get_flat_index({rate.region, Status(s), AgeGroup(age)})]; + } } // accumulate influences ScalarType influences = 0.0; for (size_t i = 0; i < rate.influences.size(); i++) { influences += - rate.influences[i].factor * x[pop.get_flat_index({rate.region, rate.influences[i].status})]; + rate.influences[i].factor * + x[pop.get_flat_index({rate.region, rate.influences[i].status, rate.influences[i].age_group})]; } return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; } @@ -87,9 +92,9 @@ class Model * @param[in] x The current state of the model. * @return Current value of the transition rate. */ - ScalarType evaluate(const TransitionRate& rate, const Eigen::VectorXd& x) const + ScalarType evaluate(const TransitionRate& rate, const Eigen::VectorXd& x) const { - const auto source = this->populations.get_flat_index({rate.from, rate.status}); + const auto source = this->populations.get_flat_index({rate.from, rate.status, rate.age_group}); return rate.factor * x[source]; } diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h index 70b8212291..cc6e734851 100644 --- a/cpp/models/smm/parameters.h +++ b/cpp/models/smm/parameters.h @@ -35,9 +35,9 @@ namespace smm * @brief A vector of AdoptionRate%s, see mio::AdoptionRate * @tparam Status An infection state enum. */ -template +template struct AdoptionRates { - using Type = std::vector>; + using Type = std::vector>; const static std::string name() { return "AdoptionRates"; @@ -48,25 +48,25 @@ struct AdoptionRates { * @brief Struct defining a possible regional transition in a Model based on Poisson Processes. * @tparam Status An infection state enum. */ -template +template struct TransitionRate { Status status; // i + AgeGroup age_group; mio::regions::Region from; // k mio::regions::Region to; // l ScalarType factor; // lambda_i^{kl} }; - -template +template struct TransitionRates { - using Type = std::vector>; + using Type = std::vector>; const static std::string name() { return "TransitionRates"; } }; -template -using ParametersBase = mio::ParameterSet, TransitionRates>; +template +using ParametersBase = mio::ParameterSet, TransitionRates>; } // namespace smm diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index b63e77fc99..4cf13abff5 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -37,12 +37,12 @@ namespace smm * @tparam regions The number of regions. * @tparam Status An infection state enum. */ -template +template class Simulation { public: public: - using Model = smm::Model; + using Model = smm::Model; /** * @brief Set up the simulation for a Stochastic Metapopulation Model. @@ -100,18 +100,24 @@ class Simulation if (next_event < adoption_rates().size()) { // perform adoption event const auto& rate = adoption_rates()[next_event]; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from})] -= 1; - m_model->populations[{rate.region, rate.from}] -= 1; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to})] += 1; - m_model->populations[{rate.region, rate.to}] += 1; + m_result + .get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from, rate.age_group})] -= + 1; + m_model->populations[{rate.region, rate.from, rate.age_group}] -= 1; + m_result + .get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to, rate.age_group})] += 1; + m_model->populations[{rate.region, rate.to, rate.age_group}] += 1; } else { // perform transition event const auto& rate = transition_rates()[next_event - adoption_rates().size()]; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status})] -= 1; - m_model->populations[{rate.from, rate.status}] -= 1; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status})] += 1; - m_model->populations[{rate.to, rate.status}] += 1; + m_result + .get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status, rate.age_group})] -= + 1; + m_model->populations[{rate.from, rate.status, rate.age_group}] -= 1; + m_result + .get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status, rate.age_group})] += 1; + m_model->populations[{rate.to, rate.status, rate.age_group}] += 1; } // update internal times for (size_t i = 0; i < m_internal_time.size(); i++) { @@ -161,17 +167,17 @@ class Simulation /** * @brief Returns the model's transition rates. */ - inline constexpr const typename smm::TransitionRates::Type& transition_rates() + inline constexpr const typename smm::TransitionRates::Type& transition_rates() { - return m_model->parameters.template get>(); + return m_model->parameters.template get>(); } /** * @brief Returns the model's adoption rates. */ - inline constexpr const typename smm::AdoptionRates::Type& adoption_rates() + inline constexpr const typename smm::AdoptionRates::Type& adoption_rates() { - return m_model->parameters.template get>(); + return m_model->parameters.template get>(); } /** From d78bd7c364ac2a0c202ee1c7026877faae97da22 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Tue, 9 Sep 2025 13:54:46 +0200 Subject: [PATCH 02/19] CHG: functional base model with template parameter list (not used yet) --- cpp/examples/CMakeLists.txt | 6 ++-- cpp/memilio/epidemiology/adoption_rate.h | 8 ++--- cpp/models/smm/model.h | 30 +++++++++---------- cpp/models/smm/parameters.h | 15 +++++----- cpp/models/smm/simulation.h | 37 ++++++++++-------------- 5 files changed, 41 insertions(+), 55 deletions(-) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 980164d324..d97c595ad7 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -148,9 +148,9 @@ add_executable(dabm_example d_abm.cpp) target_link_libraries(dabm_example PRIVATE memilio d_abm) target_compile_options(dabm_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) -add_executable(smm_example smm.cpp) -target_link_libraries(smm_example PRIVATE memilio smm) -target_compile_options(smm_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +# add_executable(smm_example smm.cpp) +# target_link_libraries(smm_example PRIVATE memilio smm) +# target_compile_options(smm_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) add_executable(temporal_hybrid_example temporal_hybrid_dabm_osecir.cpp) target_link_libraries(temporal_hybrid_example PRIVATE memilio hybrid ode_secir d_abm) diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h index 43fc44aedf..b05ac2b0a8 100644 --- a/cpp/memilio/epidemiology/adoption_rate.h +++ b/cpp/memilio/epidemiology/adoption_rate.h @@ -32,10 +32,9 @@ namespace mio * The population having "status" is multiplied with "factor." * @tparam Status An infection state enum. */ -template +template struct Influence { Status status; - AgeGroup age_group; ScalarType factor; }; @@ -46,14 +45,13 @@ struct Influence { * the sum over all "influences", which scale their "status" with the respective "factor". * @tparam Status An infection state enum. */ -template +template struct AdoptionRate { Status from; // i Status to; // j - AgeGroup age_group; mio::regions::Region region; // k ScalarType factor; // gammahat_{ij}^k - std::vector> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) + std::vector> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) }; } // namespace mio diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index 4c41732fc8..4699fb7e27 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -37,18 +37,18 @@ namespace smm * @tparam regions Number of regions. * @tparam Status An infection state enum. */ -template +template class Model : public mio::CompartmentalModel, - ParametersBase> + mio::Populations, + ParametersBase> { using Base = mio::CompartmentalModel, - ParametersBase>; + mio::Populations, + ParametersBase>; public: Model() - : Base(typename Base::Populations({static_cast(regions), Status::Count, AgeGroup::Count}), + : Base(typename Base::Populations({static_cast(regions), Status::Count}), typename Base::ParameterSet()) { } @@ -59,11 +59,10 @@ class Model : public mio::CompartmentalModel& rate, const Eigen::VectorXd& x) const + ScalarType evaluate(const AdoptionRate& rate, const Eigen::VectorXd& x) const { - const auto& pop = this->populations; - const auto source = - pop.get_flat_index({rate.region, rate.from, rate.age_group}); // Why is here rate.from used? KV + const auto& pop = this->populations; + const auto source = pop.get_flat_index({rate.region, rate.from}); // Why is here rate.from used? KV // determine order and calculate rate if (rate.influences.size() == 0) { // first order adoption return rate.factor * x[source]; @@ -71,16 +70,13 @@ class Model : public mio::CompartmentalModel(Status::Count); ++s) { - for (size_t age = 0; age < static_cast(AgeGroup::Count); ++age) { - N += x[pop.get_flat_index({rate.region, Status(s), AgeGroup(age)})]; - } + N += x[pop.get_flat_index({rate.region, Status(s)})]; } // accumulate influences ScalarType influences = 0.0; for (size_t i = 0; i < rate.influences.size(); i++) { influences += - rate.influences[i].factor * - x[pop.get_flat_index({rate.region, rate.influences[i].status, rate.influences[i].age_group})]; + rate.influences[i].factor * x[pop.get_flat_index({rate.region, rate.influences[i].status})]; } return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; } @@ -92,9 +88,9 @@ class Model : public mio::CompartmentalModel& rate, const Eigen::VectorXd& x) const + ScalarType evaluate(const TransitionRate& rate, const Eigen::VectorXd& x) const { - const auto source = this->populations.get_flat_index({rate.from, rate.status, rate.age_group}); + const auto source = this->populations.get_flat_index({rate.from, rate.status}); return rate.factor * x[source]; } diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h index cc6e734851..8a251f63af 100644 --- a/cpp/models/smm/parameters.h +++ b/cpp/models/smm/parameters.h @@ -35,9 +35,9 @@ namespace smm * @brief A vector of AdoptionRate%s, see mio::AdoptionRate * @tparam Status An infection state enum. */ -template +template struct AdoptionRates { - using Type = std::vector>; + using Type = std::vector>; const static std::string name() { return "AdoptionRates"; @@ -48,25 +48,24 @@ struct AdoptionRates { * @brief Struct defining a possible regional transition in a Model based on Poisson Processes. * @tparam Status An infection state enum. */ -template +template struct TransitionRate { Status status; // i - AgeGroup age_group; mio::regions::Region from; // k mio::regions::Region to; // l ScalarType factor; // lambda_i^{kl} }; -template +template struct TransitionRates { - using Type = std::vector>; + using Type = std::vector>; const static std::string name() { return "TransitionRates"; } }; -template -using ParametersBase = mio::ParameterSet, TransitionRates>; +template +using ParametersBase = mio::ParameterSet, TransitionRates>; } // namespace smm diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 4cf13abff5..01afaa1689 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -37,12 +37,11 @@ namespace smm * @tparam regions The number of regions. * @tparam Status An infection state enum. */ -template +template class Simulation { public: -public: - using Model = smm::Model; + using Model = smm::Model; /** * @brief Set up the simulation for a Stochastic Metapopulation Model. @@ -100,24 +99,18 @@ class Simulation if (next_event < adoption_rates().size()) { // perform adoption event const auto& rate = adoption_rates()[next_event]; - m_result - .get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from, rate.age_group})] -= - 1; - m_model->populations[{rate.region, rate.from, rate.age_group}] -= 1; - m_result - .get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to, rate.age_group})] += 1; - m_model->populations[{rate.region, rate.to, rate.age_group}] += 1; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from})] -= 1; + m_model->populations[{rate.region, rate.from}] -= 1; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to})] += 1; + m_model->populations[{rate.region, rate.to}] += 1; } else { // perform transition event const auto& rate = transition_rates()[next_event - adoption_rates().size()]; - m_result - .get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status, rate.age_group})] -= - 1; - m_model->populations[{rate.from, rate.status, rate.age_group}] -= 1; - m_result - .get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status, rate.age_group})] += 1; - m_model->populations[{rate.to, rate.status, rate.age_group}] += 1; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status})] -= 1; + m_model->populations[{rate.from, rate.status}] -= 1; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status})] += 1; + m_model->populations[{rate.to, rate.status}] += 1; } // update internal times for (size_t i = 0; i < m_internal_time.size(); i++) { @@ -167,17 +160,17 @@ class Simulation /** * @brief Returns the model's transition rates. */ - inline constexpr const typename smm::TransitionRates::Type& transition_rates() + inline constexpr const typename smm::TransitionRates::Type& transition_rates() { - return m_model->parameters.template get>(); + return m_model->parameters.template get>(); } /** * @brief Returns the model's adoption rates. */ - inline constexpr const typename smm::AdoptionRates::Type& adoption_rates() + inline constexpr const typename smm::AdoptionRates::Type& adoption_rates() { - return m_model->parameters.template get>(); + return m_model->parameters.template get>(); } /** @@ -221,7 +214,7 @@ class Simulation m_current_rates; ///< Current values of both types of rates i.e. adoption and transition rates. }; -} //namespace smm +} // namespace smm } // namespace mio #endif From f2babeee69b842cdcfffc7584144adfd0a7a3215 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Tue, 9 Sep 2025 16:58:24 +0200 Subject: [PATCH 03/19] CHG: Another compiling status that is not leading anywhere --- cpp/examples/CMakeLists.txt | 6 +- cpp/examples/smm.cpp | 72 +++++----- cpp/memilio/epidemiology/adoption_rate.h | 2 + cpp/models/smm/model.h | 100 +++++++++---- cpp/models/smm/parameters.h | 4 +- cpp/models/smm/simulation.h | 170 ++++++++++++++++------- 6 files changed, 238 insertions(+), 116 deletions(-) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index d97c595ad7..980164d324 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -148,9 +148,9 @@ add_executable(dabm_example d_abm.cpp) target_link_libraries(dabm_example PRIVATE memilio d_abm) target_compile_options(dabm_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) -# add_executable(smm_example smm.cpp) -# target_link_libraries(smm_example PRIVATE memilio smm) -# target_compile_options(smm_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +add_executable(smm_example smm.cpp) +target_link_libraries(smm_example PRIVATE memilio smm) +target_compile_options(smm_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) add_executable(temporal_hybrid_example temporal_hybrid_dabm_osecir.cpp) target_link_libraries(temporal_hybrid_example PRIVATE memilio hybrid ode_secir d_abm) diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index 5ee794e6d5..ca43ebda00 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -37,53 +37,49 @@ enum class InfectionState }; -enum class AgeGroup -{ - All, - Count -}; - int main() { //Example how to run the stochastic metapopulation models with four regions - const size_t num_regions = 4; - using Model = mio::smm::Model; + const size_t num_regions = 4; + const size_t num_age_groups = 1; + using Model = mio::smm::Model; double numE = 12, numC = 4, numI = 12, numR = 0, numD = 0; Model model; //Population are distributed uniformly to the four regions for (size_t r = 0; r < num_regions; ++r) { - model.populations[{mio::regions::Region(r), InfectionState::S, AgeGroup::All}] = + model.populations[{mio::regions::Region(r), InfectionState::S, mio::AgeGroup(1)}] = (1000 - numE - numC - numI - numR - numD) / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::E, AgeGroup::All}] = numE / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::C, AgeGroup::All}] = numC / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::I, AgeGroup::All}] = numI / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::R, AgeGroup::All}] = numR / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::D, AgeGroup::All}] = numD / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::E, mio::AgeGroup(1)}] = numE / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::C, mio::AgeGroup(1)}] = numC / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::I, mio::AgeGroup(1)}] = numI / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::R, mio::AgeGroup(1)}] = numR / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::D, mio::AgeGroup(1)}] = numD / num_regions; } //Set infection state adoption and spatial transition rates - std::vector> adoption_rates; - std::vector> transition_rates; + std::vector> adoption_rates; + std::vector> transition_rates; for (size_t r = 0; r < num_regions; ++r) { - adoption_rates.push_back({InfectionState::S, - InfectionState::E, - AgeGroup::All, - mio::regions::Region(r), - 0.1, - {{InfectionState::C, AgeGroup::All, 1}, {InfectionState::I, AgeGroup::All, 0.5}}}); adoption_rates.push_back( - {InfectionState::E, InfectionState::C, AgeGroup::All, mio::regions::Region(r), 1.0 / 5., {}}); + {InfectionState::S, + InfectionState::E, + mio::regions::Region(r), + 0.1, + {{InfectionState::C, 1, mio::AgeGroup(1)}, {InfectionState::I, 0.5, mio::AgeGroup(1)}}, + mio::AgeGroup(1)}); + adoption_rates.push_back( + {InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, mio::AgeGroup(1)}); adoption_rates.push_back( - {InfectionState::C, InfectionState::R, AgeGroup::All, mio::regions::Region(r), 0.2 / 3., {}}); + {InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, mio::AgeGroup(1)}); adoption_rates.push_back( - {InfectionState::C, InfectionState::I, AgeGroup::All, mio::regions::Region(r), 0.8 / 3., {}}); + {InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, mio::AgeGroup(1)}); adoption_rates.push_back( - {InfectionState::I, InfectionState::R, AgeGroup::All, mio::regions::Region(r), 0.99 / 5., {}}); + {InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, mio::AgeGroup(1)}); adoption_rates.push_back( - {InfectionState::I, InfectionState::D, AgeGroup::All, mio::regions::Region(r), 0.01 / 5., {}}); + {InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, mio::AgeGroup(1)}); } //Agents in infection state D do not transition @@ -91,16 +87,26 @@ int main() for (size_t i = 0; i < num_regions; ++i) { for (size_t j = 0; j < num_regions; ++j) if (i != j) { - transition_rates.push_back( - {InfectionState(s), AgeGroup::All, mio::regions::Region(i), mio::regions::Region(j), 0.01}); - transition_rates.push_back( - {InfectionState(s), AgeGroup::All, mio::regions::Region(j), mio::regions::Region(i), 0.01}); + transition_rates.push_back({ + InfectionState(s), + mio::regions::Region(i), + mio::regions::Region(j), + 0.01, + mio::AgeGroup(1), + }); + transition_rates.push_back({ + InfectionState(s), + mio::regions::Region(j), + mio::regions::Region(i), + 0.01, + mio::AgeGroup(1), + }); } } } - model.parameters.get>() = adoption_rates; - model.parameters.get>() = transition_rates; + model.parameters.get>() = adoption_rates; + model.parameters.get>() = transition_rates; double dt = 0.1; double tmax = 30.; diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h index b05ac2b0a8..56f1bb78f1 100644 --- a/cpp/memilio/epidemiology/adoption_rate.h +++ b/cpp/memilio/epidemiology/adoption_rate.h @@ -36,6 +36,7 @@ template struct Influence { Status status; ScalarType factor; + std::tuple group_indices{}; }; /** @@ -52,6 +53,7 @@ struct AdoptionRate { mio::regions::Region region; // k ScalarType factor; // gammahat_{ij}^k std::vector> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) + std::tuple group_indices{}; }; } // namespace mio diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index 4699fb7e27..ccd9d7ff38 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * -* Authors: René Schmieding, Julia Bicker +* Authors: René Schmieding, Julia Bicker, Kilian Volmer * * Contact: Martin J. Kuehn * @@ -22,6 +22,7 @@ #define MIO_SMM_MODEL_H #include "memilio/config.h" +#include "memilio/epidemiology/age_group.h" #include "smm/parameters.h" #include "memilio/compartments/compartmental_model.h" #include "memilio/epidemiology/populations.h" @@ -32,23 +33,29 @@ namespace mio namespace smm { +template +using age_group = T; + /** * @brief Stochastic Metapopulation Model. * @tparam regions Number of regions. * @tparam Status An infection state enum. */ -template -class Model : public mio::CompartmentalModel, - ParametersBase> +template +class Model : public mio::CompartmentalModel< + ScalarType, Status, + mio::Populations...>, + ParametersBase...>> { - using Base = mio::CompartmentalModel, - ParametersBase>; + using Base = mio::CompartmentalModel< + ScalarType, Status, + mio::Populations...>, + ParametersBase...>>; public: Model() - : Base(typename Base::Populations({static_cast(regions), Status::Count}), + : Base(typename Base::Populations( + {static_cast(regions), Status::Count, static_cast(Groups)...}), typename Base::ParameterSet()) { } @@ -59,26 +66,55 @@ class Model : public mio::CompartmentalModel& rate, const Eigen::VectorXd& x) const + ScalarType evaluate(const AdoptionRate...>& rate, + const Eigen::VectorXd& x) const { - const auto& pop = this->populations; - const auto source = pop.get_flat_index({rate.region, rate.from}); // Why is here rate.from used? KV - // determine order and calculate rate - if (rate.influences.size() == 0) { // first order adoption - return rate.factor * x[source]; + if constexpr (sizeof...(Groups) == 0) { + const auto& pop = this->populations; + const auto source = pop.get_flat_index({rate.region, rate.from}); // Why is here rate.from used? KV + // determine order and calculate rate + if (rate.influences.size() == 0) { // first order adoption + return rate.factor * x[source]; + } + else { // second order adoption + ScalarType N = 0; + for (size_t s = 0; s < static_cast(Status::Count); ++s) { + N += x[pop.get_flat_index({rate.region, Status(s)})]; + } + // accumulate influences + ScalarType influences = 0.0; + for (size_t i = 0; i < rate.influences.size(); i++) { + influences += + rate.influences[i].factor * x[pop.get_flat_index({rate.region, rate.influences[i].status})]; + } + return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; + } } - else { // second order adoption - ScalarType N = 0; - for (size_t s = 0; s < static_cast(Status::Count); ++s) { - N += x[pop.get_flat_index({rate.region, Status(s)})]; + else { + const auto& pop = this->populations; + const auto source = pop.get_flat_index({rate.region, rate.from, + std::make_from_tuple>( + rate.group_indices)}); // Why is here rate.from used? KV + // determine order and calculate rate + if (rate.influences.size() == 0) { // first order adoption + return rate.factor * x[source]; } - // accumulate influences - ScalarType influences = 0.0; - for (size_t i = 0; i < rate.influences.size(); i++) { - influences += - rate.influences[i].factor * x[pop.get_flat_index({rate.region, rate.influences[i].status})]; + else { // second order adoption + ScalarType N = 0; + for (size_t s = 0; s < static_cast(Status::Count); ++s) { + N += x[pop.get_flat_index( + {rate.region, Status(s), std::make_from_tuple>(rate.group_indices)})]; + } + // accumulate influences + ScalarType influences = 0.0; + for (size_t i = 0; i < rate.influences.size(); i++) { + influences += + rate.influences[i].factor * + x[pop.get_flat_index({rate.region, rate.influences[i].status, + std::make_from_tuple>(rate.group_indices)})]; + } + return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; } - return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; } } @@ -88,10 +124,18 @@ class Model : public mio::CompartmentalModel& rate, const Eigen::VectorXd& x) const + ScalarType evaluate(const TransitionRate...>& rate, + const Eigen::VectorXd& x) const { - const auto source = this->populations.get_flat_index({rate.from, rate.status}); - return rate.factor * x[source]; + if constexpr (sizeof...(Groups) == 0) { + const auto source = this->populations.get_flat_index({rate.from, rate.status}); + return rate.factor * x[source]; + } + else { + const auto source = this->populations.get_flat_index( + {rate.from, rate.status, std::make_from_tuple>(rate.group_indices)}); + return rate.factor * x[source]; + } } /** diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h index 8a251f63af..945cd08f12 100644 --- a/cpp/models/smm/parameters.h +++ b/cpp/models/smm/parameters.h @@ -25,6 +25,7 @@ #include "memilio/geography/regions.h" #include "memilio/utils/parameter_set.h" #include "memilio/epidemiology/adoption_rate.h" +#include namespace mio { @@ -54,6 +55,7 @@ struct TransitionRate { mio::regions::Region from; // k mio::regions::Region to; // l ScalarType factor; // lambda_i^{kl} + std::tuple group_indices{}; }; template struct TransitionRates { @@ -65,7 +67,7 @@ struct TransitionRates { }; template -using ParametersBase = mio::ParameterSet, TransitionRates>; +using ParametersBase = mio::ParameterSet, TransitionRates>; } // namespace smm diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 01afaa1689..7a7bd69590 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -22,9 +22,12 @@ #define MIO_SMM_SIMULATION_H #include "memilio/config.h" +#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/logging.h" #include "smm/model.h" #include "smm/parameters.h" #include "memilio/compartments/simulation.h" +#include namespace mio { @@ -37,7 +40,7 @@ namespace smm * @tparam regions The number of regions. * @tparam Status An infection state enum. */ -template +template class Simulation { public: @@ -79,56 +82,118 @@ class Simulation */ Eigen::Ref advance(ScalarType tmax) { - update_current_rates_and_waiting_times(); - size_t next_event = determine_next_event(); // index of the next event - ScalarType current_time = m_result.get_last_time(); - // set in the past to add a new time point immediately - ScalarType last_result_time = current_time - m_dt; - // iterate over time - while (current_time + m_waiting_times[next_event] < tmax) { - // update time - current_time += m_waiting_times[next_event]; - // regularily save current state in m_results - if (current_time > last_result_time + m_dt) { - last_result_time = current_time; - m_result.add_time_point(current_time); - // copy from the previous last value - m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; + if constexpr (sizeof...(Groups) == 0) { + update_current_rates_and_waiting_times(); + size_t next_event = determine_next_event(); // index of the next event + ScalarType current_time = m_result.get_last_time(); + // set in the past to add a new time point immediately + ScalarType last_result_time = current_time - m_dt; + // iterate over time + while (current_time + m_waiting_times[next_event] < tmax) { + // update time + current_time += m_waiting_times[next_event]; + // regularily save current state in m_results + if (current_time > last_result_time + m_dt) { + last_result_time = current_time; + m_result.add_time_point(current_time); + // copy from the previous last value + m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; + } + // decide event type by index and perform it + if (next_event < adoption_rates().size()) { + // perform adoption event + const auto& rate = adoption_rates()[next_event]; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from})] -= 1; + m_model->populations[{rate.region, rate.from}] -= 1; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to})] += 1; + m_model->populations[{rate.region, rate.to}] += 1; + } + else { + // perform transition event + const auto& rate = transition_rates()[next_event - adoption_rates().size()]; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status})] -= 1; + m_model->populations[{rate.from, rate.status}] -= 1; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status})] += 1; + m_model->populations[{rate.to, rate.status}] += 1; + } + // update internal times + for (size_t i = 0; i < m_internal_time.size(); i++) { + m_internal_time[i] += m_current_rates[i] * m_waiting_times[next_event]; + } + // draw new "next event" time for the occured event + m_tp_next_event[next_event] += + mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); + // precalculate next event + update_current_rates_and_waiting_times(); + next_event = determine_next_event(); } - // decide event type by index and perform it - if (next_event < adoption_rates().size()) { - // perform adoption event - const auto& rate = adoption_rates()[next_event]; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from})] -= 1; - m_model->populations[{rate.region, rate.from}] -= 1; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to})] += 1; - m_model->populations[{rate.region, rate.to}] += 1; + // copy last result, if no event occurs between last_result_time and tmax + if (last_result_time < tmax) { + m_result.add_time_point(tmax); + m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; } - else { - // perform transition event - const auto& rate = transition_rates()[next_event - adoption_rates().size()]; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status})] -= 1; - m_model->populations[{rate.from, rate.status}] -= 1; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status})] += 1; - m_model->populations[{rate.to, rate.status}] += 1; + return m_result.get_last_value(); + } + else { + update_current_rates_and_waiting_times(); + size_t next_event = determine_next_event(); // index of the next event + ScalarType current_time = m_result.get_last_time(); + // set in the past to add a new time point immediately + ScalarType last_result_time = current_time - m_dt; + // iterate over time + while (current_time + m_waiting_times[next_event] < tmax) { + // update time + current_time += m_waiting_times[next_event]; + // regularily save current state in m_results + if (current_time > last_result_time + m_dt) { + last_result_time = current_time; + m_result.add_time_point(current_time); + // copy from the previous last value + m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; + } + // decide event type by index and perform it + if (next_event < adoption_rates().size()) { + // perform adoption event + const auto& rate = adoption_rates()[next_event]; + m_result.get_last_value()[m_model->populations.get_flat_index( + {rate.region, rate.from, std::make_from_tuple(rate.group_indices)})] -= 1; + m_model->populations[{rate.region, rate.from, + std::make_from_tuple(rate.group_indices)}] -= 1; + m_result.get_last_value()[m_model->populations.get_flat_index( + {rate.region, rate.to, std::make_from_tuple(rate.group_indices)})] += 1; + m_model->populations[{rate.region, rate.to, + std::make_from_tuple(rate.group_indices)}] += 1; + } + else { + // perform transition event + const auto& rate = transition_rates()[next_event - adoption_rates().size()]; + m_result.get_last_value()[m_model->populations.get_flat_index( + {rate.from, rate.status, std::make_from_tuple(rate.group_indices)})] -= 1; + m_model->populations[{rate.from, rate.status, + std::make_from_tuple(rate.group_indices)}] -= 1; + m_result.get_last_value()[m_model->populations.get_flat_index( + {rate.to, rate.status, std::make_from_tuple(rate.group_indices)})] += 1; + m_model->populations[{rate.to, rate.status, + std::make_from_tuple(rate.group_indices)}] += 1; + } + // update internal times + for (size_t i = 0; i < m_internal_time.size(); i++) { + m_internal_time[i] += m_current_rates[i] * m_waiting_times[next_event]; + } + // draw new "next event" time for the occured event + m_tp_next_event[next_event] += + mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); + // precalculate next event + update_current_rates_and_waiting_times(); + next_event = determine_next_event(); } - // update internal times - for (size_t i = 0; i < m_internal_time.size(); i++) { - m_internal_time[i] += m_current_rates[i] * m_waiting_times[next_event]; + // copy last result, if no event occurs between last_result_time and tmax + if (last_result_time < tmax) { + m_result.add_time_point(tmax); + m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; } - // draw new "next event" time for the occured event - m_tp_next_event[next_event] += - mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); - // precalculate next event - update_current_rates_and_waiting_times(); - next_event = determine_next_event(); - } - // copy last result, if no event occurs between last_result_time and tmax - if (last_result_time < tmax) { - m_result.add_time_point(tmax); - m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; + return m_result.get_last_value(); } - return m_result.get_last_value(); } /** @@ -160,21 +225,24 @@ class Simulation /** * @brief Returns the model's transition rates. */ - inline constexpr const typename smm::TransitionRates::Type& transition_rates() + inline constexpr const typename smm::TransitionRates...>::Type& + transition_rates() { - return m_model->parameters.template get>(); + return m_model->parameters.template get...>>(); } /** * @brief Returns the model's adoption rates. */ - inline constexpr const typename smm::AdoptionRates::Type& adoption_rates() + inline constexpr const typename smm::AdoptionRates...>::Type& + adoption_rates() { - return m_model->parameters.template get>(); + return m_model->parameters.template get...>>(); } /** - * @brief Calculate current values for m_current_rates and m_waiting_times. + * @brief + * */ inline void update_current_rates_and_waiting_times() { From bf4b839195421dab3d290e06dbfc6ab9d76c2a21 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Tue, 9 Sep 2025 17:42:24 +0200 Subject: [PATCH 04/19] CHG: Compiling example with multiple age groups --- cpp/examples/smm.cpp | 45 +++++++++++++++++++----------------- cpp/models/smm/model.h | 38 ++++++++++++++++++++---------- cpp/models/smm/simulation.h | 46 ++++++++++++++++++++++++------------- 3 files changed, 80 insertions(+), 49 deletions(-) diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index ca43ebda00..6eca53be93 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -34,52 +34,55 @@ enum class InfectionState R, D, Count - }; +using Age = mio::AgeGroup; +using Species = mio::AgeGroup; + int main() { //Example how to run the stochastic metapopulation models with four regions const size_t num_regions = 4; const size_t num_age_groups = 1; - using Model = mio::smm::Model; + const size_t num_groups = 1; + using Model = mio::smm::Model; double numE = 12, numC = 4, numI = 12, numR = 0, numD = 0; Model model; //Population are distributed uniformly to the four regions for (size_t r = 0; r < num_regions; ++r) { - model.populations[{mio::regions::Region(r), InfectionState::S, mio::AgeGroup(1)}] = + model.populations[{mio::regions::Region(r), InfectionState::S, Age(1), Species(1)}] = (1000 - numE - numC - numI - numR - numD) / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::E, mio::AgeGroup(1)}] = numE / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::C, mio::AgeGroup(1)}] = numC / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::I, mio::AgeGroup(1)}] = numI / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::R, mio::AgeGroup(1)}] = numR / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::D, mio::AgeGroup(1)}] = numD / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::E, Age(1), Species(1)}] = numE / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::C, Age(1), Species(1)}] = numC / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::I, Age(1), Species(1)}] = numI / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::R, Age(1), Species(1)}] = numR / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::D, Age(1), Species(1)}] = numD / num_regions; } //Set infection state adoption and spatial transition rates - std::vector> adoption_rates; - std::vector> transition_rates; + std::vector> adoption_rates; + std::vector> transition_rates; for (size_t r = 0; r < num_regions; ++r) { adoption_rates.push_back( {InfectionState::S, InfectionState::E, mio::regions::Region(r), 0.1, - {{InfectionState::C, 1, mio::AgeGroup(1)}, {InfectionState::I, 0.5, mio::AgeGroup(1)}}, - mio::AgeGroup(1)}); + {{InfectionState::C, 1, {Age(1), Species(1)}}, {InfectionState::I, 0.5, {Age(1), Species(1)}}}, + {Age(1), Species(1)}}); adoption_rates.push_back( - {InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, mio::AgeGroup(1)}); + {InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, {Age(1), Species(1)}}); adoption_rates.push_back( - {InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, mio::AgeGroup(1)}); + {InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, {Age(1), Species(1)}}); adoption_rates.push_back( - {InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, mio::AgeGroup(1)}); + {InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, {Age(1), Species(1)}}); adoption_rates.push_back( - {InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, mio::AgeGroup(1)}); + {InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, {Age(1), Species(1)}}); adoption_rates.push_back( - {InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, mio::AgeGroup(1)}); + {InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, {Age(1), Species(1)}}); } //Agents in infection state D do not transition @@ -92,21 +95,21 @@ int main() mio::regions::Region(i), mio::regions::Region(j), 0.01, - mio::AgeGroup(1), + {Age(1), Species(1)}, }); transition_rates.push_back({ InfectionState(s), mio::regions::Region(j), mio::regions::Region(i), 0.01, - mio::AgeGroup(1), + {Age(1), Species(1)}, }); } } } - model.parameters.get>() = adoption_rates; - model.parameters.get>() = transition_rates; + model.parameters.get>() = adoption_rates; + model.parameters.get>() = transition_rates; double dt = 0.1; double tmax = 30.; diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index ccd9d7ff38..1966d76b28 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -51,6 +51,7 @@ class Model : public mio::CompartmentalModel< ScalarType, Status, mio::Populations...>, ParametersBase...>>; + using Index = mio::Index...>; public: Model() @@ -91,10 +92,13 @@ class Model : public mio::CompartmentalModel< } } else { - const auto& pop = this->populations; - const auto source = pop.get_flat_index({rate.region, rate.from, - std::make_from_tuple>( - rate.group_indices)}); // Why is here rate.from used? KV + const auto& pop = this->populations; + const auto index_from = std::apply( + [&](auto&&... args) { + return Index{rate.region, rate.from, std::forward(args)...}; + }, + rate.group_indices); + const auto source = pop.get_flat_index(index_from); // Why is here rate.from used? KV // determine order and calculate rate if (rate.influences.size() == 0) { // first order adoption return rate.factor * x[source]; @@ -102,16 +106,22 @@ class Model : public mio::CompartmentalModel< else { // second order adoption ScalarType N = 0; for (size_t s = 0; s < static_cast(Status::Count); ++s) { - N += x[pop.get_flat_index( - {rate.region, Status(s), std::make_from_tuple>(rate.group_indices)})]; + const auto index = std::apply( + [&](auto&&... args) { + return Index{rate.region, Status(s), std::forward(args)...}; + }, + rate.group_indices); + N += x[pop.get_flat_index(index)]; } // accumulate influences ScalarType influences = 0.0; for (size_t i = 0; i < rate.influences.size(); i++) { - influences += - rate.influences[i].factor * - x[pop.get_flat_index({rate.region, rate.influences[i].status, - std::make_from_tuple>(rate.group_indices)})]; + const auto index = std::apply( + [&](auto&&... args) { + return Index{rate.region, rate.influences[i].status, std::forward(args)...}; + }, + rate.group_indices); + influences += rate.influences[i].factor * x[pop.get_flat_index(index)]; } return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; } @@ -132,8 +142,12 @@ class Model : public mio::CompartmentalModel< return rate.factor * x[source]; } else { - const auto source = this->populations.get_flat_index( - {rate.from, rate.status, std::make_from_tuple>(rate.group_indices)}); + auto index = std::apply( + [&](auto&&... args) { + return Index{rate.from, rate.status, std::forward(args)...}; + }, + rate.group_indices); + const auto source = this->populations.get_flat_index(index); return rate.factor * x[source]; } } diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 7a7bd69590..177e83218b 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -23,6 +23,7 @@ #include "memilio/config.h" #include "memilio/epidemiology/age_group.h" +#include "memilio/utils/index.h" #include "memilio/utils/logging.h" #include "smm/model.h" #include "smm/parameters.h" @@ -45,6 +46,7 @@ class Simulation { public: using Model = smm::Model; + using Index = mio::Index...>; /** * @brief Set up the simulation for a Stochastic Metapopulation Model. @@ -155,26 +157,38 @@ class Simulation if (next_event < adoption_rates().size()) { // perform adoption event const auto& rate = adoption_rates()[next_event]; - m_result.get_last_value()[m_model->populations.get_flat_index( - {rate.region, rate.from, std::make_from_tuple(rate.group_indices)})] -= 1; - m_model->populations[{rate.region, rate.from, - std::make_from_tuple(rate.group_indices)}] -= 1; - m_result.get_last_value()[m_model->populations.get_flat_index( - {rate.region, rate.to, std::make_from_tuple(rate.group_indices)})] += 1; - m_model->populations[{rate.region, rate.to, - std::make_from_tuple(rate.group_indices)}] += 1; + auto index_from = std::apply( + [&](auto&&... args) { + return Index{rate.region, rate.from, std::forward(args)...}; + }, + rate.group_indices); + m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1; + m_model->populations[index_from] -= 1; + auto index_to = std::apply( + [&](auto&&... args) { + return Index{rate.region, rate.to, std::forward(args)...}; + }, + rate.group_indices); + m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1; + m_model->populations[index_to] += 1; } else { // perform transition event const auto& rate = transition_rates()[next_event - adoption_rates().size()]; - m_result.get_last_value()[m_model->populations.get_flat_index( - {rate.from, rate.status, std::make_from_tuple(rate.group_indices)})] -= 1; - m_model->populations[{rate.from, rate.status, - std::make_from_tuple(rate.group_indices)}] -= 1; - m_result.get_last_value()[m_model->populations.get_flat_index( - {rate.to, rate.status, std::make_from_tuple(rate.group_indices)})] += 1; - m_model->populations[{rate.to, rate.status, - std::make_from_tuple(rate.group_indices)}] += 1; + auto index_from = std::apply( + [&](auto&&... args) { + return Index{rate.from, rate.status, std::forward(args)...}; + }, + rate.group_indices); + m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1; + m_model->populations[index_from] -= 1; + auto index_to = std::apply( + [&](auto&&... args) { + return Index{rate.to, rate.status, std::forward(args)...}; + }, + rate.group_indices); + m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1; + m_model->populations[index_to] += 1; } // update internal times for (size_t i = 0; i < m_internal_time.size(); i++) { From cdba05668a28cb026e44833c5e133db57902ac75 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Tue, 9 Sep 2025 17:49:13 +0200 Subject: [PATCH 05/19] Compiling solution without code duplicates --- cpp/models/smm/model.h | 101 +++++++------------- cpp/models/smm/simulation.h | 185 +++++++++++++----------------------- 2 files changed, 102 insertions(+), 184 deletions(-) diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index 1966d76b28..637b2f2aba 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -70,61 +70,38 @@ class Model : public mio::CompartmentalModel< ScalarType evaluate(const AdoptionRate...>& rate, const Eigen::VectorXd& x) const { - if constexpr (sizeof...(Groups) == 0) { - const auto& pop = this->populations; - const auto source = pop.get_flat_index({rate.region, rate.from}); // Why is here rate.from used? KV - // determine order and calculate rate - if (rate.influences.size() == 0) { // first order adoption - return rate.factor * x[source]; - } - else { // second order adoption - ScalarType N = 0; - for (size_t s = 0; s < static_cast(Status::Count); ++s) { - N += x[pop.get_flat_index({rate.region, Status(s)})]; - } - // accumulate influences - ScalarType influences = 0.0; - for (size_t i = 0; i < rate.influences.size(); i++) { - influences += - rate.influences[i].factor * x[pop.get_flat_index({rate.region, rate.influences[i].status})]; - } - return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; - } + const auto& pop = this->populations; + const auto index_from = std::apply( + [&](auto&&... args) { + return Index{rate.region, rate.from, std::forward(args)...}; + }, + rate.group_indices); + const auto source = pop.get_flat_index(index_from); // Why is here rate.from used? KV + // determine order and calculate rate + if (rate.influences.size() == 0) { // first order adoption + return rate.factor * x[source]; } - else { - const auto& pop = this->populations; - const auto index_from = std::apply( - [&](auto&&... args) { - return Index{rate.region, rate.from, std::forward(args)...}; - }, - rate.group_indices); - const auto source = pop.get_flat_index(index_from); // Why is here rate.from used? KV - // determine order and calculate rate - if (rate.influences.size() == 0) { // first order adoption - return rate.factor * x[source]; + else { // second order adoption + ScalarType N = 0; + for (size_t s = 0; s < static_cast(Status::Count); ++s) { + const auto index = std::apply( + [&](auto&&... args) { + return Index{rate.region, Status(s), std::forward(args)...}; + }, + rate.group_indices); + N += x[pop.get_flat_index(index)]; } - else { // second order adoption - ScalarType N = 0; - for (size_t s = 0; s < static_cast(Status::Count); ++s) { - const auto index = std::apply( - [&](auto&&... args) { - return Index{rate.region, Status(s), std::forward(args)...}; - }, - rate.group_indices); - N += x[pop.get_flat_index(index)]; - } - // accumulate influences - ScalarType influences = 0.0; - for (size_t i = 0; i < rate.influences.size(); i++) { - const auto index = std::apply( - [&](auto&&... args) { - return Index{rate.region, rate.influences[i].status, std::forward(args)...}; - }, - rate.group_indices); - influences += rate.influences[i].factor * x[pop.get_flat_index(index)]; - } - return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; + // accumulate influences + ScalarType influences = 0.0; + for (size_t i = 0; i < rate.influences.size(); i++) { + const auto index = std::apply( + [&](auto&&... args) { + return Index{rate.region, rate.influences[i].status, std::forward(args)...}; + }, + rate.group_indices); + influences += rate.influences[i].factor * x[pop.get_flat_index(index)]; } + return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; } } @@ -137,19 +114,13 @@ class Model : public mio::CompartmentalModel< ScalarType evaluate(const TransitionRate...>& rate, const Eigen::VectorXd& x) const { - if constexpr (sizeof...(Groups) == 0) { - const auto source = this->populations.get_flat_index({rate.from, rate.status}); - return rate.factor * x[source]; - } - else { - auto index = std::apply( - [&](auto&&... args) { - return Index{rate.from, rate.status, std::forward(args)...}; - }, - rate.group_indices); - const auto source = this->populations.get_flat_index(index); - return rate.factor * x[source]; - } + auto index = std::apply( + [&](auto&&... args) { + return Index{rate.from, rate.status, std::forward(args)...}; + }, + rate.group_indices); + const auto source = this->populations.get_flat_index(index); + return rate.factor * x[source]; } /** diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 177e83218b..140afe070b 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -84,130 +84,77 @@ class Simulation */ Eigen::Ref advance(ScalarType tmax) { - if constexpr (sizeof...(Groups) == 0) { - update_current_rates_and_waiting_times(); - size_t next_event = determine_next_event(); // index of the next event - ScalarType current_time = m_result.get_last_time(); - // set in the past to add a new time point immediately - ScalarType last_result_time = current_time - m_dt; - // iterate over time - while (current_time + m_waiting_times[next_event] < tmax) { - // update time - current_time += m_waiting_times[next_event]; - // regularily save current state in m_results - if (current_time > last_result_time + m_dt) { - last_result_time = current_time; - m_result.add_time_point(current_time); - // copy from the previous last value - m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; - } - // decide event type by index and perform it - if (next_event < adoption_rates().size()) { - // perform adoption event - const auto& rate = adoption_rates()[next_event]; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from})] -= 1; - m_model->populations[{rate.region, rate.from}] -= 1; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to})] += 1; - m_model->populations[{rate.region, rate.to}] += 1; - } - else { - // perform transition event - const auto& rate = transition_rates()[next_event - adoption_rates().size()]; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status})] -= 1; - m_model->populations[{rate.from, rate.status}] -= 1; - m_result.get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status})] += 1; - m_model->populations[{rate.to, rate.status}] += 1; - } - // update internal times - for (size_t i = 0; i < m_internal_time.size(); i++) { - m_internal_time[i] += m_current_rates[i] * m_waiting_times[next_event]; - } - // draw new "next event" time for the occured event - m_tp_next_event[next_event] += - mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); - // precalculate next event - update_current_rates_and_waiting_times(); - next_event = determine_next_event(); - } - // copy last result, if no event occurs between last_result_time and tmax - if (last_result_time < tmax) { - m_result.add_time_point(tmax); + update_current_rates_and_waiting_times(); + size_t next_event = determine_next_event(); // index of the next event + ScalarType current_time = m_result.get_last_time(); + // set in the past to add a new time point immediately + ScalarType last_result_time = current_time - m_dt; + // iterate over time + while (current_time + m_waiting_times[next_event] < tmax) { + // update time + current_time += m_waiting_times[next_event]; + // regularily save current state in m_results + if (current_time > last_result_time + m_dt) { + last_result_time = current_time; + m_result.add_time_point(current_time); + // copy from the previous last value m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; } - return m_result.get_last_value(); - } - else { - update_current_rates_and_waiting_times(); - size_t next_event = determine_next_event(); // index of the next event - ScalarType current_time = m_result.get_last_time(); - // set in the past to add a new time point immediately - ScalarType last_result_time = current_time - m_dt; - // iterate over time - while (current_time + m_waiting_times[next_event] < tmax) { - // update time - current_time += m_waiting_times[next_event]; - // regularily save current state in m_results - if (current_time > last_result_time + m_dt) { - last_result_time = current_time; - m_result.add_time_point(current_time); - // copy from the previous last value - m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; - } - // decide event type by index and perform it - if (next_event < adoption_rates().size()) { - // perform adoption event - const auto& rate = adoption_rates()[next_event]; - auto index_from = std::apply( - [&](auto&&... args) { - return Index{rate.region, rate.from, std::forward(args)...}; - }, - rate.group_indices); - m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1; - m_model->populations[index_from] -= 1; - auto index_to = std::apply( - [&](auto&&... args) { - return Index{rate.region, rate.to, std::forward(args)...}; - }, - rate.group_indices); - m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1; - m_model->populations[index_to] += 1; - } - else { - // perform transition event - const auto& rate = transition_rates()[next_event - adoption_rates().size()]; - auto index_from = std::apply( - [&](auto&&... args) { - return Index{rate.from, rate.status, std::forward(args)...}; - }, - rate.group_indices); - m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1; - m_model->populations[index_from] -= 1; - auto index_to = std::apply( - [&](auto&&... args) { - return Index{rate.to, rate.status, std::forward(args)...}; - }, - rate.group_indices); - m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1; - m_model->populations[index_to] += 1; - } - // update internal times - for (size_t i = 0; i < m_internal_time.size(); i++) { - m_internal_time[i] += m_current_rates[i] * m_waiting_times[next_event]; - } - // draw new "next event" time for the occured event - m_tp_next_event[next_event] += - mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); - // precalculate next event - update_current_rates_and_waiting_times(); - next_event = determine_next_event(); + // decide event type by index and perform it + if (next_event < adoption_rates().size()) { + // perform adoption event + const auto& rate = adoption_rates()[next_event]; + auto index_from = std::apply( + [&](auto&&... args) { + return Index{rate.region, rate.from, std::forward(args)...}; + }, + rate.group_indices); + m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1; + m_model->populations[index_from] -= 1; + auto index_to = std::apply( + [&](auto&&... args) { + return Index{rate.region, rate.to, std::forward(args)...}; + }, + rate.group_indices); + m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1; + m_model->populations[index_to] += 1; } - // copy last result, if no event occurs between last_result_time and tmax - if (last_result_time < tmax) { - m_result.add_time_point(tmax); - m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; + else { + // perform transition event + const auto& rate = transition_rates()[next_event - adoption_rates().size()]; + auto index_from = std::apply( + [&](auto&&... args) { + return Index{rate.from, rate.status, std::forward(args)...}; + }, + rate.group_indices); + m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1; + m_model->populations[index_from] -= 1; + auto index_to = std::apply( + [&](auto&&... args) { + return Index{rate.to, rate.status, std::forward(args)...}; + }, + rate.group_indices); + m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1; + m_model->populations[index_to] += 1; + } + // update internal times + for (size_t i = 0; i < m_internal_time.size(); i++) { + m_internal_time[i] += m_current_rates[i] * m_waiting_times[next_event]; } - return m_result.get_last_value(); + // draw new "next event" time for the occured event + m_tp_next_event[next_event] += + mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); + // precalculate next event + update_current_rates_and_waiting_times(); + next_event = determine_next_event(); + } + // copy last result, if no event occurs between last_result_time and tmax + if (last_result_time < tmax) { + m_result.add_time_point(tmax); + m_result.get_last_value() = m_result[m_result.get_num_time_points() - 2]; } + return m_result.get_last_value(); + // } } /** From c0cb41b6be31e5fd5b51a1c56ed9a0a6b6e70f07 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Fri, 12 Sep 2025 20:09:54 +0200 Subject: [PATCH 06/19] CHG: semifunctional state --- cpp/examples/smm.cpp | 50 ++++++++++++------------ cpp/memilio/epidemiology/adoption_rate.h | 47 ++++++++++++++++++++++ cpp/models/smm/model.h | 2 +- cpp/models/smm/parameters.h | 3 +- cpp/models/smm/simulation.h | 4 +- 5 files changed, 76 insertions(+), 30 deletions(-) diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index 6eca53be93..77d72d7f7c 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -66,23 +66,23 @@ int main() std::vector> adoption_rates; std::vector> transition_rates; for (size_t r = 0; r < num_regions; ++r) { + adoption_rates.push_back({InfectionState::S, + InfectionState::E, + mio::regions::Region(r), + 0.1, + {{InfectionState::C, 1, mio::regions::Region(1), {Age(1), Species(1)}}, + {InfectionState::I, 0.5, mio::regions::Region(1), {Age(1), Species(1)}}}, + {Age(1), Species(1)}}); adoption_rates.push_back( - {InfectionState::S, - InfectionState::E, - mio::regions::Region(r), - 0.1, - {{InfectionState::C, 1, {Age(1), Species(1)}}, {InfectionState::I, 0.5, {Age(1), Species(1)}}}, - {Age(1), Species(1)}}); + {{InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, {Age(1), Species(1)}}}); adoption_rates.push_back( - {InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, {Age(1), Species(1)}}); + {{InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, {Age(1), Species(1)}}}); adoption_rates.push_back( - {InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, {Age(1), Species(1)}}); + {{InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, {Age(1), Species(1)}}}); adoption_rates.push_back( - {InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, {Age(1), Species(1)}}); + {{InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, {Age(1), Species(1)}}}); adoption_rates.push_back( - {InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, {Age(1), Species(1)}}); - adoption_rates.push_back( - {InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, {Age(1), Species(1)}}); + {{InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, {Age(1), Species(1)}}}); } //Agents in infection state D do not transition @@ -90,20 +90,18 @@ int main() for (size_t i = 0; i < num_regions; ++i) { for (size_t j = 0; j < num_regions; ++j) if (i != j) { - transition_rates.push_back({ - InfectionState(s), - mio::regions::Region(i), - mio::regions::Region(j), - 0.01, - {Age(1), Species(1)}, - }); - transition_rates.push_back({ - InfectionState(s), - mio::regions::Region(j), - mio::regions::Region(i), - 0.01, - {Age(1), Species(1)}, - }); + transition_rates.push_back({InfectionState(s), + mio::regions::Region(i), + mio::regions::Region(j), + 0.01, + {Age(1), Species(1)}, + {Age(1), Species(1)}}); + transition_rates.push_back({InfectionState(s), + mio::regions::Region(j), + mio::regions::Region(i), + 0.01, + {Age(1), Species(1)}, + {Age(1), Species(1)}}); } } } diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h index 56f1bb78f1..dfe659e922 100644 --- a/cpp/memilio/epidemiology/adoption_rate.h +++ b/cpp/memilio/epidemiology/adoption_rate.h @@ -31,11 +31,13 @@ namespace mio * @brief Struct defining an influence for a second-order adoption. * The population having "status" is multiplied with "factor." * @tparam Status An infection state enum. + * @tparam Groups Additional grouping indices. */ template struct Influence { Status status; ScalarType factor; + mio::regions::Region region; std::tuple group_indices{}; }; @@ -45,6 +47,7 @@ struct Influence { * In the d_abm and smm simulations, "from" is implicitly an influence, scaled by "factor". This is multiplied by * the sum over all "influences", which scale their "status" with the respective "factor". * @tparam Status An infection state enum. + * @tparam Groups Additional grouping indices. */ template struct AdoptionRate { @@ -54,6 +57,50 @@ struct AdoptionRate { ScalarType factor; // gammahat_{ij}^k std::vector> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) std::tuple group_indices{}; + + template + AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, + ScalarType initializer_factor, std::vector> second_order_influences, + T... Ts) + : from(initializer_from) + , to(initializer_to) + , region(initializer_region) + , factor(initializer_factor) + , group_indices{} + { + for (const auto& [status, scalar] : second_order_influences) { + influences.emplace_back(Influence{status, scalar, initializer_region}); + } + mio::unused(Ts...); + } + template + AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, + ScalarType initializer_factor, + std::vector>> + second_order_influences, + T... Ts) + : from(initializer_from) + , to(initializer_to) + , region(initializer_region) + , factor(initializer_factor) + , group_indices{} + { + for (const auto& [status, scalar, influence_region, groups] : second_order_influences) { + influences.emplace_back(Influence{status, scalar, influence_region, groups}); + } + mio::unused(Ts...); + } + AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, + ScalarType initializer_factor, std::tuple<> second_order_influences) + : from(initializer_from) + , to(initializer_to) + , region(initializer_region) + , factor(initializer_factor) + , influences{} + , group_indices{} + { + mio::unused(second_order_influences); + } }; } // namespace mio diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index 637b2f2aba..f2e1a8d904 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -118,7 +118,7 @@ class Model : public mio::CompartmentalModel< [&](auto&&... args) { return Index{rate.from, rate.status, std::forward(args)...}; }, - rate.group_indices); + rate.group_indices_from); const auto source = this->populations.get_flat_index(index); return rate.factor * x[source]; } diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h index 945cd08f12..71c1c90b94 100644 --- a/cpp/models/smm/parameters.h +++ b/cpp/models/smm/parameters.h @@ -55,7 +55,8 @@ struct TransitionRate { mio::regions::Region from; // k mio::regions::Region to; // l ScalarType factor; // lambda_i^{kl} - std::tuple group_indices{}; + std::tuple group_indices_from{}; + std::tuple group_indices_to{}; }; template struct TransitionRates { diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 140afe070b..924036736d 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -126,14 +126,14 @@ class Simulation [&](auto&&... args) { return Index{rate.from, rate.status, std::forward(args)...}; }, - rate.group_indices); + rate.group_indices_from); m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1; m_model->populations[index_from] -= 1; auto index_to = std::apply( [&](auto&&... args) { return Index{rate.to, rate.status, std::forward(args)...}; }, - rate.group_indices); + rate.group_indices_to); m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1; m_model->populations[index_to] += 1; } From ee5143b4963199c4e0a2faebe0db8a1779a1dff8 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Fri, 26 Sep 2025 13:35:22 +0200 Subject: [PATCH 07/19] CHG: Change numbers in rates --- cpp/examples/smm.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index 77d72d7f7c..ad331bee19 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -74,15 +74,15 @@ int main() {InfectionState::I, 0.5, mio::regions::Region(1), {Age(1), Species(1)}}}, {Age(1), Species(1)}}); adoption_rates.push_back( - {{InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, {Age(1), Species(1)}}}); + {InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, {Age(1), Species(1)}}); adoption_rates.push_back( - {{InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, {Age(1), Species(1)}}}); + {InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, {Age(1), Species(1)}}); adoption_rates.push_back( - {{InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, {Age(1), Species(1)}}}); + {InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, {Age(1), Species(1)}}); adoption_rates.push_back( - {{InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, {Age(1), Species(1)}}}); + {InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, {Age(1), Species(1)}}); adoption_rates.push_back( - {{InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, {Age(1), Species(1)}}}); + {InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, {Age(1), Species(1)}}); } //Agents in infection state D do not transition From 4823827079d993d9b3ec2c80f8a4d3b4d674100c Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Fri, 26 Sep 2025 13:35:40 +0200 Subject: [PATCH 08/19] CHG: Move Influence into AdoptionRate --- cpp/memilio/epidemiology/adoption_rate.h | 159 +++++++++++++++-------- 1 file changed, 106 insertions(+), 53 deletions(-) diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h index dfe659e922..4cf13f3970 100644 --- a/cpp/memilio/epidemiology/adoption_rate.h +++ b/cpp/memilio/epidemiology/adoption_rate.h @@ -31,15 +31,17 @@ namespace mio * @brief Struct defining an influence for a second-order adoption. * The population having "status" is multiplied with "factor." * @tparam Status An infection state enum. - * @tparam Groups Additional grouping indices. - */ -template -struct Influence { - Status status; - ScalarType factor; - mio::regions::Region region; - std::tuple group_indices{}; -}; +// * @tparam Groups Additional grouping indices. +// */ +// template +// struct Influence { +// Status status; +// ScalarType factor; +// mio::regions::Region region{0}; +// std::tuple group_indices{}; +// }; +//try out to define Influence struct in AdoptionRate struct +// otherwise Influence struct in Adoption Rate class /** * @brief Struct defining a possible status adoption in a Model based on Poisson Processes. @@ -51,56 +53,107 @@ struct Influence { */ template struct AdoptionRate { + struct Influence { + Status status; + ScalarType factor; + mio::regions::Region region = region; + std::tuple group_indices{}; + }; + Status from; // i Status to; // j mio::regions::Region region; // k ScalarType factor; // gammahat_{ij}^k - std::vector> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) + std::vector influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) std::tuple group_indices{}; - template - AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, - ScalarType initializer_factor, std::vector> second_order_influences, - T... Ts) - : from(initializer_from) - , to(initializer_to) - , region(initializer_region) - , factor(initializer_factor) - , group_indices{} - { - for (const auto& [status, scalar] : second_order_influences) { - influences.emplace_back(Influence{status, scalar, initializer_region}); - } - mio::unused(Ts...); - } - template - AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, - ScalarType initializer_factor, - std::vector>> - second_order_influences, - T... Ts) - : from(initializer_from) - , to(initializer_to) - , region(initializer_region) - , factor(initializer_factor) - , group_indices{} - { - for (const auto& [status, scalar, influence_region, groups] : second_order_influences) { - influences.emplace_back(Influence{status, scalar, influence_region, groups}); - } - mio::unused(Ts...); - } - AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, - ScalarType initializer_factor, std::tuple<> second_order_influences) - : from(initializer_from) - , to(initializer_to) - , region(initializer_region) - , factor(initializer_factor) - , influences{} - , group_indices{} - { - mio::unused(second_order_influences); - } + // template + // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, + // ScalarType initializer_factor, std::vector> second_order_influences, + // T... Ts) + // : from(initializer_from) + // , to(initializer_to) + // , region(initializer_region) + // , factor(initializer_factor) + // , group_indices{} + // { + // for (const auto& [status, scalar] : second_order_influences) { + // influences.emplace_back(Influence{status, scalar, initializer_region}); + // } + // mio::unused(Ts...); + // } + // template + // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, + // ScalarType initializer_factor, + // std::vector>> + // second_order_influences, + // T... Ts) + // : from(initializer_from) + // , to(initializer_to) + // , region(initializer_region) + // , factor(initializer_factor) + // , group_indices{} + // { + // for (const auto& [status, scalar, influence_region, groups] : second_order_influences) { + // influences.emplace_back(Influence{status, scalar, influence_region, groups}); + // } + // mio::unused(Ts...); + // } + // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, + // ScalarType initializer_factor, std::tuple<> second_order_influences) + // : from(initializer_from) + // , to(initializer_to) + // , region(initializer_region) + // , factor(initializer_factor) + // , influences{} + // , group_indices{} + // { + // mio::unused(second_order_influences); + // } + // template + // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, + // ScalarType initializer_factor, std::vector> second_order_influences, + // std::tuple initializer_groups, T... Ts) + // : from(initializer_from) + // , to(initializer_to) + // , region(initializer_region) + // , factor(initializer_factor) + // , group_indices{initializer_groups} + // { + // for (const auto& [status, scalar] : second_order_influences) { + // influences.emplace_back(Influence{status, scalar, initializer_region}); + // } + // mio::unused(Ts...); + // } + // template + // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, + // ScalarType initializer_factor, + // std::vector>> + // second_order_influences, + // std::tuple initializer_groups, T... Ts) + // : from(initializer_from) + // , to(initializer_to) + // , region(initializer_region) + // , factor(initializer_factor) + // , group_indices{initializer_groups} + // { + // for (const auto& [status, scalar, influence_region, groups] : second_order_influences) { + // influences.emplace_back(Influence{status, scalar, influence_region, groups}); + // } + // mio::unused(Ts...); + // } + // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, + // ScalarType initializer_factor, std::tuple<> second_order_influences, + // std::tuple initializer_groups) + // : from(initializer_from) + // , to(initializer_to) + // , region(initializer_region) + // , factor(initializer_factor) + // , influences{} + // , group_indices{initializer_groups} + // { + // mio::unused(second_order_influences); + // } }; } // namespace mio From 893214413513219c59a974f81baf8d82bc148e58 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Wed, 1 Oct 2025 13:18:47 +0200 Subject: [PATCH 09/19] CHG: Update influences --- cpp/examples/smm.cpp | 36 +++---- cpp/memilio/epidemiology/adoption_rate.h | 119 +++-------------------- cpp/models/smm/model.h | 5 +- 3 files changed, 34 insertions(+), 126 deletions(-) diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index ad331bee19..cb9c0f62fc 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -53,13 +53,13 @@ int main() Model model; //Population are distributed uniformly to the four regions for (size_t r = 0; r < num_regions; ++r) { - model.populations[{mio::regions::Region(r), InfectionState::S, Age(1), Species(1)}] = + model.populations[{mio::regions::Region(r), InfectionState::S, Age(0), Species(0)}] = (1000 - numE - numC - numI - numR - numD) / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::E, Age(1), Species(1)}] = numE / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::C, Age(1), Species(1)}] = numC / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::I, Age(1), Species(1)}] = numI / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::R, Age(1), Species(1)}] = numR / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::D, Age(1), Species(1)}] = numD / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::E, Age(0), Species(0)}] = numE / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::C, Age(0), Species(0)}] = numC / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::I, Age(0), Species(0)}] = numI / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::R, Age(0), Species(0)}] = numR / num_regions; + model.populations[{mio::regions::Region(r), InfectionState::D, Age(0), Species(0)}] = numD / num_regions; } //Set infection state adoption and spatial transition rates @@ -70,19 +70,19 @@ int main() InfectionState::E, mio::regions::Region(r), 0.1, - {{InfectionState::C, 1, mio::regions::Region(1), {Age(1), Species(1)}}, - {InfectionState::I, 0.5, mio::regions::Region(1), {Age(1), Species(1)}}}, - {Age(1), Species(1)}}); + {{InfectionState::C, 1, mio::regions::Region(3), {Age(0), Species(0)}}, + {InfectionState::I, 0.5, mio::regions::Region(1), {Age(0), Species(0)}}}, + {Age(0), Species(0)}}); adoption_rates.push_back( - {InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, {Age(1), Species(1)}}); + {InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, {Age(0), Species(0)}}); adoption_rates.push_back( - {InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, {Age(1), Species(1)}}); + {InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, {Age(0), Species(0)}}); adoption_rates.push_back( - {InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, {Age(1), Species(1)}}); + {InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, {Age(0), Species(0)}}); adoption_rates.push_back( - {InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, {Age(1), Species(1)}}); + {InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, {Age(0), Species(0)}}); adoption_rates.push_back( - {InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, {Age(1), Species(1)}}); + {InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, {Age(0), Species(0)}}); } //Agents in infection state D do not transition @@ -94,14 +94,14 @@ int main() mio::regions::Region(i), mio::regions::Region(j), 0.01, - {Age(1), Species(1)}, - {Age(1), Species(1)}}); + {Age(0), Species(0)}, + {Age(0), Species(0)}}); transition_rates.push_back({InfectionState(s), mio::regions::Region(j), mio::regions::Region(i), 0.01, - {Age(1), Species(1)}, - {Age(1), Species(1)}}); + {Age(0), Species(0)}, + {Age(0), Species(0)}}); } } } diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h index 4cf13f3970..2917d841d3 100644 --- a/cpp/memilio/epidemiology/adoption_rate.h +++ b/cpp/memilio/epidemiology/adoption_rate.h @@ -23,26 +23,12 @@ #include "memilio/utils/index.h" #include "memilio/config.h" #include "memilio/geography/regions.h" - +#include +#include +#include namespace mio { -/** - * @brief Struct defining an influence for a second-order adoption. - * The population having "status" is multiplied with "factor." - * @tparam Status An infection state enum. -// * @tparam Groups Additional grouping indices. -// */ -// template -// struct Influence { -// Status status; -// ScalarType factor; -// mio::regions::Region region{0}; -// std::tuple group_indices{}; -// }; -//try out to define Influence struct in AdoptionRate struct -// otherwise Influence struct in Adoption Rate class - /** * @brief Struct defining a possible status adoption in a Model based on Poisson Processes. * The AdoptionRate is considered to be of second-order if there are any "influences". @@ -53,10 +39,19 @@ namespace mio */ template struct AdoptionRate { + + /** + * @brief Struct defining an influence for a second-order adoption. + * The population having "status" is multiplied with "factor." + * @tparam status An infection state enum. + * @tparam factor Scaling factor for the influence. + * @tparam + * @tparam Groups Additional grouping indices. + */ struct Influence { Status status; ScalarType factor; - mio::regions::Region region = region; + std::optional region = std::nullopt; std::tuple group_indices{}; }; @@ -66,94 +61,6 @@ struct AdoptionRate { ScalarType factor; // gammahat_{ij}^k std::vector influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) std::tuple group_indices{}; - - // template - // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, - // ScalarType initializer_factor, std::vector> second_order_influences, - // T... Ts) - // : from(initializer_from) - // , to(initializer_to) - // , region(initializer_region) - // , factor(initializer_factor) - // , group_indices{} - // { - // for (const auto& [status, scalar] : second_order_influences) { - // influences.emplace_back(Influence{status, scalar, initializer_region}); - // } - // mio::unused(Ts...); - // } - // template - // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, - // ScalarType initializer_factor, - // std::vector>> - // second_order_influences, - // T... Ts) - // : from(initializer_from) - // , to(initializer_to) - // , region(initializer_region) - // , factor(initializer_factor) - // , group_indices{} - // { - // for (const auto& [status, scalar, influence_region, groups] : second_order_influences) { - // influences.emplace_back(Influence{status, scalar, influence_region, groups}); - // } - // mio::unused(Ts...); - // } - // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, - // ScalarType initializer_factor, std::tuple<> second_order_influences) - // : from(initializer_from) - // , to(initializer_to) - // , region(initializer_region) - // , factor(initializer_factor) - // , influences{} - // , group_indices{} - // { - // mio::unused(second_order_influences); - // } - // template - // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, - // ScalarType initializer_factor, std::vector> second_order_influences, - // std::tuple initializer_groups, T... Ts) - // : from(initializer_from) - // , to(initializer_to) - // , region(initializer_region) - // , factor(initializer_factor) - // , group_indices{initializer_groups} - // { - // for (const auto& [status, scalar] : second_order_influences) { - // influences.emplace_back(Influence{status, scalar, initializer_region}); - // } - // mio::unused(Ts...); - // } - // template - // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, - // ScalarType initializer_factor, - // std::vector>> - // second_order_influences, - // std::tuple initializer_groups, T... Ts) - // : from(initializer_from) - // , to(initializer_to) - // , region(initializer_region) - // , factor(initializer_factor) - // , group_indices{initializer_groups} - // { - // for (const auto& [status, scalar, influence_region, groups] : second_order_influences) { - // influences.emplace_back(Influence{status, scalar, influence_region, groups}); - // } - // mio::unused(Ts...); - // } - // AdoptionRate(Status initializer_from, Status initializer_to, mio::regions::Region initializer_region, - // ScalarType initializer_factor, std::tuple<> second_order_influences, - // std::tuple initializer_groups) - // : from(initializer_from) - // , to(initializer_to) - // , region(initializer_region) - // , factor(initializer_factor) - // , influences{} - // , group_indices{initializer_groups} - // { - // mio::unused(second_order_influences); - // } }; } // namespace mio diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index f2e1a8d904..878066552c 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -96,9 +96,10 @@ class Model : public mio::CompartmentalModel< for (size_t i = 0; i < rate.influences.size(); i++) { const auto index = std::apply( [&](auto&&... args) { - return Index{rate.region, rate.influences[i].status, std::forward(args)...}; + return Index{rate.influences[i].region.value_or(rate.region), rate.influences[i].status, + std::forward(args)...}; }, - rate.group_indices); + rate.influences[i].group_indices); influences += rate.influences[i].factor * x[pop.get_flat_index(index)]; } return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; From 098e1ff182c500d2e7fc6afb505b114cadd9e86c Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Wed, 1 Oct 2025 15:42:47 +0200 Subject: [PATCH 10/19] CHG: Change location of N, modify tests --- cpp/models/smm/model.h | 26 +++++----- cpp/tests/test_smm_model.cpp | 93 +++++++++++++++++++++++++++++------- 2 files changed, 91 insertions(+), 28 deletions(-) diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index 878066552c..1c3e3f8123 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -82,27 +82,31 @@ class Model : public mio::CompartmentalModel< return rate.factor * x[source]; } else { // second order adoption - ScalarType N = 0; - for (size_t s = 0; s < static_cast(Status::Count); ++s) { - const auto index = std::apply( - [&](auto&&... args) { - return Index{rate.region, Status(s), std::forward(args)...}; - }, - rate.group_indices); - N += x[pop.get_flat_index(index)]; - } // accumulate influences ScalarType influences = 0.0; for (size_t i = 0; i < rate.influences.size(); i++) { + ScalarType N = 0; // Welches N brauchen wir hier?? + + for (size_t s = 0; s < static_cast(Status::Count); ++s) { + const auto index = std::apply( + [&](auto&&... args) { + return Index{rate.influences[i].region.value_or(rate.region), Status(s), + std::forward(args)...}; + }, + rate.influences[i].group_indices); + N += x[pop.get_flat_index(index)]; + } const auto index = std::apply( [&](auto&&... args) { return Index{rate.influences[i].region.value_or(rate.region), rate.influences[i].status, std::forward(args)...}; }, rate.influences[i].group_indices); - influences += rate.influences[i].factor * x[pop.get_flat_index(index)]; + if (N > 0) { + influences += rate.influences[i].factor * x[pop.get_flat_index(index)] / N; + } } - return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; + return rate.factor * x[source] * influences; } } diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp index 39fe57451e..85a4a4b85c 100644 --- a/cpp/tests/test_smm_model.cpp +++ b/cpp/tests/test_smm_model.cpp @@ -23,7 +23,9 @@ #include #include #include +#include "memilio/epidemiology/age_group.h" #include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/logging.h" #include "smm/model.h" #include "smm/parameters.h" #include "smm/simulation.h" @@ -45,6 +47,14 @@ enum class InfectionState }; +enum class Infections +{ + S, + I, + R, + Count +}; + TEST(TestSMM, evaluateAdoptionRate) { //Test whether the adoption rates are evaluated correctly. @@ -80,31 +90,80 @@ TEST(TestSMM, evaluateAdoptionRate) EXPECT_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 2.); } +TEST(TestSMM, evaluateinterregionalAdoptionRate) +{ + using Model = mio::smm::Model<2, Infections, 2>; + using Age = mio::AgeGroup; + + std::vector> adoption_rates; + //Second-order adoption + adoption_rates.push_back({Infections::S, + Infections::I, + mio::regions::Region(0), + 0.1, + {{Infections::I, 0.1, mio::regions::Region(0), {Age(1)}}, + {Infections::I, 0.2, mio::regions::Region(1), {Age(1)}}}, + Age(0)}); + //First-order adoption + adoption_rates.push_back({Infections::I, Infections::R, mio::regions::Region(0), 0.2, {}, {Age(1)}}); + Model model; + + model.populations[{mio::regions::Region(0), Infections::S, Age(0)}] = 50; + model.populations[{mio::regions::Region(0), Infections::I, Age(0)}] = 10; + model.populations[{mio::regions::Region(0), Infections::R, Age(0)}] = 0; + + model.populations[{mio::regions::Region(0), Infections::S, Age(1)}] = 100; + model.populations[{mio::regions::Region(0), Infections::I, Age(1)}] = 20; + model.populations[{mio::regions::Region(0), Infections::R, Age(1)}] = 0; + + model.populations[{mio::regions::Region(1), Infections::S, Age(0)}] = 40; + model.populations[{mio::regions::Region(1), Infections::I, Age(0)}] = 80; + model.populations[{mio::regions::Region(1), Infections::R, Age(0)}] = 0; + + model.populations[{mio::regions::Region(1), Infections::S, Age(1)}] = 80; + model.populations[{mio::regions::Region(1), Infections::I, Age(1)}] = 16; + model.populations[{mio::regions::Region(1), Infections::R, Age(1)}] = 0; + + EXPECT_FLOAT_EQ(model.evaluate(adoption_rates[0], model.populations.get_compartments()), + 0.1 * 50. * (0.1 * 20. * 1. / 120. + 0.2 * 16 * 1. / 96.)); + EXPECT_FLOAT_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 4.); +} + TEST(TestSMM, evaluateTransitionRate) { //Same test as 'evaluateAdoptionRate' only for spatial transition rates. //Transition rates are given by: rate.factor * N(rate.status, rate.from) - using Model = mio::smm::Model<2, InfectionState>; + using Model = mio::smm::Model<2, InfectionState, 2>; Model model; //Initialize model populations - model.populations[{mio::regions::Region(0), InfectionState::S}] = 50; - model.populations[{mio::regions::Region(0), InfectionState::E}] = 10; - model.populations[{mio::regions::Region(0), InfectionState::C}] = 5; - model.populations[{mio::regions::Region(0), InfectionState::I}] = 0; - model.populations[{mio::regions::Region(0), InfectionState::R}] = 0; - model.populations[{mio::regions::Region(0), InfectionState::D}] = 0; - - model.populations[{mio::regions::Region(1), InfectionState::S}] = 55; - model.populations[{mio::regions::Region(1), InfectionState::E}] = 10; - model.populations[{mio::regions::Region(1), InfectionState::C}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::I}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::R}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::D}] = 0; + model.populations[{mio::regions::Region(0), InfectionState::S, mio::AgeGroup(0)}] = 50; + model.populations[{mio::regions::Region(0), InfectionState::E, mio::AgeGroup(0)}] = 10; + model.populations[{mio::regions::Region(0), InfectionState::C, mio::AgeGroup(0)}] = 5; + model.populations[{mio::regions::Region(0), InfectionState::I, mio::AgeGroup(0)}] = 0; + model.populations[{mio::regions::Region(0), InfectionState::R, mio::AgeGroup(0)}] = 0; + model.populations[{mio::regions::Region(0), InfectionState::D, mio::AgeGroup(0)}] = 0; + + model.populations[{mio::regions::Region(1), InfectionState::S, mio::AgeGroup(0)}] = 55; + model.populations[{mio::regions::Region(1), InfectionState::E, mio::AgeGroup(0)}] = 10; + model.populations[{mio::regions::Region(1), InfectionState::C, mio::AgeGroup(0)}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::I, mio::AgeGroup(0)}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::R, mio::AgeGroup(0)}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::D, mio::AgeGroup(0)}] = 0; //Set transition rates - std::vector> transition_rates; - transition_rates.push_back({InfectionState::S, mio::regions::Region(0), mio::regions::Region(1), 0.01}); - transition_rates.push_back({InfectionState::E, mio::regions::Region(1), mio::regions::Region(0), 0.1}); + std::vector> transition_rates; + transition_rates.push_back({InfectionState::S, + mio::regions::Region(0), + mio::regions::Region(1), + 0.01, + {mio::AgeGroup(0)}, + {mio::AgeGroup(0)}}); + transition_rates.push_back({InfectionState::E, + mio::regions::Region(1), + mio::regions::Region(0), + 0.1, + {mio::AgeGroup(0)}, + {mio::AgeGroup(0)}}); EXPECT_EQ(model.evaluate(transition_rates[0], model.populations.get_compartments()), 0.5); EXPECT_EQ(model.evaluate(transition_rates[1], model.populations.get_compartments()), 1.); From dcda987cdb2d5728ad65675893d1cf0110af4389 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Wed, 1 Oct 2025 15:48:42 +0200 Subject: [PATCH 11/19] [skip ci] CHG: Update Header --- cpp/examples/smm.cpp | 2 +- cpp/memilio/epidemiology/adoption_rate.h | 3 ++- cpp/tests/test_smm_model.cpp | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index cb9c0f62fc..cf457e537e 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) * -* Authors: Julia Bicker, René Schmieding +* Authors: Julia Bicker, René Schmieding, Kilian Volmer * * Contact: Martin J. Kuehn * diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h index 2917d841d3..86addd7544 100644 --- a/cpp/memilio/epidemiology/adoption_rate.h +++ b/cpp/memilio/epidemiology/adoption_rate.h @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2025 MEmilio * -* Authors: René Schmieding, Julia Bicker +* Authors: René Schmieding, Julia Bicker, Kilian Volmer * * Contact: Martin J. Kuehn * @@ -26,6 +26,7 @@ #include #include #include + namespace mio { diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp index 85a4a4b85c..555d660cc0 100644 --- a/cpp/tests/test_smm_model.cpp +++ b/cpp/tests/test_smm_model.cpp @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) * -* Authors: Julia Bicker +* Authors: Julia Bicker, Kilian Volmer * * Contact: Martin J. Kuehn * From c58b74ae408e89e7d501b4d90bb11fc4163537b6 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Thu, 2 Oct 2025 09:20:06 +0200 Subject: [PATCH 12/19] [skip ci] Update copyright information --- cpp/examples/d_abm.cpp | 2 +- cpp/examples/graph_abm.cpp | 2 +- cpp/examples/smm.cpp | 2 +- cpp/models/d_abm/model.cpp | 2 +- cpp/models/d_abm/model.h | 2 +- cpp/models/d_abm/parameters.h | 2 +- cpp/models/d_abm/quad_well.h | 2 +- cpp/models/d_abm/simulation.cpp | 2 +- cpp/models/d_abm/simulation.h | 2 +- cpp/models/d_abm/single_well.h | 2 +- cpp/models/graph_abm/graph_abm_mobility.cpp | 2 +- cpp/models/graph_abm/graph_abm_mobility.h | 2 +- cpp/models/graph_abm/graph_abmodel.h | 2 +- cpp/models/hybrid/temporal_hybrid_model.cpp | 2 +- cpp/models/hybrid/temporal_hybrid_model.h | 2 +- cpp/models/smm/model.cpp | 2 +- cpp/models/smm/model.h | 2 +- cpp/models/smm/parameters.h | 2 +- cpp/models/smm/simulation.h | 2 +- cpp/tests/test_d_abm_model.cpp | 2 +- cpp/tests/test_smm_model.cpp | 2 +- pycode/memilio-epidata/memilio/epidata/getJHData.py | 2 +- .../memilio/generation/template/template_cpp.txt | 2 +- .../memilio/generation/template/template_py.txt | 2 +- .../memilio/generation_test/test_data/test_oseir.cpp.txt | 2 +- .../memilio/generation_test/test_data/test_oseir.py.txt | 2 +- 26 files changed, 26 insertions(+), 26 deletions(-) diff --git a/cpp/examples/d_abm.cpp b/cpp/examples/d_abm.cpp index 09ffa04d5c..4f91f65de0 100644 --- a/cpp/examples/d_abm.cpp +++ b/cpp/examples/d_abm.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker, René Schmieding * diff --git a/cpp/examples/graph_abm.cpp b/cpp/examples/graph_abm.cpp index 3e7443fbc9..dedc28b60e 100644 --- a/cpp/examples/graph_abm.cpp +++ b/cpp/examples/graph_abm.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 MEmilio +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker * diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index 9eae9129a4..cf13a0c7be 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker, René Schmieding, Kilian Volmer * diff --git a/cpp/models/d_abm/model.cpp b/cpp/models/d_abm/model.cpp index f52bd1c298..528fa1a4f4 100644 --- a/cpp/models/d_abm/model.cpp +++ b/cpp/models/d_abm/model.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: René Schmieding * diff --git a/cpp/models/d_abm/model.h b/cpp/models/d_abm/model.h index b0983c5496..5b858a4efa 100644 --- a/cpp/models/d_abm/model.h +++ b/cpp/models/d_abm/model.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/d_abm/parameters.h b/cpp/models/d_abm/parameters.h index bd02befd9b..5b16eb433f 100644 --- a/cpp/models/d_abm/parameters.h +++ b/cpp/models/d_abm/parameters.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/d_abm/quad_well.h b/cpp/models/d_abm/quad_well.h index 305617304e..ef1d0e021e 100644 --- a/cpp/models/d_abm/quad_well.h +++ b/cpp/models/d_abm/quad_well.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/d_abm/simulation.cpp b/cpp/models/d_abm/simulation.cpp index 3b538490b8..3915b84300 100644 --- a/cpp/models/d_abm/simulation.cpp +++ b/cpp/models/d_abm/simulation.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 MEmilio +* Copyright (C) 2020-2025 MEmilio * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/d_abm/simulation.h b/cpp/models/d_abm/simulation.h index 54ccad8b65..e20f40f595 100644 --- a/cpp/models/d_abm/simulation.h +++ b/cpp/models/d_abm/simulation.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/d_abm/single_well.h b/cpp/models/d_abm/single_well.h index 62ef7671c6..7bd8cb82e1 100644 --- a/cpp/models/d_abm/single_well.h +++ b/cpp/models/d_abm/single_well.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker, René Schmieding * diff --git a/cpp/models/graph_abm/graph_abm_mobility.cpp b/cpp/models/graph_abm/graph_abm_mobility.cpp index af3fae0c89..c8482f2e48 100644 --- a/cpp/models/graph_abm/graph_abm_mobility.cpp +++ b/cpp/models/graph_abm/graph_abm_mobility.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 MEmilio +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker * diff --git a/cpp/models/graph_abm/graph_abm_mobility.h b/cpp/models/graph_abm/graph_abm_mobility.h index 2d7a749edf..1af88d311a 100644 --- a/cpp/models/graph_abm/graph_abm_mobility.h +++ b/cpp/models/graph_abm/graph_abm_mobility.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 MEmilio +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker * diff --git a/cpp/models/graph_abm/graph_abmodel.h b/cpp/models/graph_abm/graph_abmodel.h index 0be1a7f222..040cd8c2da 100644 --- a/cpp/models/graph_abm/graph_abmodel.h +++ b/cpp/models/graph_abm/graph_abmodel.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 MEmilio +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker * diff --git a/cpp/models/hybrid/temporal_hybrid_model.cpp b/cpp/models/hybrid/temporal_hybrid_model.cpp index f356bd2321..337b8da3de 100644 --- a/cpp/models/hybrid/temporal_hybrid_model.cpp +++ b/cpp/models/hybrid/temporal_hybrid_model.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker, René Schmieding * diff --git a/cpp/models/hybrid/temporal_hybrid_model.h b/cpp/models/hybrid/temporal_hybrid_model.h index 25ea35d8d0..3ba096c115 100644 --- a/cpp/models/hybrid/temporal_hybrid_model.h +++ b/cpp/models/hybrid/temporal_hybrid_model.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker, René Schmieding * diff --git a/cpp/models/smm/model.cpp b/cpp/models/smm/model.cpp index af0d58eb73..d613003ac2 100644 --- a/cpp/models/smm/model.cpp +++ b/cpp/models/smm/model.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: René Schmieding * diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index 7d6583bd08..b5942bc251 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: René Schmieding, Julia Bicker, Kilian Volmer * diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h index c7759d8953..fe27c35526 100644 --- a/cpp/models/smm/parameters.h +++ b/cpp/models/smm/parameters.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 8364236f28..8fb1423571 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/tests/test_d_abm_model.cpp b/cpp/tests/test_d_abm_model.cpp index c31b4aa5e4..79397d8e65 100644 --- a/cpp/tests/test_d_abm_model.cpp +++ b/cpp/tests/test_d_abm_model.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker * diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp index 3d8e07f2db..2351e33ef6 100644 --- a/cpp/tests/test_smm_model.cpp +++ b/cpp/tests/test_smm_model.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker, Kilian Volmer * diff --git a/pycode/memilio-epidata/memilio/epidata/getJHData.py b/pycode/memilio-epidata/memilio/epidata/getJHData.py index 19ec2b67f9..dd0a3bb873 100644 --- a/pycode/memilio-epidata/memilio/epidata/getJHData.py +++ b/pycode/memilio-epidata/memilio/epidata/getJHData.py @@ -1,5 +1,5 @@ ############################################################################# -# Copyright (C) 2020-2024 MEmilio +# Copyright (C) 2020-2025 MEmilio # # Authors: Kathrin Rack # diff --git a/pycode/memilio-generation/memilio/generation/template/template_cpp.txt b/pycode/memilio-generation/memilio/generation/template/template_cpp.txt index 28d9f80042..e33668d1ee 100644 --- a/pycode/memilio-generation/memilio/generation/template/template_cpp.txt +++ b/pycode/memilio-generation/memilio/generation/template/template_cpp.txt @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 MEmilio +* Copyright (C) 2020-2025 MEmilio * * Authors: Martin Siggel, Daniel Abele, Martin J. Kuehn, Jan Kleinert * diff --git a/pycode/memilio-generation/memilio/generation/template/template_py.txt b/pycode/memilio-generation/memilio/generation/template/template_py.txt index 2ae2f78e5a..d1cc9710c7 100644 --- a/pycode/memilio-generation/memilio/generation/template/template_py.txt +++ b/pycode/memilio-generation/memilio/generation/template/template_py.txt @@ -1,5 +1,5 @@ ############################################################################# -# Copyright (C) 2020-2024 MEmilio +# Copyright (C) 2020-2025 MEmilio # # Authors: Daniel Abele # diff --git a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt index b419020ba3..edacafb623 100644 --- a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt +++ b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 MEmilio +* Copyright (C) 2020-2025 MEmilio * * Authors: Martin Siggel, Daniel Abele, Martin J. Kuehn, Jan Kleinert * diff --git a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt index 3f289a1251..1fd1861e46 100644 --- a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt +++ b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt @@ -1,5 +1,5 @@ ############################################################################# -# Copyright (C) 2020-2024 MEmilio +# Copyright (C) 2020-2025 MEmilio # # Authors: Daniel Abele # From 17cd0b9caf021035f34211ad5ccc5166f975ebaa Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Thu, 2 Oct 2025 09:24:28 +0200 Subject: [PATCH 13/19] CHG: REplace FLOAT test by double test --- cpp/tests/test_smm_model.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp index 2351e33ef6..9e7be226f4 100644 --- a/cpp/tests/test_smm_model.cpp +++ b/cpp/tests/test_smm_model.cpp @@ -124,9 +124,9 @@ TEST(TestSMM, evaluateinterregionalAdoptionRate) model.populations[{mio::regions::Region(1), Infections::I, Age(1)}] = 16; model.populations[{mio::regions::Region(1), Infections::R, Age(1)}] = 0; - EXPECT_FLOAT_EQ(model.evaluate(adoption_rates[0], model.populations.get_compartments()), - 0.1 * 50. * (0.1 * 20. * 1. / 120. + 0.2 * 16 * 1. / 96.)); - EXPECT_FLOAT_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 4.); + EXPECT_DOUBLE_EQ(model.evaluate(adoption_rates[0], model.populations.get_compartments()), + 0.1 * 50. * (0.1 * 20. * 1. / 120. + 0.2 * 16 * 1. / 96.)); + EXPECT_DOUBLE_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 4.); } TEST(TestSMM, evaluateTransitionRate) From f37fbe1ae453a31db174709da3edb12c78916f7e Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Thu, 9 Oct 2025 09:49:39 +0200 Subject: [PATCH 14/19] Revert "[skip ci] Update copyright information" This reverts commit c58b74ae408e89e7d501b4d90bb11fc4163537b6. --- cpp/examples/d_abm.cpp | 2 +- cpp/examples/graph_abm.cpp | 2 +- cpp/examples/smm.cpp | 2 +- cpp/models/d_abm/model.cpp | 2 +- cpp/models/d_abm/model.h | 2 +- cpp/models/d_abm/parameters.h | 2 +- cpp/models/d_abm/quad_well.h | 2 +- cpp/models/d_abm/simulation.cpp | 2 +- cpp/models/d_abm/simulation.h | 2 +- cpp/models/d_abm/single_well.h | 2 +- cpp/models/graph_abm/graph_abm_mobility.cpp | 2 +- cpp/models/graph_abm/graph_abm_mobility.h | 2 +- cpp/models/graph_abm/graph_abmodel.h | 2 +- cpp/models/hybrid/temporal_hybrid_model.cpp | 2 +- cpp/models/hybrid/temporal_hybrid_model.h | 2 +- cpp/models/smm/model.cpp | 2 +- cpp/models/smm/model.h | 2 +- cpp/models/smm/parameters.h | 2 +- cpp/models/smm/simulation.h | 2 +- cpp/tests/test_d_abm_model.cpp | 2 +- cpp/tests/test_smm_model.cpp | 2 +- pycode/memilio-epidata/memilio/epidata/getJHData.py | 2 +- .../memilio/generation/template/template_cpp.txt | 2 +- .../memilio/generation/template/template_py.txt | 2 +- .../memilio/generation_test/test_data/test_oseir.cpp.txt | 2 +- .../memilio/generation_test/test_data/test_oseir.py.txt | 2 +- 26 files changed, 26 insertions(+), 26 deletions(-) diff --git a/cpp/examples/d_abm.cpp b/cpp/examples/d_abm.cpp index 4f91f65de0..09ffa04d5c 100644 --- a/cpp/examples/d_abm.cpp +++ b/cpp/examples/d_abm.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) * * Authors: Julia Bicker, René Schmieding * diff --git a/cpp/examples/graph_abm.cpp b/cpp/examples/graph_abm.cpp index dedc28b60e..3e7443fbc9 100644 --- a/cpp/examples/graph_abm.cpp +++ b/cpp/examples/graph_abm.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 MEmilio * * Authors: Julia Bicker * diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index cf13a0c7be..9eae9129a4 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) * * Authors: Julia Bicker, René Schmieding, Kilian Volmer * diff --git a/cpp/models/d_abm/model.cpp b/cpp/models/d_abm/model.cpp index 528fa1a4f4..f52bd1c298 100644 --- a/cpp/models/d_abm/model.cpp +++ b/cpp/models/d_abm/model.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: René Schmieding * diff --git a/cpp/models/d_abm/model.h b/cpp/models/d_abm/model.h index 5b858a4efa..b0983c5496 100644 --- a/cpp/models/d_abm/model.h +++ b/cpp/models/d_abm/model.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/d_abm/parameters.h b/cpp/models/d_abm/parameters.h index 5b16eb433f..bd02befd9b 100644 --- a/cpp/models/d_abm/parameters.h +++ b/cpp/models/d_abm/parameters.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/d_abm/quad_well.h b/cpp/models/d_abm/quad_well.h index ef1d0e021e..305617304e 100644 --- a/cpp/models/d_abm/quad_well.h +++ b/cpp/models/d_abm/quad_well.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/d_abm/simulation.cpp b/cpp/models/d_abm/simulation.cpp index 3915b84300..3b538490b8 100644 --- a/cpp/models/d_abm/simulation.cpp +++ b/cpp/models/d_abm/simulation.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 MEmilio * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/d_abm/simulation.h b/cpp/models/d_abm/simulation.h index e20f40f595..54ccad8b65 100644 --- a/cpp/models/d_abm/simulation.h +++ b/cpp/models/d_abm/simulation.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/d_abm/single_well.h b/cpp/models/d_abm/single_well.h index 7bd8cb82e1..62ef7671c6 100644 --- a/cpp/models/d_abm/single_well.h +++ b/cpp/models/d_abm/single_well.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: Julia Bicker, René Schmieding * diff --git a/cpp/models/graph_abm/graph_abm_mobility.cpp b/cpp/models/graph_abm/graph_abm_mobility.cpp index c8482f2e48..af3fae0c89 100644 --- a/cpp/models/graph_abm/graph_abm_mobility.cpp +++ b/cpp/models/graph_abm/graph_abm_mobility.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 MEmilio * * Authors: Julia Bicker * diff --git a/cpp/models/graph_abm/graph_abm_mobility.h b/cpp/models/graph_abm/graph_abm_mobility.h index 1af88d311a..2d7a749edf 100644 --- a/cpp/models/graph_abm/graph_abm_mobility.h +++ b/cpp/models/graph_abm/graph_abm_mobility.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 MEmilio * * Authors: Julia Bicker * diff --git a/cpp/models/graph_abm/graph_abmodel.h b/cpp/models/graph_abm/graph_abmodel.h index 040cd8c2da..0be1a7f222 100644 --- a/cpp/models/graph_abm/graph_abmodel.h +++ b/cpp/models/graph_abm/graph_abmodel.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 MEmilio * * Authors: Julia Bicker * diff --git a/cpp/models/hybrid/temporal_hybrid_model.cpp b/cpp/models/hybrid/temporal_hybrid_model.cpp index 337b8da3de..f356bd2321 100644 --- a/cpp/models/hybrid/temporal_hybrid_model.cpp +++ b/cpp/models/hybrid/temporal_hybrid_model.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: Julia Bicker, René Schmieding * diff --git a/cpp/models/hybrid/temporal_hybrid_model.h b/cpp/models/hybrid/temporal_hybrid_model.h index 3ba096c115..25ea35d8d0 100644 --- a/cpp/models/hybrid/temporal_hybrid_model.h +++ b/cpp/models/hybrid/temporal_hybrid_model.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: Julia Bicker, René Schmieding * diff --git a/cpp/models/smm/model.cpp b/cpp/models/smm/model.cpp index d613003ac2..af0d58eb73 100644 --- a/cpp/models/smm/model.cpp +++ b/cpp/models/smm/model.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: René Schmieding * diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index b5942bc251..7d6583bd08 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: René Schmieding, Julia Bicker, Kilian Volmer * diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h index fe27c35526..c7759d8953 100644 --- a/cpp/models/smm/parameters.h +++ b/cpp/models/smm/parameters.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 8fb1423571..8364236f28 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC) * * Authors: René Schmieding, Julia Bicker * diff --git a/cpp/tests/test_d_abm_model.cpp b/cpp/tests/test_d_abm_model.cpp index 79397d8e65..c31b4aa5e4 100644 --- a/cpp/tests/test_d_abm_model.cpp +++ b/cpp/tests/test_d_abm_model.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) * * Authors: Julia Bicker * diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp index 9e7be226f4..6e8901e37e 100644 --- a/cpp/tests/test_smm_model.cpp +++ b/cpp/tests/test_smm_model.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC) * * Authors: Julia Bicker, Kilian Volmer * diff --git a/pycode/memilio-epidata/memilio/epidata/getJHData.py b/pycode/memilio-epidata/memilio/epidata/getJHData.py index dd0a3bb873..19ec2b67f9 100644 --- a/pycode/memilio-epidata/memilio/epidata/getJHData.py +++ b/pycode/memilio-epidata/memilio/epidata/getJHData.py @@ -1,5 +1,5 @@ ############################################################################# -# Copyright (C) 2020-2025 MEmilio +# Copyright (C) 2020-2024 MEmilio # # Authors: Kathrin Rack # diff --git a/pycode/memilio-generation/memilio/generation/template/template_cpp.txt b/pycode/memilio-generation/memilio/generation/template/template_cpp.txt index e33668d1ee..28d9f80042 100644 --- a/pycode/memilio-generation/memilio/generation/template/template_cpp.txt +++ b/pycode/memilio-generation/memilio/generation/template/template_cpp.txt @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 MEmilio * * Authors: Martin Siggel, Daniel Abele, Martin J. Kuehn, Jan Kleinert * diff --git a/pycode/memilio-generation/memilio/generation/template/template_py.txt b/pycode/memilio-generation/memilio/generation/template/template_py.txt index d1cc9710c7..2ae2f78e5a 100644 --- a/pycode/memilio-generation/memilio/generation/template/template_py.txt +++ b/pycode/memilio-generation/memilio/generation/template/template_py.txt @@ -1,5 +1,5 @@ ############################################################################# -# Copyright (C) 2020-2025 MEmilio +# Copyright (C) 2020-2024 MEmilio # # Authors: Daniel Abele # diff --git a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt index edacafb623..b419020ba3 100644 --- a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt +++ b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.cpp.txt @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2025 MEmilio +* Copyright (C) 2020-2024 MEmilio * * Authors: Martin Siggel, Daniel Abele, Martin J. Kuehn, Jan Kleinert * diff --git a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt index 1fd1861e46..3f289a1251 100644 --- a/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt +++ b/pycode/memilio-generation/memilio/generation_test/test_data/test_oseir.py.txt @@ -1,5 +1,5 @@ ############################################################################# -# Copyright (C) 2020-2025 MEmilio +# Copyright (C) 2020-2024 MEmilio # # Authors: Daniel Abele # From d367215611613ac696cb6b5eea2a68dfb1c8fa01 Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Fri, 21 Nov 2025 16:56:02 +0100 Subject: [PATCH 15/19] generalize Region and State --- cpp/examples/smm.cpp | 85 +++++------ cpp/memilio/epidemiology/adoption_rate.h | 40 ++---- cpp/memilio/epidemiology/populations.h | 6 + cpp/memilio/utils/index.h | 159 +++++++++++++++------ cpp/models/hybrid/conversion_functions.cpp | 4 +- cpp/models/hybrid/conversion_functions.h | 7 +- cpp/models/smm/model.h | 80 ++++------- cpp/models/smm/parameters.h | 22 ++- cpp/models/smm/simulation.h | 64 +++------ cpp/tests/test_smm_model.cpp | 147 +++++++++---------- cpp/tests/test_temporal_hybrid_model.cpp | 17 +-- 11 files changed, 311 insertions(+), 320 deletions(-) diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index cf13a0c7be..89dd7185b3 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -18,12 +18,12 @@ * limitations under the License. */ +#include "memilio/config.h" #include "memilio/epidemiology/age_group.h" #include "smm/simulation.h" #include "smm/model.h" #include "smm/parameters.h" #include "memilio/data/analyze_result.h" -#include "memilio/epidemiology/adoption_rate.h" enum class InfectionState { @@ -36,78 +36,71 @@ enum class InfectionState Count }; -using Age = mio::AgeGroup; -using Species = mio::AgeGroup; +struct Species : public mio::Index { + Species(size_t val) + : Index(val) + { + } +}; int main() { + using Age = mio::AgeGroup; + using Status = mio::Index; + using mio::regions::Region; + using enum InfectionState; //Example how to run the stochastic metapopulation models with four regions const size_t num_regions = 4; const size_t num_age_groups = 1; - const size_t num_groups = 1; - using Model = mio::smm::Model; + const size_t num_species = 1; + using Model = mio::smm::Model; ScalarType numE = 12, numC = 4, numI = 12, numR = 0, numD = 0; - Model model; + Model model(Status{Count, Age(num_age_groups), Species(num_species)}, Region(num_regions)); //Population are distributed uniformly to the four regions for (size_t r = 0; r < num_regions; ++r) { - model.populations[{mio::regions::Region(r), InfectionState::S, Age(0), Species(0)}] = - (1000 - numE - numC - numI - numR - numD) / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::E, Age(0), Species(0)}] = numE / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::C, Age(0), Species(0)}] = numC / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::I, Age(0), Species(0)}] = numI / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::R, Age(0), Species(0)}] = numR / num_regions; - model.populations[{mio::regions::Region(r), InfectionState::D, Age(0), Species(0)}] = numD / num_regions; + model.populations[{Region(r), S, Age(0), Species(0)}] = (1000 - numE - numC - numI - numR - numD) / num_regions; + model.populations[{Region(r), E, Age(0), Species(0)}] = numE / num_regions; + model.populations[{Region(r), C, Age(0), Species(0)}] = numC / num_regions; + model.populations[{Region(r), I, Age(0), Species(0)}] = numI / num_regions; + model.populations[{Region(r), R, Age(0), Species(0)}] = numR / num_regions; + model.populations[{Region(r), D, Age(0), Species(0)}] = numD / num_regions; } + using AR = mio::smm::AdoptionRates; + using TR = mio::smm::TransitionRates; + //Set infection state adoption and spatial transition rates - std::vector> adoption_rates; - std::vector> transition_rates; + AR::Type adoption_rates; + TR::Type transition_rates; for (size_t r = 0; r < num_regions; ++r) { - adoption_rates.push_back({InfectionState::S, - InfectionState::E, - mio::regions::Region(r), + adoption_rates.push_back({{S, Age(0), Species(0)}, + {E, Age(0), Species(0)}, + Region(r), 0.1, - {{InfectionState::C, 1, mio::regions::Region(3), {Age(0), Species(0)}}, - {InfectionState::I, 0.5, mio::regions::Region(1), {Age(0), Species(0)}}}, - {Age(0), Species(0)}}); - adoption_rates.push_back( - {InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}, {Age(0), Species(0)}}); - adoption_rates.push_back( - {InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}, {Age(0), Species(0)}}); - adoption_rates.push_back( - {InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}, {Age(0), Species(0)}}); - adoption_rates.push_back( - {InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}, {Age(0), Species(0)}}); - adoption_rates.push_back( - {InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}, {Age(0), Species(0)}}); + {{{C, Age(0), Species(0)}, 1}, {{I, Age(0), Species(0)}, 0.5}}}); + adoption_rates.push_back({{C, Age(0), Species(0)}, {R, Age(0), Species(0)}, Region(r), 0.2 / 3., {}}); + adoption_rates.push_back({{E, Age(0), Species(0)}, {C, Age(0), Species(0)}, Region(r), 1.0 / 5., {}}); + adoption_rates.push_back({{C, Age(0), Species(0)}, {I, Age(0), Species(0)}, Region(r), 0.8 / 3., {}}); + adoption_rates.push_back({{I, Age(0), Species(0)}, {R, Age(0), Species(0)}, Region(r), 0.99 / 5., {}}); + adoption_rates.push_back({{I, Age(0), Species(0)}, {D, Age(0), Species(0)}, Region(r), 0.01 / 5., {}}); } //Agents in infection state D do not transition - for (size_t s = 0; s < static_cast(InfectionState::D); ++s) { + for (size_t s = 0; s < static_cast(D); ++s) { for (size_t i = 0; i < num_regions; ++i) { for (size_t j = 0; j < num_regions; ++j) if (i != j) { - transition_rates.push_back({InfectionState(s), - mio::regions::Region(i), - mio::regions::Region(j), - 0.01, - {Age(0), Species(0)}, - {Age(0), Species(0)}}); - transition_rates.push_back({InfectionState(s), - mio::regions::Region(j), - mio::regions::Region(i), - 0.01, - {Age(0), Species(0)}, - {Age(0), Species(0)}}); + transition_rates.push_back({{InfectionState(s), Age(0), Species(0)}, Region(i), Region(j), 0.01}); + transition_rates.push_back({{InfectionState(s), Age(0), Species(0)}, Region(j), Region(i), 0.01}); } } } - model.parameters.get>() = adoption_rates; - model.parameters.get>() = transition_rates; + model.parameters.get() = adoption_rates; + model.parameters.get() = transition_rates; ScalarType dt = 0.1; ScalarType tmax = 30.0; diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h index f1a89a8d2c..5a2a953ca1 100644 --- a/cpp/memilio/epidemiology/adoption_rate.h +++ b/cpp/memilio/epidemiology/adoption_rate.h @@ -20,16 +20,23 @@ #ifndef MIO_EPI_ADOPTIONRATE_H #define MIO_EPI_ADOPTIONRATE_H -#include "memilio/utils/index.h" -#include "memilio/config.h" #include "memilio/geography/regions.h" -#include -#include -#include +#include namespace mio { +/** + * @brief Struct defining an influence for a second-order adoption. + * The population having "status" is multiplied with "factor." + * @tparam Status An infection state enum. + */ +template +struct Influence { + Status status; + FP factor; +}; + /** * @brief Struct defining a possible status adoption in a Model based on Poisson Processes. * The AdoptionRate is considered to be of second-order if there are any "influences". @@ -38,30 +45,13 @@ namespace mio * @tparam Status An infection state enum. * @tparam Groups Additional grouping indices. */ -template +template struct AdoptionRate { - - /** - * @brief Struct defining an influence for a second-order adoption. - * The population having "status" is multiplied with "factor." - * @tparam status An infection state enum. - * @tparam factor Scaling factor for the influence. - * @tparam - * @tparam Groups Additional grouping indices. - */ - struct Influence { - Status status; - FP factor; - std::optional region = std::nullopt; - std::tuple group_indices{}; - }; - Status from; // i Status to; // j - mio::regions::Region region; // k + Region region; // k FP factor; // gammahat_{ij}^k - std::vector influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} ) - std::tuple group_indices{}; + std::vector> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} )}; }; } // namespace mio diff --git a/cpp/memilio/epidemiology/populations.h b/cpp/memilio/epidemiology/populations.h index f2da6092ad..24246b4471 100644 --- a/cpp/memilio/epidemiology/populations.h +++ b/cpp/memilio/epidemiology/populations.h @@ -287,6 +287,12 @@ class Populations : public CustomIndexArray, Categories...> } }; +template +class Populations> : public Populations +{ + using Populations::Populations; +}; + } // namespace mio #endif // MIO_EPI_POPULATIONS_H diff --git a/cpp/memilio/utils/index.h b/cpp/memilio/utils/index.h index cba56b1d57..70c5f75d38 100644 --- a/cpp/memilio/utils/index.h +++ b/cpp/memilio/utils/index.h @@ -20,15 +20,69 @@ #ifndef INDEX_H #define INDEX_H +#include "memilio/io/io.h" #include "memilio/utils/compiler_diagnostics.h" #include "memilio/utils/type_safe.h" +#include +#include +#include +#include + namespace mio { template class Index; +namespace details +{ + +/// @brief Function definition that accepts a MultiIndex, used for the definition of IsMultiIndex. +template +void is_multi_index_impl(Index); + +} // namespace details + +/// @brief A MultiIndex is an Index any number of categories. Does accept empty or single category indices. +template +concept IsMultiIndex = requires(Index i) { details::is_multi_index_impl(i); }; + +namespace details +{ + +/// @brief Obtain a tuple of singular indices from a Index or MultiIndex. +template +std::tuple...> get_tuple(const Index& i) +{ + if constexpr (sizeof...(Tags) == 1) { + return std::tuple(i); + } + else { + return i.indices; + } +} + +/// @brief Obtain a tuple of one index from an enum value. +template +std::tuple> get_tuple(Enum i) + requires std::is_enum::value +{ + return std::tuple(Index(i)); +} + +template + requires((std::is_enum_v || IsMultiIndex) && ...) +decltype(auto) merge_indices_impl(IndexArgs&&... args) +{ + return std::tuple_cat(details::get_tuple(args)...); +} + +template +Index tuple_to_index(std::tuple...>); + +} // namespace details + /** * @brief An Index with a single template parameter is a typesafe wrapper for size_t * that is associated with a Tag. It is used to index into a CustomIndexArray @@ -72,8 +126,8 @@ class MEMILIO_ENABLE_EBO Index : public TypeSafe::value, void>* = nullptr> - Index(Dummy val) + Index(CategoryTag val) + requires std::is_enum_v : TypeSafe>((size_t)val) { } @@ -132,6 +186,13 @@ class Index { } + template + requires(sizeof...(IndexArgs) > 1) + Index(IndexArgs&&... subindices) + : indices(details::merge_indices_impl(std::forward(subindices)...)) + { + } + private: Index(const std::tuple...>& _indices) : indices(_indices) @@ -197,70 +258,69 @@ struct index_of_type, Index> { static constexpr std::size_t value = 0; }; -// retrieves the Index at the Ith position for a Index with more than one Tag -template 1), void>* = nullptr> +// retrieves the Index at the Ith position for a Index +template constexpr typename std::tuple_element...>>::type& get(Index& i) noexcept { - return std::get(i.indices); -} - -// retrieves the Index at the Ith position for a Index with one Tag, equals identity function -template * = nullptr> -constexpr typename std::tuple_element...>>::type& -get(Index& i) noexcept -{ - static_assert(I == 0, "I must be equal to zero for an Index with just one template parameter"); - return i; -} - -// retrieves the Index at the Ith position for a Index with more than one Tag const version -template 1), void>* = nullptr> -constexpr typename std::tuple_element...>>::type const& -get(Index const& i) noexcept -{ - return std::get(i.indices); + if constexpr (sizeof...(CategoryTags) == 1) { + static_assert(I == 0, "I must be equal to zero for an Index with just one template parameter"); + return i; + } + else { + return std::get(i.indices); + } } -// retrieves the Index at the Ith position for a Index with one Tag, equals identity function const version -template * = nullptr> +// retrieves the Index at the Ith position for a Index, const version +template constexpr typename std::tuple_element...>>::type const& get(Index const& i) noexcept { - static_assert(I == 0, "I must be equal to zero for an Index with just one template parameter"); - return i; + if constexpr (sizeof...(CategoryTags) == 1) { + static_assert(I == 0, "I must be equal to zero for an Index with just one template parameter"); + return i; + } + else { + return std::get(i.indices); + } } // retrieves the Index for the tag Tag of a Index with more than one Tag -template 1), void>* = nullptr> +template constexpr Index& get(Index& i) noexcept { - return std::get>(i.indices); -} - -// retrieves the Index for the tag Tag of a Index with one Tag, equals identity function -template * = nullptr> -constexpr Index& get(Index& i) noexcept -{ - using IndexTag = std::tuple_element_t<0, std::tuple>; - static_assert(std::is_same::value, "Tags must match for an Index with just one template parameter"); - return i; + if constexpr (sizeof...(CategoryTags) == 1) { + using IndexTag = std::tuple_element_t<0, std::tuple>; + static_assert(std::is_same::value, + "Tags must match for an Index with just one template parameter"); + return i; + } + else { + return std::get>(i.indices); + } } -// retrieves the Index for the tag Tag for a Index with more than one Tag const version -template 1), void>* = nullptr> +// retrieves the Index for the tag Tag of a Index with more than one Tag +template constexpr Index const& get(Index const& i) noexcept { - return std::get>(i.indices); + if constexpr (sizeof...(CategoryTags) == 1) { + using IndexTag = std::tuple_element_t<0, std::tuple>; + static_assert(std::is_same::value, + "Tags must match for an Index with just one template parameter"); + return i; + } + else { + return std::get>(i.indices); + } } -// retrieves the Index for the tag Tag for a Index with one Tag, equals identity function const version -template * = nullptr> -constexpr Index const& get(Index const& i) noexcept +template +decltype(auto) merge_indices(IndexArgs&&... args) { - using IndexTag = std::tuple_element_t<0, std::tuple>; - static_assert(std::is_same::value, "Tags must match for an Index with just one template parameter"); - return i; + using MergedIndex = decltype(details::tuple_to_index(details::merge_indices_impl(std::declval()...))); + return MergedIndex(std::forward(args)...); } namespace details @@ -317,6 +377,13 @@ SubIndex reduce_index(const SuperIndex& index) return details::reduce_index_impl(index, mio::Tag{}); } +template + requires std::is_enum_v +Index reduce_index(const SuperIndex& index) +{ + return details::reduce_index_impl(index, mio::Tag>{}); +} + /** * @brief Create a SuperIndex by copying values from SubIndex, filling new categories with fill_value. * If a type T is contained multiple times in SubIndex, only the first occurance of T is used. diff --git a/cpp/models/hybrid/conversion_functions.cpp b/cpp/models/hybrid/conversion_functions.cpp index 8fa4b047b4..7fa8d4831a 100644 --- a/cpp/models/hybrid/conversion_functions.cpp +++ b/cpp/models/hybrid/conversion_functions.cpp @@ -36,7 +36,7 @@ namespace hybrid { template <> void convert_model(const dabm::Simulation>& current_model, - smm::Simulation& target_model) + smm::Simulation& target_model) { auto& current_result = current_model.get_result(); auto& target_result = target_model.get_result(); @@ -58,7 +58,7 @@ void convert_model(const dabm::Simulation -void convert_model(const smm::Simulation& current_model, +void convert_model(const smm::Simulation& current_model, dabm::Simulation>& target_model) { auto& current_result = current_model.get_result(); diff --git a/cpp/models/hybrid/conversion_functions.h b/cpp/models/hybrid/conversion_functions.h index 9d0e316543..2d15f6d638 100644 --- a/cpp/models/hybrid/conversion_functions.h +++ b/cpp/models/hybrid/conversion_functions.h @@ -37,11 +37,12 @@ namespace hybrid template <> void convert_model(const dabm::Simulation>& current_model, - smm::Simulation& target_model); + smm::Simulation>& target_model); template <> -void convert_model(const smm::Simulation& current_model, - dabm::Simulation>& target_model); +void convert_model( + const smm::Simulation>& current_model, + dabm::Simulation>& target_model); template <> void convert_model(const dabm::Simulation>& current_model, diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index b5942bc251..cd6be1577d 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -21,12 +21,13 @@ #ifndef MIO_SMM_MODEL_H #define MIO_SMM_MODEL_H -#include "memilio/config.h" -#include "memilio/epidemiology/age_group.h" +#include "memilio/utils/index.h" +#include "memilio/utils/index_range.h" #include "smm/parameters.h" #include "memilio/compartments/compartmental_model.h" #include "memilio/epidemiology/populations.h" #include "memilio/geography/regions.h" +#include namespace mio { @@ -36,26 +37,27 @@ namespace smm template using age_group = T; +template +using PopulationIndex = decltype(merge_indices(std::declval(), std::declval())); + /** * @brief Stochastic Metapopulation Model. * @tparam regions Number of regions. * @tparam Status An infection state enum. */ -template -class Model : public mio::CompartmentalModel< - FP, Status, mio::Populations...>, - ParametersBase...>> +template +class Model : public mio::CompartmentalModel>, + ParametersBase> { - using Base = - mio::CompartmentalModel...>, - ParametersBase...>>; - using Index = mio::Index...>; + using Base = mio::CompartmentalModel>, + ParametersBase>; public: - Model() - : Base(typename Base::Populations( - {static_cast(regions), Status::Count, static_cast(Groups)...}), + using Status = StatusT; + using Region = RegionT; + + Model(Status status_dimensions, Region region_dimensions) + : Base(typename Base::Populations(merge_indices(region_dimensions, status_dimensions)), typename Base::ParameterSet()) { } @@ -66,46 +68,26 @@ class Model : public mio::CompartmentalModel< * @param[in] x The current state of the model. * @return Current value of the adoption rate. */ - FP evaluate(const AdoptionRate...>& rate, - const Eigen::VectorXd& x) const + FP evaluate(const AdoptionRate& rate, const Eigen::VectorXd& x) const { - const auto& pop = this->populations; - const auto index_from = std::apply( - [&](auto&&... args) { - return Index{rate.region, rate.from, std::forward(args)...}; - }, - rate.group_indices); - const auto source = pop.get_flat_index(index_from); // Why is here rate.from used? KV + const auto& pop = this->populations; + const auto source = pop.get_flat_index({rate.region, rate.from}); // determine order and calculate rate if (rate.influences.size() == 0) { // first order adoption return rate.factor * x[source]; } else { // second order adoption + FP N = 0; + for (auto status : make_index_range(reduce_index(this->populations.size()))) { + N += x[pop.get_flat_index({rate.region, status})]; + } // accumulate influences FP influences = 0.0; for (size_t i = 0; i < rate.influences.size(); i++) { - FP N = 0; // Welches N brauchen wir hier?? - - for (size_t s = 0; s < static_cast(Status::Count); ++s) { - const auto index = std::apply( - [&](auto&&... args) { - return Index{rate.influences[i].region.value_or(rate.region), Status(s), - std::forward(args)...}; - }, - rate.influences[i].group_indices); - N += x[pop.get_flat_index(index)]; - } - const auto index = std::apply( - [&](auto&&... args) { - return Index{rate.influences[i].region.value_or(rate.region), rate.influences[i].status, - std::forward(args)...}; - }, - rate.influences[i].group_indices); - if (N > 0) { - influences += rate.influences[i].factor * x[pop.get_flat_index(index)] / N; - } + influences += + rate.influences[i].factor * x[pop.get_flat_index({rate.region, rate.influences[i].status})]; } - return rate.factor * x[source] * influences; + return (N > 0) ? (rate.factor * x[source] * influences / N) : 0; } } @@ -115,15 +97,9 @@ class Model : public mio::CompartmentalModel< * @param[in] x The current state of the model. * @return Current value of the transition rate. */ - FP evaluate(const TransitionRate...>& rate, - const Eigen::VectorXd& x) const + FP evaluate(const TransitionRate& rate, const Eigen::VectorXd& x) const { - auto index = std::apply( - [&](auto&&... args) { - return Index{rate.from, rate.status, std::forward(args)...}; - }, - rate.group_indices_from); - const auto source = this->populations.get_flat_index(index); + const auto source = this->populations.get_flat_index({rate.from, rate.status}); return rate.factor * x[source]; } diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h index fe27c35526..b44b63e827 100644 --- a/cpp/models/smm/parameters.h +++ b/cpp/models/smm/parameters.h @@ -25,7 +25,6 @@ #include "memilio/geography/regions.h" #include "memilio/utils/parameter_set.h" #include "memilio/epidemiology/adoption_rate.h" -#include namespace mio { @@ -36,9 +35,9 @@ namespace smm * @brief A vector of AdoptionRate%s, see mio::AdoptionRate * @tparam Status An infection state enum. */ -template +template struct AdoptionRates { - using Type = std::vector>; + using Type = std::vector>; const static std::string name() { return "AdoptionRates"; @@ -49,26 +48,25 @@ struct AdoptionRates { * @brief Struct defining a possible regional transition in a Model based on Poisson Processes. * @tparam Status An infection state enum. */ -template +template struct TransitionRate { Status status; // i - mio::regions::Region from; // k - mio::regions::Region to; // l + Region from; // k + Region to; // l FP factor; // lambda_i^{kl} - std::tuple group_indices_from{}; - std::tuple group_indices_to{}; }; -template + +template struct TransitionRates { - using Type = std::vector>; + using Type = std::vector>; const static std::string name() { return "TransitionRates"; } }; -template -using ParametersBase = mio::ParameterSet, TransitionRates>; +template +using ParametersBase = mio::ParameterSet, TransitionRates>; } // namespace smm diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 53096ef1e1..3f7f922b25 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -21,14 +21,9 @@ #ifndef MIO_SMM_SIMULATION_H #define MIO_SMM_SIMULATION_H -#include "memilio/config.h" -#include "memilio/epidemiology/age_group.h" -#include "memilio/utils/index.h" -#include "memilio/utils/logging.h" #include "smm/model.h" #include "smm/parameters.h" -#include "memilio/compartments/simulation.h" -#include +#include "memilio/utils/time_series.h" namespace mio { @@ -41,12 +36,12 @@ namespace smm * @tparam regions The number of regions. * @tparam Status An infection state enum. */ -template +template class Simulation { public: - using Model = smm::Model; - using Index = mio::Index...>; + using Model = smm::Model; + using Index = typename Model::Populations::Index; /** * @brief Set up the simulation for a Stochastic Metapopulation Model. @@ -104,38 +99,18 @@ class Simulation if (next_event < adoption_rates().size()) { // perform adoption event const auto& rate = adoption_rates()[next_event]; - auto index_from = std::apply( - [&](auto&&... args) { - return Index{rate.region, rate.from, std::forward(args)...}; - }, - rate.group_indices); - m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1; - m_model->populations[index_from] -= 1; - auto index_to = std::apply( - [&](auto&&... args) { - return Index{rate.region, rate.to, std::forward(args)...}; - }, - rate.group_indices); - m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1; - m_model->populations[index_to] += 1; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.from})] -= 1; + m_model->populations[{rate.region, rate.from}] -= 1.0; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.region, rate.to})] += 1; + m_model->populations[{rate.region, rate.to}] += 1.0; } else { // perform transition event const auto& rate = transition_rates()[next_event - adoption_rates().size()]; - auto index_from = std::apply( - [&](auto&&... args) { - return Index{rate.from, rate.status, std::forward(args)...}; - }, - rate.group_indices_from); - m_result.get_last_value()[m_model->populations.get_flat_index(index_from)] -= 1; - m_model->populations[index_from] -= 1; - auto index_to = std::apply( - [&](auto&&... args) { - return Index{rate.to, rate.status, std::forward(args)...}; - }, - rate.group_indices_to); - m_result.get_last_value()[m_model->populations.get_flat_index(index_to)] += 1; - m_model->populations[index_to] += 1; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.from, rate.status})] -= 1; + m_model->populations[{rate.from, rate.status}] -= 1.0; + m_result.get_last_value()[m_model->populations.get_flat_index({rate.to, rate.status})] += 1; + m_model->populations[{rate.to, rate.status}] += 1.0; } // update internal times for (size_t i = 0; i < m_internal_time.size(); i++) { @@ -157,7 +132,6 @@ class Simulation } } return m_result.get_last_value(); - // } } /** @@ -189,25 +163,21 @@ class Simulation /** * @brief Returns the model's transition rates. */ - inline constexpr const typename smm::TransitionRates...>::Type& - transition_rates() + inline constexpr const typename smm::TransitionRates::Type& transition_rates() { - return m_model->parameters - .template get...>>(); + return m_model->parameters.template get>(); } /** * @brief Returns the model's adoption rates. */ - inline constexpr const typename smm::AdoptionRates...>::Type& - adoption_rates() + inline constexpr const typename smm::AdoptionRates::Type& adoption_rates() { - return m_model->parameters.template get...>>(); + return m_model->parameters.template get>(); } /** - * @brief - * + * @brief Calculate current values for m_current_rates and m_waiting_times. */ inline void update_current_rates_and_waiting_times() { diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp index 9e7be226f4..56386c73d0 100644 --- a/cpp/tests/test_smm_model.cpp +++ b/cpp/tests/test_smm_model.cpp @@ -18,20 +18,17 @@ * limitations under the License. */ -#include -#include -#include -#include -#include +#include "abm_helpers.h" +#include "memilio/epidemiology/adoption_rate.h" #include "memilio/epidemiology/age_group.h" -#include "memilio/utils/compiler_diagnostics.h" -#include "memilio/utils/logging.h" +#include "memilio/geography/regions.h" #include "smm/model.h" #include "smm/parameters.h" #include "smm/simulation.h" -#include "abm_helpers.h" -#include "memilio/epidemiology/adoption_rate.h" -#include + +#include + +#include #include #include @@ -63,9 +60,9 @@ TEST(TestSMM, evaluateAdoptionRate) // with N(from, region) the population in Region "region" having infection state "from" //Second-order adoption rates are given by: // rate.factor * N(rate.from, rate.region)/total_pop * sum (over all rate.influences) influence.factor * N(influence.status, rate.region) - using Model = mio::smm::Model; + using Model = mio::smm::Model; - Model model; + Model model(InfectionState::Count, mio::regions::Region(1)); //Set adoption rates std::vector> adoption_rates; @@ -92,37 +89,38 @@ TEST(TestSMM, evaluateAdoptionRate) TEST(TestSMM, evaluateinterregionalAdoptionRate) { - using Model = mio::smm::Model; - using Age = mio::AgeGroup; + using Age = mio::AgeGroup; + using mio::regions::Region; + using Status = mio::Index; + using Model = mio::smm::Model; - std::vector> adoption_rates; + std::vector> adoption_rates; //Second-order adoption - adoption_rates.push_back({Infections::S, - Infections::I, + adoption_rates.push_back({{Infections::S, Age(0), Region(0)}, + {Infections::I, Age(0), Region(0)}, mio::regions::Region(0), 0.1, - {{Infections::I, 0.1, mio::regions::Region(0), {Age(1)}}, - {Infections::I, 0.2, mio::regions::Region(1), {Age(1)}}}, - Age(0)}); + {{{Infections::I, Age(1), Region(0)}, 0.1}, {{Infections::I, Age(1), Region(1)}, 0.2}}}); //First-order adoption - adoption_rates.push_back({Infections::I, Infections::R, mio::regions::Region(0), 0.2, {}, {Age(1)}}); - Model model; + adoption_rates.push_back( + {{Infections::I, Age(1), Region(0)}, {Infections::R, Age(1), Region(0)}, Region(0), 0.2, {}}); + Model model({Infections::Count, Age(2), Region(2)}, Region(1)); - model.populations[{mio::regions::Region(0), Infections::S, Age(0)}] = 50; - model.populations[{mio::regions::Region(0), Infections::I, Age(0)}] = 10; - model.populations[{mio::regions::Region(0), Infections::R, Age(0)}] = 0; + model.populations[{Region(0), Infections::S, Age(0), Region(0)}] = 50; + model.populations[{Region(0), Infections::I, Age(0), Region(0)}] = 10; + model.populations[{Region(0), Infections::R, Age(0), Region(0)}] = 0; - model.populations[{mio::regions::Region(0), Infections::S, Age(1)}] = 100; - model.populations[{mio::regions::Region(0), Infections::I, Age(1)}] = 20; - model.populations[{mio::regions::Region(0), Infections::R, Age(1)}] = 0; + model.populations[{Region(0), Infections::S, Age(1), Region(0)}] = 100; + model.populations[{Region(0), Infections::I, Age(1), Region(0)}] = 20; + model.populations[{Region(0), Infections::R, Age(1), Region(0)}] = 0; - model.populations[{mio::regions::Region(1), Infections::S, Age(0)}] = 40; - model.populations[{mio::regions::Region(1), Infections::I, Age(0)}] = 80; - model.populations[{mio::regions::Region(1), Infections::R, Age(0)}] = 0; + model.populations[{Region(0), Infections::S, Age(0), Region(1)}] = 40; + model.populations[{Region(0), Infections::I, Age(0), Region(1)}] = 80; + model.populations[{Region(0), Infections::R, Age(0), Region(1)}] = 0; - model.populations[{mio::regions::Region(1), Infections::S, Age(1)}] = 80; - model.populations[{mio::regions::Region(1), Infections::I, Age(1)}] = 16; - model.populations[{mio::regions::Region(1), Infections::R, Age(1)}] = 0; + model.populations[{Region(0), Infections::S, Age(1), Region(1)}] = 80; + model.populations[{Region(0), Infections::I, Age(1), Region(1)}] = 16; + model.populations[{Region(0), Infections::R, Age(1), Region(1)}] = 0; EXPECT_DOUBLE_EQ(model.evaluate(adoption_rates[0], model.populations.get_compartments()), 0.1 * 50. * (0.1 * 20. * 1. / 120. + 0.2 * 16 * 1. / 96.)); @@ -133,37 +131,27 @@ TEST(TestSMM, evaluateTransitionRate) { //Same test as 'evaluateAdoptionRate' only for spatial transition rates. //Transition rates are given by: rate.factor * N(rate.status, rate.from) - using Model = mio::smm::Model; + using Model = mio::smm::Model; - Model model; + Model model(InfectionState::Count, mio::regions::Region(2)); //Initialize model populations - model.populations[{mio::regions::Region(0), InfectionState::S, mio::AgeGroup(0)}] = 50; - model.populations[{mio::regions::Region(0), InfectionState::E, mio::AgeGroup(0)}] = 10; - model.populations[{mio::regions::Region(0), InfectionState::C, mio::AgeGroup(0)}] = 5; - model.populations[{mio::regions::Region(0), InfectionState::I, mio::AgeGroup(0)}] = 0; - model.populations[{mio::regions::Region(0), InfectionState::R, mio::AgeGroup(0)}] = 0; - model.populations[{mio::regions::Region(0), InfectionState::D, mio::AgeGroup(0)}] = 0; - - model.populations[{mio::regions::Region(1), InfectionState::S, mio::AgeGroup(0)}] = 55; - model.populations[{mio::regions::Region(1), InfectionState::E, mio::AgeGroup(0)}] = 10; - model.populations[{mio::regions::Region(1), InfectionState::C, mio::AgeGroup(0)}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::I, mio::AgeGroup(0)}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::R, mio::AgeGroup(0)}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::D, mio::AgeGroup(0)}] = 0; + model.populations[{mio::regions::Region(0), InfectionState::S}] = 50; + model.populations[{mio::regions::Region(0), InfectionState::E}] = 10; + model.populations[{mio::regions::Region(0), InfectionState::C}] = 5; + model.populations[{mio::regions::Region(0), InfectionState::I}] = 0; + model.populations[{mio::regions::Region(0), InfectionState::R}] = 0; + model.populations[{mio::regions::Region(0), InfectionState::D}] = 0; + + model.populations[{mio::regions::Region(1), InfectionState::S}] = 55; + model.populations[{mio::regions::Region(1), InfectionState::E}] = 10; + model.populations[{mio::regions::Region(1), InfectionState::C}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::I}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::R}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::D}] = 0; //Set transition rates - std::vector> transition_rates; - transition_rates.push_back({InfectionState::S, - mio::regions::Region(0), - mio::regions::Region(1), - 0.01, - {mio::AgeGroup(0)}, - {mio::AgeGroup(0)}}); - transition_rates.push_back({InfectionState::E, - mio::regions::Region(1), - mio::regions::Region(0), - 0.1, - {mio::AgeGroup(0)}, - {mio::AgeGroup(0)}}); + std::vector> transition_rates; + transition_rates.push_back({InfectionState::S, mio::regions::Region(0), mio::regions::Region(1), 0.01}); + transition_rates.push_back({InfectionState::E, mio::regions::Region(1), mio::regions::Region(0), 0.1}); EXPECT_EQ(model.evaluate(transition_rates[0], model.populations.get_compartments()), 0.5); EXPECT_EQ(model.evaluate(transition_rates[1], model.populations.get_compartments()), 1.); @@ -173,9 +161,9 @@ TEST(TestSMMSimulation, advance) { //Test whether Gillespie algorithm calculates events in the correct order using testing::Return; - using Model = mio::smm::Model; + using Model = mio::smm::Model; - Model model; + Model model(InfectionState::Count, mio::regions::Region(2)); //Initialize model populations model.populations[{mio::regions::Region(0), InfectionState::S}] = 1; model.populations[{mio::regions::Region(0), InfectionState::E}] = 0; @@ -210,8 +198,8 @@ TEST(TestSMMSimulation, advance) //Spatial transition transition_rates.push_back({InfectionState::R, mio::regions::Region(1), mio::regions::Region(0), 0.01}); - model.parameters.get>() = adoption_rates; - model.parameters.get>() = transition_rates; + model.parameters.get>() = adoption_rates; + model.parameters.get>() = transition_rates; //Mock exponential distribution to control the normalized waiting times that are drawn ScopedMockDistribution>>> @@ -254,9 +242,9 @@ TEST(TestSMMSimulation, stopsAtTmax) { //Test whether simulation stops at tmax and whether system state at tmax is saved if there are no adoptions/transitions using testing::Return; - using Model = mio::smm::Model; + using Model = mio::smm::Model; - Model model; + Model model(InfectionState::Count, mio::regions::Region(2)); //Set adoption and spatial transition rates std::vector> adoption_rates; @@ -270,8 +258,8 @@ TEST(TestSMMSimulation, stopsAtTmax) transition_rates.push_back({InfectionState::R, mio::regions::Region(1), mio::regions::Region(0), 0.01}); - model.parameters.get>() = adoption_rates; - model.parameters.get>() = transition_rates; + model.parameters.get>() = adoption_rates; + model.parameters.get>() = transition_rates; //As populations are not set they have value 0 i.e. no events will happen //Simulate for 30 days @@ -287,7 +275,7 @@ TEST(TestSMMSimulation, convergence) //Test whether the mean number of transitions in one unit-time step corresponds to the expected number of transitions // given by rate * pop up to some tolerance using testing::Return; - using Model = mio::smm::Model; + using Model = mio::smm::Model; double pop = 1000; size_t num_runs = 100; @@ -302,7 +290,7 @@ TEST(TestSMMSimulation, convergence) //Do 100 unit-time step simulations with 1000 agents for (size_t n = 0; n < num_runs; ++n) { - Model model; + Model model(InfectionState::Count, mio::regions::Region(2)); model.populations[{mio::regions::Region(0), InfectionState::S}] = pop; model.populations[{mio::regions::Region(0), InfectionState::E}] = 0; @@ -311,13 +299,14 @@ TEST(TestSMMSimulation, convergence) model.populations[{mio::regions::Region(0), InfectionState::R}] = 0; model.populations[{mio::regions::Region(0), InfectionState::D}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::S}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::E}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::C}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::I}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::R}] = 0; - model.populations[{mio::regions::Region(1), InfectionState::D}] = 0; - model.parameters.get>() = transition_rates; + model.populations[{mio::regions::Region(1), InfectionState::S}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::E}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::C}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::I}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::R}] = 0; + model.populations[{mio::regions::Region(1), InfectionState::D}] = 0; + model.parameters.get>() = + transition_rates; auto sim = mio::smm::Simulation(model, 0.0, 1.0); sim.advance(1.); diff --git a/cpp/tests/test_temporal_hybrid_model.cpp b/cpp/tests/test_temporal_hybrid_model.cpp index 0ae31660ec..5f0ff64891 100644 --- a/cpp/tests/test_temporal_hybrid_model.cpp +++ b/cpp/tests/test_temporal_hybrid_model.cpp @@ -151,7 +151,7 @@ TEST(TestTemporalHybrid, test_advance) TEST(TestTemporalHybrid, test_conversion_dabm_smm) { using Model1 = mio::dabm::Model>; - using Model2 = mio::smm::Model; + using Model2 = mio::smm::Model; //Initialize agents for dabm SingleWell::Agent a1{Eigen::Vector2d{-0.5, 0}, @@ -162,13 +162,14 @@ TEST(TestTemporalHybrid, test_conversion_dabm_smm) mio::osecir::InfectionState::InfectedSymptoms}; Model1 model1({a1, a2, a3}, {}); - Model2 model2; - model2.parameters.get>().push_back( - {mio::osecir::InfectionState::Susceptible, - mio::osecir::InfectionState::Exposed, - mio::regions::Region(0), - 0.1, - {{mio::osecir::InfectionState::InfectedNoSymptoms, 1}, {mio::osecir::InfectionState::InfectedSymptoms, 0.5}}}); + Model2 model2(mio::osecir::InfectionState::Count, mio::regions::Region(1)); + model2.parameters.get>() + .push_back({mio::osecir::InfectionState::Susceptible, + mio::osecir::InfectionState::Exposed, + mio::regions::Region(0), + 0.1, + {{mio::osecir::InfectionState::InfectedNoSymptoms, 1}, + {mio::osecir::InfectionState::InfectedSymptoms, 0.5}}}); //Parameters for simulation double t0 = 0; From 72a33049421be6e083d2953faab09bb5d265c29c Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Tue, 25 Nov 2025 14:59:59 +0100 Subject: [PATCH 16/19] fix asserts for debug builds --- cpp/memilio/utils/index.h | 39 ++++++++++++++++++++++++++++++++----- cpp/models/smm/model.h | 1 + cpp/models/smm/simulation.h | 15 +++++++------- 3 files changed, 43 insertions(+), 12 deletions(-) diff --git a/cpp/memilio/utils/index.h b/cpp/memilio/utils/index.h index 70c5f75d38..8388aaba8b 100644 --- a/cpp/memilio/utils/index.h +++ b/cpp/memilio/utils/index.h @@ -22,6 +22,7 @@ #include "memilio/io/io.h" #include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/metaprogramming.h" #include "memilio/utils/type_safe.h" #include @@ -116,9 +117,10 @@ class MEMILIO_ENABLE_EBO Index : public TypeSafe>::TypeSafe; - static size_t constexpr size = 1; + static constexpr size_t size = 1; + static constexpr bool has_duplicates = false; - static Index constexpr Zero() + static constexpr Index Zero() { return Index((size_t)0); } @@ -173,7 +175,8 @@ template class Index { public: - static size_t constexpr size = sizeof...(CategoryTag); + static constexpr size_t size = sizeof...(CategoryTag); + static constexpr bool has_duplicates = has_duplicates_v; static Index constexpr Zero() { @@ -211,6 +214,25 @@ class Index return !(this == other); } + bool operator<(Index const& other) const + { + // use apply to unfold both tuples, then use a fold expression to evaluate a pairwise less + return std::apply( + [&other](auto&&... indices_) { + return std::apply( + [&](auto&&... other_indices_) { + return ((indices_ < other_indices_) && ...); + }, + other.indices); + }, + indices); + } + + bool operator<=(Index const& other) const + { + return (*this == other) || (*this < other); + } + /** * serialize this. * @see mio::serialize @@ -372,9 +394,16 @@ inline Index extend_index_impl(const Index& i, const * @return A (sub)index with the given categories and values from index. */ template -SubIndex reduce_index(const SuperIndex& index) +decltype(auto) reduce_index(const SuperIndex& index) { - return details::reduce_index_impl(index, mio::Tag{}); + if constexpr (SubIndex::size == 1 && std::is_base_of_v, SubIndex>) { + // this case handles reducing from e.g. Index directly to AgeGroup + // the default case would instead reduce to Index, which may cause conversion errors + return details::reduce_index_impl(index, mio::Tag>{}); + } + else { + return details::reduce_index_impl(index, mio::Tag{}); + } } template diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index cd6be1577d..0fc0eca2f6 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -51,6 +51,7 @@ class Model : public mio::CompartmentalModel>, ParametersBase>; + static_assert(!Base::Populations::Index::has_duplicates, "Do not use the same Index tag multiple times!"); public: using Status = StatusT; diff --git a/cpp/models/smm/simulation.h b/cpp/models/smm/simulation.h index 3f7f922b25..0ed56ee13c 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -41,7 +41,6 @@ class Simulation { public: using Model = smm::Model; - using Index = typename Model::Populations::Index; /** * @brief Set up the simulation for a Stochastic Metapopulation Model. @@ -60,12 +59,14 @@ class Simulation { assert(dt > 0); assert(m_waiting_times.size() > 0); - assert(std::all_of(adoption_rates().begin(), adoption_rates().end(), [](auto&& r) { - return static_cast(r.region) < regions; - })); - assert(std::all_of(transition_rates().begin(), transition_rates().end(), [](auto&& r) { - return static_cast(r.from) < regions && static_cast(r.to) < regions; - })); + assert(std::all_of(adoption_rates().begin(), adoption_rates().end(), + [regions = reduce_index(model.populations.size())](auto&& r) { + return r.region < regions; + })); + assert(std::all_of(transition_rates().begin(), transition_rates().end(), + [regions = reduce_index(model.populations.size())](auto&& r) { + return r.from < regions && r.to < regions; + })); // initialize (internal) next event times by random values for (size_t i = 0; i < m_tp_next_event.size(); i++) { m_tp_next_event[i] += mio::ExponentialDistribution::get_instance()(m_model->get_rng(), 1.0); From d47625cd468c888825a40d06d5557d8dda0f93bb Mon Sep 17 00:00:00 2001 From: reneSchm <49305466+reneSchm@users.noreply.github.com> Date: Tue, 25 Nov 2025 15:00:27 +0100 Subject: [PATCH 17/19] change test to use region as part of status --- cpp/tests/test_smm_model.cpp | 37 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp index 56386c73d0..0e906dc13f 100644 --- a/cpp/tests/test_smm_model.cpp +++ b/cpp/tests/test_smm_model.cpp @@ -92,38 +92,37 @@ TEST(TestSMM, evaluateinterregionalAdoptionRate) using Age = mio::AgeGroup; using mio::regions::Region; using Status = mio::Index; - using Model = mio::smm::Model; + using Model = mio::smm::Model>; - std::vector> adoption_rates; + std::vector>> adoption_rates; //Second-order adoption adoption_rates.push_back({{Infections::S, Age(0), Region(0)}, {Infections::I, Age(0), Region(0)}, - mio::regions::Region(0), + {}, 0.1, {{{Infections::I, Age(1), Region(0)}, 0.1}, {{Infections::I, Age(1), Region(1)}, 0.2}}}); //First-order adoption - adoption_rates.push_back( - {{Infections::I, Age(1), Region(0)}, {Infections::R, Age(1), Region(0)}, Region(0), 0.2, {}}); - Model model({Infections::Count, Age(2), Region(2)}, Region(1)); + adoption_rates.push_back({{Infections::I, Age(1), Region(0)}, {Infections::R, Age(1), Region(0)}, {}, 0.2, {}}); + Model model({Infections::Count, Age(2), Region(2)}, {}); - model.populations[{Region(0), Infections::S, Age(0), Region(0)}] = 50; - model.populations[{Region(0), Infections::I, Age(0), Region(0)}] = 10; - model.populations[{Region(0), Infections::R, Age(0), Region(0)}] = 0; + model.populations[{Infections::S, Age(0), Region(0)}] = 50; + model.populations[{Infections::I, Age(0), Region(0)}] = 10; + model.populations[{Infections::R, Age(0), Region(0)}] = 0; - model.populations[{Region(0), Infections::S, Age(1), Region(0)}] = 100; - model.populations[{Region(0), Infections::I, Age(1), Region(0)}] = 20; - model.populations[{Region(0), Infections::R, Age(1), Region(0)}] = 0; + model.populations[{Infections::S, Age(1), Region(0)}] = 100; + model.populations[{Infections::I, Age(1), Region(0)}] = 20; + model.populations[{Infections::R, Age(1), Region(0)}] = 0; - model.populations[{Region(0), Infections::S, Age(0), Region(1)}] = 40; - model.populations[{Region(0), Infections::I, Age(0), Region(1)}] = 80; - model.populations[{Region(0), Infections::R, Age(0), Region(1)}] = 0; + model.populations[{Infections::S, Age(0), Region(1)}] = 40; + model.populations[{Infections::I, Age(0), Region(1)}] = 80; + model.populations[{Infections::R, Age(0), Region(1)}] = 0; - model.populations[{Region(0), Infections::S, Age(1), Region(1)}] = 80; - model.populations[{Region(0), Infections::I, Age(1), Region(1)}] = 16; - model.populations[{Region(0), Infections::R, Age(1), Region(1)}] = 0; + model.populations[{Infections::S, Age(1), Region(1)}] = 80; + model.populations[{Infections::I, Age(1), Region(1)}] = 16; + model.populations[{Infections::R, Age(1), Region(1)}] = 0; EXPECT_DOUBLE_EQ(model.evaluate(adoption_rates[0], model.populations.get_compartments()), - 0.1 * 50. * (0.1 * 20. * 1. / 120. + 0.2 * 16 * 1. / 96.)); + 0.1 * 50. * (0.1 * 20. * 1. + 0.2 * 16 * 1.) / (60 + 120 + 120 + 96)); EXPECT_DOUBLE_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 4.); } From d2deb8fc145b7799c67e6235d7589c32fd337f43 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Mon, 15 Dec 2025 10:56:03 +0100 Subject: [PATCH 18/19] CHG: improve docstrings --- cpp/memilio/epidemiology/adoption_rate.h | 2 +- cpp/memilio/epidemiology/populations.h | 7 +++++++ cpp/models/smm/model.h | 1 + cpp/models/smm/parameters.h | 13 +++++++++++++ 4 files changed, 22 insertions(+), 1 deletion(-) diff --git a/cpp/memilio/epidemiology/adoption_rate.h b/cpp/memilio/epidemiology/adoption_rate.h index 5a2a953ca1..9ad4340d9f 100644 --- a/cpp/memilio/epidemiology/adoption_rate.h +++ b/cpp/memilio/epidemiology/adoption_rate.h @@ -43,7 +43,7 @@ struct Influence { * In the d_abm and smm simulations, "from" is implicitly an influence, scaled by "factor". This is multiplied by * the sum over all "influences", which scale their "status" with the respective "factor". * @tparam Status An infection state enum. - * @tparam Groups Additional grouping indices. + * @tparam Region An (multi)-index. */ template struct AdoptionRate { diff --git a/cpp/memilio/epidemiology/populations.h b/cpp/memilio/epidemiology/populations.h index 24246b4471..6bdefaa841 100644 --- a/cpp/memilio/epidemiology/populations.h +++ b/cpp/memilio/epidemiology/populations.h @@ -287,6 +287,13 @@ class Populations : public CustomIndexArray, Categories...> } }; +/** + * @brief Population template specialization for Index types. + * + * @tparam FP Floating point type + * @tparam Categories Index categories + */ + template class Populations> : public Populations { diff --git a/cpp/models/smm/model.h b/cpp/models/smm/model.h index 0fc0eca2f6..ac4cd4272e 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.h @@ -44,6 +44,7 @@ using PopulationIndex = decltype(merge_indices(std::declval(), std::decl * @brief Stochastic Metapopulation Model. * @tparam regions Number of regions. * @tparam Status An infection state enum. + * @tparam Region An (multi)-index, default is @ref mio::regions::Region. */ template class Model : public mio::CompartmentalModel>, diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h index b44b63e827..b27e06cb39 100644 --- a/cpp/models/smm/parameters.h +++ b/cpp/models/smm/parameters.h @@ -33,7 +33,10 @@ namespace smm /** * @brief A vector of AdoptionRate%s, see mio::AdoptionRate + * + * @tparam FP Floating point type * @tparam Status An infection state enum. + * @tparam Region An (multi)-index. */ template struct AdoptionRates { @@ -46,7 +49,10 @@ struct AdoptionRates { /** * @brief Struct defining a possible regional transition in a Model based on Poisson Processes. + * + * @tparam FP Floating point type * @tparam Status An infection state enum. + * @tparam Region An (multi)-index. */ template struct TransitionRate { @@ -56,6 +62,13 @@ struct TransitionRate { FP factor; // lambda_i^{kl} }; +/** + * @brief A vector of TransitionRate%s, see mio::TransitionRate + * + * @tparam FP Floating point type + * @tparam Status An infection state enum. + * @tparam Region An (multi)-index. + */ template struct TransitionRates { using Type = std::vector>; From 168116b6927e16667886ffc97d081bd4517e4310 Mon Sep 17 00:00:00 2001 From: Kilian Volmer <13285635+kilianvolmer@users.noreply.github.com> Date: Mon, 15 Dec 2025 15:45:18 +0100 Subject: [PATCH 19/19] CHG: Update smm.rst --- docs/source/cpp/smm.rst | 128 ++++++++++++++++++++++++++++++---------- 1 file changed, 98 insertions(+), 30 deletions(-) diff --git a/docs/source/cpp/smm.rst b/docs/source/cpp/smm.rst index e2fb413f75..b1a1901f3c 100644 --- a/docs/source/cpp/smm.rst +++ b/docs/source/cpp/smm.rst @@ -1,11 +1,20 @@ Stochastic metapopulation model =============================== -The stochastic metapopulation model uses a Markov process to simulate disease dynamics. Similar to the `Diffusive Agent-based Model` the Markov process is given by location and infection state changes. However, in contrast to the diffusive ABM, location changes are not given by agent-dependent diffusion processes, but by stochastic jumps between regions with the requirement that the domain is split into disjoint regions. Hence, there is no further spatial resolution within one region and locations or positions are only given by the region index. The evolution of the system state is determined by the following master equation +The stochastic metapopulation model uses a Markov process to simulate disease dynamics. Similar to the :doc:`diffusive_abm` the +Markov process is given by location and infection state changes. However, in contrast to the diffusive ABM, location changes are +not given by agent-dependent diffusion processes, but by stochastic jumps between regions with the requirement that the domain +is split into disjoint regions. Hence, there is no further spatial resolution within one region and locations or positions are +only given by the region index. The evolution of the system state is determined by the following master equation :math:`\partial_t p(X,Z;t) = G p(X,Z;t) + L p(X,Z;t)`. -The operator :math:`G` defines the infection state adoptions and only acts on :math:`Z`, the vector containing all subpopulations stratified by infection state. :math:`L` defines location changes, only acting on :math:`X`, the vector containing all subpopulations stratified by region. Infection state adoptions are modeled as stochastic jumps with independent Poisson processes given by adoption rate functions. Similar to the infection state dynamics, spatial transitions between regions are also modeled as stochastic jumps with independent Poisson processes given by transition rate functions. Gillespie's direct method (stochastic simulation algorithm) is used for simulation. +The operator :math:`G` defines the infection state adoptions and only acts on :math:`Z`, the vector containing all subpopulations +stratified by infection state. :math:`L` defines location changes, only acting on :math:`X`, the vector containing all subpopulations +stratified by region. Infection state adoptions are modeled as stochastic jumps with independent Poisson processes given by adoption +rate functions. Similar to the infection state dynamics, spatial transitions between regions are also modeled as stochastic jumps +with independent Poisson processes given by transition rate functions. Gillespie's direct method (stochastic simulation algorithm) +is used for the simulation. In the following, we present more details of the stochastic metapopulation model, including code examples. An overview of nonstandard but often used data types can be found under :doc:`data_types`. @@ -13,7 +22,8 @@ An overview of nonstandard but often used data types can be found under :doc:`da Infection states ---------------- -The model does not have fixed infection states, but gets an enum class of infection states as template argument. Thus it can be used with any set of infection states. +The model does not have fixed infection states, but gets an enum class of infection states as template argument. Thus it +can be used with any set of infection states. Using the infection states Susceptible (S), Exposed (E), Carrier (C), Infected (I), Recovered (R) and Dead (D), this can be done as follows: .. code-block:: cpp @@ -32,41 +42,44 @@ Using the infection states Susceptible (S), Exposed (E), Carrier (C), Infected ( const size_t num_regions = 2; - using Model = mio::smm::Model; + using Status = mio::Index; + using Model = mio::smm::Model; Model model; Infection state transitions --------------------------- -The infection state transitions are explicitly given by the adoption rates and are therefore subject to user input. Adoption rates always depend on their source infection state. If an adoption event requires interaction of agents (e.g. disease transmission), the corresponding rate depends not only on the source infection state, but also on other infection states, the **Influence**\s. An adoption rate that only depends on the source infection state, e.g. recovery or worsening of disease symptoms, is called `first-order` adoption rate and an adoption rate that has influences is called `second-order` adoption rate. Adoption rates are region-dependent; therefore it is possible to have different rates in two regions for the same infection state transition which can be useful when having e.g. region-dependent interventions or contact behavior. +The infection state transitions are explicitly given by the adoption rates and are therefore subject to user input. +Adoption rates always depend on their source infection state. If an adoption event requires interaction of agents (e.g. +disease transmission), the corresponding rate depends not only on the source infection state, but also on other infection +states, the **Influence**\s. An adoption rate that only depends on the source infection state, e.g. recovery or worsening +of disease symptoms, is called `first-order` adoption rate and an adoption rate that has influences is called `second-order` +adoption rate. Adoption rates are region-dependent; therefore it is possible to have different rates in two regions for +the same infection state transition which can be useful when having e.g. region-dependent interventions or contact behavior. -Using the infection states from above and two regions, there are five first-order adoption rates per region and one second-order adoption rate per region. In the example below, the second-order adoption rate (transition from S to E) differs between the regions: +Using the infection states from above and two regions, there are five first-order adoption rates per region and one second-order +adoption rate per region. In the example below, the second-order adoption rate (transition from S to E) differs between the regions: .. code-block:: cpp - std::vector> adoption_rates; + std::vector> adoption_rates; //Set first-order adoption rates for both regions for (size_t r = 0; r < num_regions; ++r) { - adoption_rates.push_back({InfectionState::E, InfectionState::C, mio::regions::Region(r), 1.0 / 5., {}}); - adoption_rates.push_back({InfectionState::C, InfectionState::R, mio::regions::Region(r), 0.2 / 3., {}}); - adoption_rates.push_back({InfectionState::C, InfectionState::I, mio::regions::Region(r), 0.8 / 3., {}}); - adoption_rates.push_back({InfectionState::I, InfectionState::R, mio::regions::Region(r), 0.99 / 5., {}}); - adoption_rates.push_back({InfectionState::I, InfectionState::D, mio::regions::Region(r), 0.01 / 5., {}}); + adoption_rates.push_back({{InfectionState::E}, {InfectionState::C}, mio::regions::Region(r), 1.0 / 5., {}}); + adoption_rates.push_back({{InfectionState::C}, {InfectionState::R}, mio::regions::Region(r), 0.2 / 3., {}}); + adoption_rates.push_back({{InfectionState::C}, {InfectionState::I}, mio::regions::Region(r), 0.8 / 3., {}}); + adoption_rates.push_back({{InfectionState::I}, {InfectionState::R}, mio::regions::Region(r), 0.99 / 5., {}}); + adoption_rates.push_back({{InfectionState::I}, {InfectionState::D}, mio::regions::Region(r), 0.01 / 5., {}}); } //Set second-order adoption rate different for the two regions // adoption rate has the form {i, j, k, c_{i,j}, {{tau1.state, tau1.factor}, {tau2.state, tau2.factor}}}, see the equation below - adoption_rates.push_back({InfectionState::S, InfectionState::E, mio::regions::Region(0), 0.1, {{InfectionState::C, 1}, {InfectionState::I, 0.5}}}); - adoption_rates.push_back({InfectionState::S, InfectionState::E, mio::regions::Region(1), 0.2, {{InfectionState::C, 1}, {InfectionState::I, 0.5}}}); + adoption_rates.push_back({{InfectionState::S}, {InfectionState::E}, mio::regions::Region(0), 0.1, {{{InfectionState::C}, 1}, {{InfectionState::I}, 0.5}}}); + adoption_rates.push_back({{InfectionState::S}, {InfectionState::E}, mio::regions::Region(1), 0.2, {{{InfectionState::C}, 1}, {{InfectionState::I}, 0.5}}}); //Initialize model parameter - model.parameters.get>() = adoption_rates; - -Sociodemographic stratification -------------------------------- - -Sociodemographic stratification e.g. by age, gender or immunity can be incorporated by stratifying the set of infection states passed as template to the model. + model.parameters.get>() = adoption_rates; Parameters ---------- @@ -110,8 +123,9 @@ with :math:`i^{(k)}` the population in region :math:`k` having infection state : Initial conditions ------------------ -Before running a simulation with the stochastic metapopulation model, the initial populations i.e. the number of agents per infection state for every region have to be set. -These populations have the class type **Populations** and can be set via: +Before running a simulation with the stochastic metapopulation model, the initial populations i.e. the number of agents +per infection state for every region have to be set. +These populations have the class type ``Populations`` and can be set via: .. code-block:: cpp @@ -128,7 +142,8 @@ These populations have the class type **Populations** and can be set via: } If individuals should transition between regions, the spatial transition rates of the model have to be initialized as well. -As the spatial transition rates are dependent on infection state, region changes for specific infection states can be prevented. Below, symmetric spatial transition rates are set for every region: +As the spatial transition rates are dependent on infection state, region changes for specific infection states can be +prevented. Below, symmetric spatial transition rates are set for every region: .. code-block:: cpp @@ -140,27 +155,32 @@ As the spatial transition rates are dependent on infection state, region changes if (i != j) { // transition rate has the form {i, k, l, \lambda^{(k,l)}_{i}} transition_rates.push_back( - {InfectionState(s), mio::regions::Region(i), mio::regions::Region(j), 0.01}); + {{InfectionState(s)}, mio::regions::Region(i), mio::regions::Region(j), 0.01}); transition_rates.push_back( - {InfectionState(s), mio::regions::Region(j), mio::regions::Region(i), 0.01}); + {{InfectionState(s)}, mio::regions::Region(j), mio::regions::Region(i), 0.01}); } } } //Initialize model parameter - model.parameters.get>() = transition_rates; + model.parameters.get>() = transition_rates; Nonpharmaceutical interventions -------------------------------- -There are no nonpharmaceutical interventions (NPIs) explicitly implemented in the model. However, NPIs influencing the adoption or spatial transition rates can be realized by resetting the corresponding model parameters. +There are no nonpharmaceutical interventions (NPIs) explicitly implemented in the model. However, NPIs influencing the +adoption or spatial transition rates can be realized by resetting the corresponding model parameters. Simulation ---------- -At the beginning of the simulation, the waiting times for all events (infection state adoptions and spatial transitions) are drawn. Then the time is advanced until the time point of the next event - which can be a spatial transition or an infection state adoption - and the event takes places. The waiting times of the other events are updated and a new waiting time for the event that just happened is drawn. The simulation saves the system state in discrete time steps. +At the beginning of the simulation, the waiting times for all events (infection state adoptions and spatial transitions) +are drawn. Then the time is advanced until the time point of the next event - which can be a spatial transition or an +infection state adoption - and the event takes places. The waiting times of the other events are updated and a new waiting +time for the event that just happened is drawn. The simulation saves the system state in discrete time steps. -To simulate the model from `t0` to `tmax` with given step size `dt`, a **Simulation** has to be created and advanced until `tmax`. The step size is only used to regularly save the system state during the simulation. +To simulate the model from ``t0`` to ``tmax`` with given step size ``dt``, a **Simulation** has to be created and advanced +until ``tmax``. The step size is only used to regularly save the system state during the simulation. .. code-block:: cpp @@ -193,10 +213,58 @@ If one wants to interpolate the aggregated results to a ``mio::TimeSeries`` cont auto interpolated_results = mio::interpolate_simulation_result(sim.get_result()); + +Demographic Stratification +-------------------------- + +It is possible to add multiple indices to the ``Status`` to differentiate multiple groups on the same region. For example this +could represent the human and the mosquito populations in a specific region. To use this feature, one first of all has to +create a new index: + +.. code-block:: cpp + + struct Species : public mio::Index { + Species(size_t val): Index(val) + { + } + }; + +Then, one has to create the multiindex, where we reuse the ``InfectionState`` defined in the first example: + +.. code-block:: cpp + + using Status = mio::Index; + +Define the size for each index dimension that is not an enum class: + +.. code-block:: cpp + + const size_t num_species = 2; + +We can define a model: + +.. code-block:: cpp + + using Model = mio::smm::Model + Model model(Status{Count, Species(num_species)}, Region(num_regions)); + +Now, for accessing the population, all indexes need to be given: + +-- code-block:: cpp + + model.populations[{Region(r), InfectionState::S, Species(0)}] = 100; + // ... + adoption_rates.push_back({{InfectionState::S, Species(0)}, {InfectionState::E, Species(0)}, Region(r), 0.1, {{InfectionState::I, Species(1)}, 1}); + // ... + transition_rates.push_back({{InfectionState::S, Species(0)}, Region(i), Region(j), 0.01}); + + Examples -------- -An example of the stochastic metapopulation model with four regions can be found at: `examples/smm.cpp `_ +A full example with ``Status`` distributed according to ``InfectionState``, ``Age`` and ``Species`` can be found at +`examples/smm.cpp `_ + Overview of the ``smm`` namespace: