diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index 570ce17c6b..89dd7185b3 100644 --- a/cpp/examples/smm.cpp +++ b/cpp/examples/smm.cpp @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2025 MEmilio * -* Authors: Julia Bicker, René Schmieding +* Authors: Julia Bicker, René Schmieding, Kilian Volmer * * Contact: Martin J. Kuehn * @@ -18,11 +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 { @@ -33,61 +34,73 @@ enum class InfectionState R, D, Count +}; +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; - using Model = mio::smm::Model; + const size_t num_regions = 4; + const size_t num_age_groups = 1; + 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}] = - (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[{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}, {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., {}}); + {{{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}); - transition_rates.push_back( - {InfectionState(s), mio::regions::Region(j), mio::regions::Region(i), 0.01}); + 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 705f729025..8f72d650d2 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 * @@ -20,9 +20,8 @@ #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 namespace mio { @@ -30,7 +29,7 @@ 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 Status An infection state enum or MultiIndex. */ template struct Influence { @@ -43,15 +42,16 @@ struct Influence { * The AdoptionRate is considered to be of second-order if there are any "influences". * 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 Status An infection state enum or MultiIndex. + * @tparam Region A MultiIndex. */ -template +template struct AdoptionRate { 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::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..b702333ad9 100644 --- a/cpp/memilio/epidemiology/populations.h +++ b/cpp/memilio/epidemiology/populations.h @@ -287,6 +287,18 @@ class Populations : public CustomIndexArray, Categories...> } }; +/** + * @brief Population template specialization, forwarding categories from a MultiIndex to the Population. + * + * @tparam FP A floating point type, e.g., double. + * @tparam Categories Index 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..dcc84d267d 100644 --- a/cpp/memilio/utils/index.h +++ b/cpp/memilio/utils/index.h @@ -20,15 +20,75 @@ #ifndef INDEX_H #define INDEX_H +#include "memilio/io/io.h" #include "memilio/utils/compiler_diagnostics.h" +#include "memilio/utils/metaprogramming.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 with 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 single-category 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 single-category index from an enum value. +template +std::tuple> get_tuple(Enum i) + requires std::is_enum::value +{ + return std::tuple(Index(i)); +} + +/// @brief Merge a series of enums or MultIndex%s into a tuple of single-category indices. +template + requires((std::is_enum_v || IsMultiIndex) && ...) +decltype(auto) concatenate_indices_impl(IndexArgs&&... args) +{ + return std::tuple_cat(details::get_tuple(args)...); +} + +/** + * @brief Function declaration that allows type conversion from a tuple of single-category indices to MultiIndex. + * Used together with concatenate_indices_impl, this allows combining categories of multiple args into a single MultiIndex. + */ +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 @@ -62,9 +122,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); } @@ -72,8 +133,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) { } @@ -119,20 +180,31 @@ 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; + /// @brief Construct an Index filled with zeroes. static Index constexpr Zero() { return Index(Index::Zero()...); } - // constructor from Indices + /// @brief Constructor from individual Indices. Index(Index const&... _indices) : indices{_indices...} { } + /// @brief Constructor from mixed Indices and MultiIndices. + template + requires(sizeof...(IndexArgs) > 1) + Index(IndexArgs&&... subindices) + : indices(details::concatenate_indices_impl(std::forward(subindices)...)) + { + } + private: + /// @brief Internal constructor from a tuple. Index(const std::tuple...>& _indices) : indices(_indices) { @@ -150,6 +222,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 @@ -197,70 +288,75 @@ 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> +/// @brief Retrieves the Index (by reference) at the Ith position of a MultiIndex. +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> +/// @brief Retrieves the Index (by const reference) at the Ith position of a MultiIndex. +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; -} - -// retrieves the Index for the tag Tag of a Index with more than one Tag -template 1), void>* = nullptr> -constexpr Index& get(Index& 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 for the tag Tag of a Index with one Tag, equals identity function -template * = nullptr> +/// @brief Retrieves the Index (by reference) by its Tag in a MultiIndex. Requires unique tags. +template 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> +/// @brief Retrieves the Index (by const reference) by its Tag in a MultiIndex. Requires unique tags. +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 +/** + * @brief Combine several Index%s into one MultiIndex. + * @param args Either enum or MultiIndex values. + * @return A MultiIndex with all categories and values of the given Indices concatonated. + */ +template +decltype(auto) concatenate_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::concatenate_indices_impl(std::declval()...))); + return MergedIndex(std::forward(args)...); } namespace details @@ -310,12 +406,27 @@ inline Index extend_index_impl(const Index& i, const * @tparam SubIndex An Index that contains a subset of the categories from SuperIndex. * @tparam SuperIndex Any Index. * @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) +{ + 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 + requires std::is_enum_v +Index reduce_index(const SuperIndex& index) { - return details::reduce_index_impl(index, mio::Tag{}); + return details::reduce_index_impl(index, mio::Tag>{}); } +/** @} */ /** * @brief Create a SuperIndex by copying values from SubIndex, filling new categories with fill_value. 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 56d9b56251..056878cfb7 100644 --- a/cpp/models/smm/model.h +++ b/cpp/models/smm/model.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 * @@ -21,32 +21,46 @@ #ifndef MIO_SMM_MODEL_H #define MIO_SMM_MODEL_H -#include "memilio/config.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 { namespace smm { +/// @brief The Index type used to define the SMM population. +template +using PopulationIndex = decltype(concatenate_indices(std::declval(), std::declval())); + /** * @brief Stochastic Metapopulation Model. - * @tparam regions Number of regions. - * @tparam Status An infection state enum. + * The stratification of the population of this model is split between "Status" and "Region". This split is mostly + * arbitrary, with the important distinction, that for second order rates the reference population + * (i.e., the N in S' = S * I / N) is calculated by accumulating subpopulations only over the Status. + * @tparam Comp An enum representing the infection states. Must also be contained in Status + * @tparam Status A MultiIndex allowing to further stratify infection state adoptions. + * @tparam Region A MultiIndex for "spatially" separate subpopulations, default is @ref mio::regions::Region. */ -template -class Model : public mio::CompartmentalModel, - ParametersBase> +template +class Model : public mio::CompartmentalModel>, + ParametersBase> { - using Base = mio::CompartmentalModel, - ParametersBase>; + using Base = mio::CompartmentalModel>, + ParametersBase>; + static_assert(!Base::Populations::Index::has_duplicates, "Do not use the same Index tag multiple times!"); public: - Model() - : Base(typename Base::Populations({static_cast(regions), Status::Count}, 0.0), + using Status = StatusT; + using Region = RegionT; + + Model(Status status_dimensions, Region region_dimensions) + : Base(typename Base::Populations(concatenate_indices(region_dimensions, status_dimensions)), typename Base::ParameterSet()) { } @@ -57,7 +71,7 @@ class Model : public mio::CompartmentalModel& rate, const Eigen::VectorX& x) const + FP 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}); @@ -67,8 +81,8 @@ class Model : public mio::CompartmentalModel(Status::Count); ++s) { - N += x[pop.get_flat_index({rate.region, Status(s)})]; + 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; @@ -86,7 +100,7 @@ class Model : public mio::CompartmentalModel& rate, const Eigen::VectorX& x) const + FP 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]; diff --git a/cpp/models/smm/parameters.h b/cpp/models/smm/parameters.h index d0a6011f6e..652d8891d8 100644 --- a/cpp/models/smm/parameters.h +++ b/cpp/models/smm/parameters.h @@ -33,11 +33,13 @@ namespace smm /** * @brief A vector of AdoptionRate%s, see mio::AdoptionRate - * @tparam Status An infection state enum. + * @tparam FP A floating point type, e.g., double. + * @tparam Status A MultiIndex, containing the infection state enum. + * @tparam Region A MultiIndex for spatial stratification. */ -template +template struct AdoptionRates { - using Type = std::vector>; + using Type = std::vector>; const static std::string name() { return "AdoptionRates"; @@ -46,27 +48,35 @@ struct AdoptionRates { /** * @brief Struct defining a possible regional transition in a Model based on Poisson Processes. - * @tparam Status An infection state enum. + * @tparam FP A floating point type, e.g., double. + * @tparam Status A MultiIndex, containing the infection state enum. + * @tparam Region A MultiIndex for spatial stratification. */ -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} }; -template +/** + * @brief A vector of TransitionRate%s, see mio::TransitionRate + * @tparam FP A floating point type, e.g., double. + * @tparam Status A MultiIndex, containing the infection state enum. + * @tparam Region A MultiIndex for spatial stratification. + */ +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 49d4dcbb8f..0ed56ee13c 100644 --- a/cpp/models/smm/simulation.h +++ b/cpp/models/smm/simulation.h @@ -21,10 +21,9 @@ #ifndef MIO_SMM_SIMULATION_H #define MIO_SMM_SIMULATION_H -#include "memilio/config.h" #include "smm/model.h" #include "smm/parameters.h" -#include "memilio/compartments/simulation.h" +#include "memilio/utils/time_series.h" namespace mio { @@ -37,12 +36,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. @@ -61,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); @@ -164,17 +164,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>(); } /** @@ -217,7 +217,7 @@ class Simulation std::vector m_current_rates; ///< Current values of both types of rates i.e. adoption and transition rates. }; -} //namespace smm +} // namespace smm } // namespace mio #endif diff --git a/cpp/tests/test_smm_model.cpp b/cpp/tests/test_smm_model.cpp index b671cfa2d0..0e906dc13f 100644 --- a/cpp/tests/test_smm_model.cpp +++ b/cpp/tests/test_smm_model.cpp @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2025 MEmilio * -* Authors: Julia Bicker +* Authors: Julia Bicker, Kilian Volmer * * Contact: Martin J. Kuehn * @@ -18,18 +18,17 @@ * limitations under the License. */ -#include -#include -#include -#include -#include -#include "memilio/utils/compiler_diagnostics.h" +#include "abm_helpers.h" +#include "memilio/epidemiology/adoption_rate.h" +#include "memilio/epidemiology/age_group.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 @@ -45,6 +44,14 @@ enum class InfectionState }; +enum class Infections +{ + S, + I, + R, + Count +}; + TEST(TestSMM, evaluateAdoptionRate) { //Test whether the adoption rates are evaluated correctly. @@ -53,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; @@ -80,13 +87,52 @@ TEST(TestSMM, evaluateAdoptionRate) EXPECT_EQ(model.evaluate(adoption_rates[1], model.populations.get_compartments()), 2.); } +TEST(TestSMM, evaluateinterregionalAdoptionRate) +{ + using Age = mio::AgeGroup; + using mio::regions::Region; + using Status = mio::Index; + using Model = mio::smm::Model>; + + std::vector>> adoption_rates; + //Second-order adoption + adoption_rates.push_back({{Infections::S, Age(0), Region(0)}, + {Infections::I, Age(0), 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)}, {}, 0.2, {}}); + Model model({Infections::Count, Age(2), Region(2)}, {}); + + 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[{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[{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[{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. + 0.2 * 16 * 1.) / (60 + 120 + 120 + 96)); + EXPECT_DOUBLE_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; + 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}] = 50; model.populations[{mio::regions::Region(0), InfectionState::E}] = 10; @@ -114,9 +160,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; @@ -151,8 +197,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>>> @@ -195,9 +241,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; @@ -211,8 +257,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 @@ -228,7 +274,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; @@ -243,7 +289,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; @@ -252,13 +298,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; 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: