1341 SMM Model with generic indices

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 <Martin.Kuehn@DLR.de>
 *
@@ -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,75 @@ enum class InfectionState
     R,
     D,
     Count
+};
 
+struct Species : public mio::Index<Species> {
+    Species(size_t val)
+        : Index<Species>(val)
+    {
+    }
 };
 
 int main()
 {
+    using Age    = mio::AgeGroup;
+    using Status = mio::Index<InfectionState, Age, Species>;
+    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<ScalarType, num_regions, InfectionState>;
+    /* Example how to run the stochastic metapopulation models with four regions. Within each region we differentiate by
+       age groups, species and infection states. The infection states are S, E, C, I, R, D. For the number of age groups
+       and species we choose: */
+    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, InfectionState, Status, Region>;
 
     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;
     }
 
-    //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;
+    using AR = mio::smm::AdoptionRates<ScalarType, Status, Region>;
+    using TR = mio::smm::TransitionRates<ScalarType, Status, Region>;
+
+    //Set infection state adoption rates. Adoptions only happen within a region.
+    AR::Type adoption_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<size_t>(InfectionState::D); ++s) {
+    //Set transition rates between regions. Agents in infection state D do not transition.
+    TR::Type transition_rates;
+    for (size_t s = 0; s < static_cast<size_t>(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<mio::smm::AdoptionRates<ScalarType, InfectionState>>()   = adoption_rates;
-    model.parameters.get<mio::smm::TransitionRates<ScalarType, InfectionState>>() = transition_rates;
+    model.parameters.get<AR>() = adoption_rates;
+    model.parameters.get<TR>() = 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 <Martin.Kuehn@DLR.de>
 *
@@ -20,17 +20,16 @@
 #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 <vector>
 
 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 <typename FP, class Status>
 struct Influence {
@@ -43,13 +42,14 @@ 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 <typename FP, class Status>
+template <typename FP, class Status, class Region = mio::regions::Region>
 struct AdoptionRate {
     Status from; // i
     Status to; // j
-    mio::regions::Region region; // k
+    Region region; // k
     FP factor; // gammahat_{ij}^k
     std::vector<Influence<FP, Status>> influences; // influences[tau] = ( Psi_{i,j,tau} , gamma_{i,j,tau} )
 };

cpp/memilio/epidemiology/populations.h

@@ -287,6 +287,18 @@ class Populations : public CustomIndexArray<UncertainValue<FP>, 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 <typename FP, class... Categories>
+class Populations<FP, Index<Categories...>> : public Populations<FP, Categories...>
+{
+    using Populations<FP, Categories...>::Populations;
+};
+
 } // namespace mio
 
 #endif // MIO_EPI_POPULATIONS_H