1341 smm model with more indices

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 <[email protected]>
 *
@@ -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"
@@ -33,61 +34,80 @@ 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;
-    using Model              = mio::smm::Model<ScalarType, num_regions, InfectionState>;
+    const size_t num_regions    = 4;
+    const size_t num_age_groups = 1;
+    const size_t num_groups     = 1;
+    using Model                 = mio::smm::Model<ScalarType, num_regions, InfectionState, num_age_groups, num_groups>;
 
     ScalarType 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, Age(0), Species(0)}] =
             (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, 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
-    std::vector<mio::AdoptionRate<ScalarType, InfectionState>> adoption_rates;
-    std::vector<mio::smm::TransitionRate<ScalarType, InfectionState>> transition_rates;
+    std::vector<mio::AdoptionRate<ScalarType, InfectionState, Age, Species>> adoption_rates;
+    std::vector<mio::smm::TransitionRate<ScalarType, InfectionState, Age, Species>> 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}, {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, 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)}});
     }
 
     //Agents in infection state D do not transition
     for (size_t s = 0; s < static_cast<size_t>(InfectionState::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),
+                                                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)}});
                 }
         }
     }
 
-    model.parameters.get<mio::smm::AdoptionRates<ScalarType, InfectionState>>()   = adoption_rates;
-    model.parameters.get<mio::smm::TransitionRates<ScalarType, InfectionState>>() = transition_rates;
+    model.parameters.get<mio::smm::AdoptionRates<ScalarType, InfectionState, Age, Species>>()   = adoption_rates;
+    model.parameters.get<mio::smm::TransitionRates<ScalarType, InfectionState, Age, Species>>() = transition_rates;
 
     ScalarType dt   = 0.1;
     ScalarType tmax = 30.0;

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 <[email protected]>
 *
@@ -23,35 +23,45 @@
 #include "memilio/utils/index.h"
 #include "memilio/config.h"
 #include "memilio/geography/regions.h"
+#include <limits>
+#include <tuple>
+#include <optional>
 
 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 <typename FP, class Status>
-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".
  * 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 <typename FP, class Status>
+template <typename FP, class Status, class... Groups>
 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<mio::regions::Region> region = std::nullopt;
+        std::tuple<Groups...> group_indices{};
+    };
+
     Status from; // i
     Status to; // j
     mio::regions::Region region; // k
     FP factor; // gammahat_{ij}^k
-    std::vector<Influence<FP, Status>> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} )
+    std::vector<Influence> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} )
+    std::tuple<Groups...> group_indices{};
 };
 
 } // namespace mio

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 <[email protected]>
 *
@@ -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,21 +33,29 @@ namespace mio
 namespace smm
 {
 
+template <class T, std::size_t>
+using age_group = T;
+
 /**
  * @brief Stochastic Metapopulation Model.
  * @tparam regions Number of regions.
  * @tparam Status An infection state enum.
  */
-template <typename FP, size_t regions, class Status>
-class Model : public mio::CompartmentalModel<FP, Status, mio::Populations<FP, mio::regions::Region, Status>,
-                                             ParametersBase<FP, Status>>
+template <typename FP, size_t regions, class Status, size_t... Groups>
+class Model : public mio::CompartmentalModel<
+                  FP, Status, mio::Populations<FP, mio::regions::Region, Status, age_group<mio::AgeGroup, Groups>...>,
+                  ParametersBase<FP, Status, age_group<mio::AgeGroup, Groups>...>>
 {
-    using Base = mio::CompartmentalModel<FP, Status, mio::Populations<FP, mio::regions::Region, Status>,
-                                         ParametersBase<FP, Status>>;
+    using Base =
+        mio::CompartmentalModel<FP, Status,
+                                mio::Populations<FP, mio::regions::Region, Status, age_group<mio::AgeGroup, Groups>...>,
+                                ParametersBase<FP, Status, age_group<mio::AgeGroup, Groups>...>>;
+    using Index = mio::Index<mio::regions::Region, Status, age_group<mio::AgeGroup, Groups>...>;
 
 public:
     Model()
-        : Base(typename Base::Populations({static_cast<mio::regions::Region>(regions), Status::Count}, 0.0),
+        : Base(typename Base::Populations(
+                   {static_cast<mio::regions::Region>(regions), Status::Count, static_cast<mio::AgeGroup>(Groups)...}),
                typename Base::ParameterSet())
     {
     }
@@ -57,26 +66,46 @@ class Model : public mio::CompartmentalModel<FP, Status, mio::Populations<FP, mi
      * @param[in] x The current state of the model.
      * @return Current value of the adoption rate.
      */
-    FP evaluate(const AdoptionRate<FP, Status>& rate, const Eigen::VectorX<FP>& x) const
+    FP evaluate(const AdoptionRate<FP, Status, age_group<mio::AgeGroup, Groups>...>& 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 index_from = std::apply(
+            [&](auto&&... args) {
+                return Index{rate.region, rate.from, std::forward<decltype(args)>(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
-            FP N = 0;
-            for (size_t s = 0; s < static_cast<size_t>(Status::Count); ++s) {
-                N += x[pop.get_flat_index({rate.region, Status(s)})];
-            }
             // accumulate influences
             FP 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})];
+                FP N = 0; // Welches N brauchen wir hier??
+
+                for (size_t s = 0; s < static_cast<size_t>(Status::Count); ++s) {
+                    const auto index = std::apply(
+                        [&](auto&&... args) {
+                            return Index{rate.influences[i].region.value_or(rate.region), Status(s),
+                                         std::forward<decltype(args)>(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<decltype(args)>(args)...};
+                    },
+                    rate.influences[i].group_indices);
+                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;
         }
     }
 
@@ -86,9 +115,15 @@ class Model : public mio::CompartmentalModel<FP, Status, mio::Populations<FP, mi
      * @param[in] x The current state of the model.
      * @return Current value of the transition rate.
      */
-    FP evaluate(const TransitionRate<FP, Status>& rate, const Eigen::VectorX<FP>& x) const
+    FP evaluate(const TransitionRate<FP, Status, age_group<mio::AgeGroup, Groups>...>& rate,
+                const Eigen::VectorXd& x) const
     {
-        const auto source = this->populations.get_flat_index({rate.from, rate.status});
+        auto index = std::apply(
+            [&](auto&&... args) {
+                return Index{rate.from, rate.status, std::forward<decltype(args)>(args)...};
+            },
+            rate.group_indices_from);
+        const auto source = this->populations.get_flat_index(index);
         return rate.factor * x[source];
     }
 

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 <tuple>
 
 namespace mio
 {
@@ -35,9 +36,9 @@ namespace smm
  * @brief A vector of AdoptionRate%s, see mio::AdoptionRate
  * @tparam Status An infection state enum.
  */
-template <typename FP, class Status>
+template <typename FP, class Status, class... Groups>
 struct AdoptionRates {
-    using Type = std::vector<AdoptionRate<FP, Status>>;
+    using Type = std::vector<AdoptionRate<FP, Status, Groups...>>;
     const static std::string name()
     {
         return "AdoptionRates";
@@ -48,25 +49,26 @@ struct AdoptionRates {
  * @brief Struct defining a possible regional transition in a Model based on Poisson Processes.
  * @tparam Status An infection state enum.
  */
-template <typename FP, class Status>
+template <typename FP, class Status, class... Groups>
 struct TransitionRate {
     Status status; // i
     mio::regions::Region from; // k
     mio::regions::Region to; // l
     FP factor; // lambda_i^{kl}
+    std::tuple<Groups...> group_indices_from{};
+    std::tuple<Groups...> group_indices_to{};
 };
-
-template <typename FP, class Status>
+template <typename FP, class Status, class... Groups>
 struct TransitionRates {
-    using Type = std::vector<TransitionRate<FP, Status>>;
+    using Type = std::vector<TransitionRate<FP, Status, Groups...>>;
     const static std::string name()
     {
         return "TransitionRates";
     }
 };
 
-template <typename FP, class Status>
-using ParametersBase = mio::ParameterSet<AdoptionRates<FP, Status>, TransitionRates<FP, Status>>;
+template <typename FP, class Status, class... Groups>
+using ParametersBase = mio::ParameterSet<AdoptionRates<FP, Status, Groups...>, TransitionRates<FP, Status, Groups...>>;
 
 } // namespace smm