896 add age resolved lct model (#1120)

- Added the possibility to use groups (e.g. age groups) in the LCT-SECIR model. The parameters and the number of subcompartments can be set differently for each group.
- Adapted functionality to initialize the model.
- Added tests.

Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com>
Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com>

Author: 3 people

cpp/examples/lct_secir.cpp

@@ -19,10 +19,12 @@
 */
 
 #include "lct_secir/model.h"
+#include "lct_secir/infection_state.h"
 #include "lct_secir/initializer_flows.h"
 #include "memilio/config.h"
 #include "memilio/utils/time_series.h"
 #include "memilio/epidemiology/uncertain_matrix.h"
+#include "memilio/epidemiology/lct_infection_state.h"
 #include "memilio/math/eigen.h"
 #include "memilio/utils/logging.h"
 #include "memilio/compartments/simulation.h"
@@ -33,14 +35,14 @@
 int main()
 {
     // Simple example to demonstrate how to run a simulation using an LCT-SECIR model.
+    // One single AgeGroup/Category member is used here.
     // Parameters, initial values and the number of subcompartments are not meant to represent a realistic scenario.
     constexpr size_t NumExposed = 2, NumInfectedNoSymptoms = 3, NumInfectedSymptoms = 1, NumInfectedSevere = 1,
                      NumInfectedCritical = 5;
-    using Model          = mio::lsecir::Model<NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms, NumInfectedSevere,
-                                     NumInfectedCritical>;
-    using LctState       = Model::LctState;
-    using InfectionState = LctState::InfectionState;
-
+    using InfState                       = mio::lsecir::InfectionState;
+    using LctState = mio::LctInfectionState<InfState, 1, NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms,
+                                            NumInfectedSevere, NumInfectedCritical, 1, 1>;
+    using Model    = mio::lsecir::Model<LctState>;
     Model model;
 
     // Variable defines whether the class Initializer is used to define an initial vector from flows or whether a manually
@@ -50,34 +52,33 @@ int main()
     ScalarType tmax = 20;
 
     // Set Parameters.
-    model.parameters.get<mio::lsecir::TimeExposed>()            = 3.2;
-    model.parameters.get<mio::lsecir::TimeInfectedNoSymptoms>() = 2.;
-    model.parameters.get<mio::lsecir::TimeInfectedSymptoms>()   = 5.8;
-    model.parameters.get<mio::lsecir::TimeInfectedSevere>()     = 9.5;
-    // It is also possible to change values with the set function.
-    model.parameters.set<mio::lsecir::TimeInfectedCritical>(7.1);
+    model.parameters.get<mio::lsecir::TimeExposed>()[0]            = 3.2;
+    model.parameters.get<mio::lsecir::TimeInfectedNoSymptoms>()[0] = 2.;
+    model.parameters.get<mio::lsecir::TimeInfectedSymptoms>()[0]   = 5.8;
+    model.parameters.get<mio::lsecir::TimeInfectedSevere>()[0]     = 9.5;
+    model.parameters.get<mio::lsecir::TimeInfectedCritical>()[0]   = 7.1;
 
-    model.parameters.get<mio::lsecir::TransmissionProbabilityOnContact>() = 0.05;
+    model.parameters.get<mio::lsecir::TransmissionProbabilityOnContact>()[0] = 0.05;
 
     mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::lsecir::ContactPatterns>();
     contact_matrix[0]                       = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10));
     // From SimulationTime 5, the contact pattern is reduced to 30% of the initial value.
     contact_matrix[0].add_damping(0.7, mio::SimulationTime(5.));
 
-    model.parameters.get<mio::lsecir::RelativeTransmissionNoSymptoms>() = 0.7;
-    model.parameters.get<mio::lsecir::RiskOfInfectionFromSymptomatic>() = 0.25;
-    model.parameters.get<mio::lsecir::RecoveredPerInfectedNoSymptoms>() = 0.09;
-    model.parameters.get<mio::lsecir::SeverePerInfectedSymptoms>()      = 0.2;
-    model.parameters.get<mio::lsecir::CriticalPerSevere>()              = 0.25;
-    model.parameters.set<mio::lsecir::DeathsPerCritical>(0.3);
+    model.parameters.get<mio::lsecir::RelativeTransmissionNoSymptoms>()[0] = 0.7;
+    model.parameters.get<mio::lsecir::RiskOfInfectionFromSymptomatic>()[0] = 0.25;
+    model.parameters.get<mio::lsecir::RecoveredPerInfectedNoSymptoms>()[0] = 0.09;
+    model.parameters.get<mio::lsecir::SeverePerInfectedSymptoms>()[0]      = 0.2;
+    model.parameters.get<mio::lsecir::CriticalPerSevere>()[0]              = 0.25;
+    model.parameters.get<mio::lsecir::DeathsPerCritical>()[0]              = 0.3;
 
     if (use_initializer_flows) {
         // Example how to use the class Initializer for the definition of an initial vector for the LCT model.
 
-        ScalarType dt                    = 0.001;
-        ScalarType total_population      = 1000000;
-        ScalarType deaths                = 10;
-        ScalarType total_confirmed_cases = 16000;
+        ScalarType dt                                 = 0.001;
+        mio::Vector<ScalarType> total_population      = mio::Vector<ScalarType>::Constant(1, 1000000.);
+        mio::Vector<ScalarType> deaths                = mio::Vector<ScalarType>::Constant(1, 10.);
+        mio::Vector<ScalarType> total_confirmed_cases = mio::Vector<ScalarType>::Constant(1, 16000.);
 
         // Create TimeSeries with num_transitions elements.
         int num_transitions = (int)mio::lsecir::InfectionTransition::Count;
@@ -120,22 +121,22 @@ int main()
                                                                     {50},  {10, 10, 5, 3, 2}, {20},         {10}};
 
         // Assert that initial_populations has the right shape.
-        if (initial_populations.size() != (size_t)InfectionState::Count) {
+        if (initial_populations.size() != (size_t)InfState::Count) {
             mio::log_error(
                 "The number of vectors in initial_populations does not match the number of InfectionStates.");
             return 1;
         }
-        if ((initial_populations[(size_t)InfectionState::Susceptible].size() !=
-             LctState::get_num_subcompartments<InfectionState::Susceptible>()) ||
-            (initial_populations[(size_t)InfectionState::Exposed].size() != NumExposed) ||
-            (initial_populations[(size_t)InfectionState::InfectedNoSymptoms].size() != NumInfectedNoSymptoms) ||
-            (initial_populations[(size_t)InfectionState::InfectedSymptoms].size() != NumInfectedSymptoms) ||
-            (initial_populations[(size_t)InfectionState::InfectedSevere].size() != NumInfectedSevere) ||
-            (initial_populations[(size_t)InfectionState::InfectedCritical].size() != NumInfectedCritical) ||
-            (initial_populations[(size_t)InfectionState::Recovered].size() !=
-             LctState::get_num_subcompartments<InfectionState::Recovered>()) ||
-            (initial_populations[(size_t)InfectionState::Dead].size() !=
-             LctState::get_num_subcompartments<InfectionState::Dead>())) {
+        if ((initial_populations[(size_t)InfState::Susceptible].size() !=
+             LctState::get_num_subcompartments<InfState::Susceptible>()) ||
+            (initial_populations[(size_t)InfState::Exposed].size() != NumExposed) ||
+            (initial_populations[(size_t)InfState::InfectedNoSymptoms].size() != NumInfectedNoSymptoms) ||
+            (initial_populations[(size_t)InfState::InfectedSymptoms].size() != NumInfectedSymptoms) ||
+            (initial_populations[(size_t)InfState::InfectedSevere].size() != NumInfectedSevere) ||
+            (initial_populations[(size_t)InfState::InfectedCritical].size() != NumInfectedCritical) ||
+            (initial_populations[(size_t)InfState::Recovered].size() !=
+             LctState::get_num_subcompartments<InfState::Recovered>()) ||
+            (initial_populations[(size_t)InfState::Dead].size() !=
+             LctState::get_num_subcompartments<InfState::Dead>())) {
             mio::log_error(
                 "The length of at least one vector in initial_populations does not match the related number of "
                 "subcompartments.");
@@ -148,15 +149,16 @@ int main()
             flat_initial_populations.insert(flat_initial_populations.end(), vec.begin(), vec.end());
         }
         for (size_t i = 0; i < LctState::Count; i++) {
-            model.populations[mio::Index<LctState>(i)] = flat_initial_populations[i];
+            model.populations[i] = flat_initial_populations[i];
         }
     }
 
     // Perform a simulation.
     mio::TimeSeries<ScalarType> result = mio::simulate<ScalarType, Model>(0, tmax, 0.5, model);
+    //result.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
     // The simulation result is divided by subcompartments.
     // We call the function calculate_comparttments to get a result according to the InfectionStates.
     mio::TimeSeries<ScalarType> population_no_subcompartments = model.calculate_compartments(result);
     auto interpolated_results = mio::interpolate_simulation_result(population_no_subcompartments);
     interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
-}
+}

cpp/memilio/CMakeLists.txt

@@ -17,6 +17,7 @@ add_library(memilio
     epidemiology/dynamic_npis.h
     epidemiology/dynamic_npis.cpp
     epidemiology/lct_infection_state.h
+    epidemiology/lct_populations.h
     geography/regions.h
     geography/regions.cpp
     epidemiology/simulation_day.h

cpp/memilio/epidemiology/lct_infection_state.h

@@ -26,10 +26,13 @@ namespace mio
 {
 /**
  * @brief Provides the functionality to be able to work with subcompartments in an LCT model.
+
+ * This class just stores the number of subcompartments for each InfectionState and not the number of individuals in
+ * each subcompartment. 
  *
  * @tparam InfectionStates An enum class that defines the basic infection states.
  * @tparam Ns Number of subcompartments for each infection state defined in InfectionState. 
- *      The number of given template arguments must be equal to the entry Count from InfectionState.
+ *      The number of given template arguments must be equal to the entry Count from InfectionStates.
  */
 template <class InfectionStates, size_t... Ns>
 class LctInfectionState
@@ -84,4 +87,4 @@ class LctInfectionState
 
 } // namespace mio
 
-#endif
+#endif

cpp/memilio/epidemiology/lct_populations.h

@@ -0,0 +1,216 @@
+/* 
+* Copyright (C) 2020-2024 MEmilio
+*
+* Authors: Lena Ploetzke
+*
+* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+#ifndef MIO_EPI_LCT_POPULATIONS_H
+#define MIO_EPI_LCT_POPULATIONS_H
+
+#include "boost/type_traits/make_void.hpp"
+#include "memilio/config.h"
+#include "memilio/utils/uncertain_value.h"
+#include "memilio/math/eigen.h"
+#include "memilio/epidemiology/lct_infection_state.h"
+#include "memilio/utils/type_list.h"
+#include "memilio/utils/metaprogramming.h"
+
+namespace mio
+{
+
+/**
+ * @brief A class template for compartment populations of LCT models.
+ *
+ * Populations can be split up into different categories, e.g. by age group, yearly income group, gender etc. 
+ * In LCT models, we want to use different numbers of subcompartments, i.e. different LctStates, 
+ * for each group of a category.
+ * (Therefore, we can't use the normal Populations class because it expects the same InfectionStates for each group.)
+ * 
+ * This template must be instantiated with one LctState for each group of a category. 
+ * The purpose of the LctStates is to define the number of subcompartments for each InfectionState.
+ * The number of LctStates also determines the number of groups. 
+ * If you want to use more than one category, e.g. age and gender, you have to define num_age_groups * num_genders 
+ * LctStates, because the number of subcompartments can be different 
+ * even for (female, A05-A14) and (female, A80+) or (male, A05-A14).
+ *
+ * The class created from this template contains a "flat array" of compartment populations and some functions for 
+ * retrieving or setting the populations. The order in the flat array is: First, all (sub-)compartments of the 
+ * first group, afterwards all (sub-)compartments of the second group and so on.
+ *
+ */
+
+template <typename FP = ScalarType, class... LctStates>
+class LctPopulations
+{
+public:
+    using Type                         = UncertainValue<FP>;
+    using InternalArrayType            = Eigen::Array<Type, Eigen::Dynamic, 1>;
+    using LctStatesGroups              = TypeList<LctStates...>;
+    static size_t constexpr num_groups = sizeof...(LctStates); ///< Number of groups.
+    static_assert(num_groups >= 1, "The number of LctStates provided should be at least one.");
+
+    /// @brief Default constructor.
+    LctPopulations()
+    {
+        set_count();
+        m_y = InternalArrayType::Constant(m_count, 1, 0.);
+    }
+
+    /**
+     * @brief get_num_compartments Returns the number of compartments.
+     * @return Number of compartments which equals the flat array size.
+     */
+    size_t get_num_compartments() const
+    {
+        return m_count;
+    }
+
+    /**
+     * @brief Returns a reference to the internally stored flat array.
+     * @return Const reference to the InternalArrayType instance.
+     */
+    auto const& array() const
+    {
+        return m_y;
+    }
+    auto& array()
+    {
+        return m_y;
+    }
+
+    /**
+     * @brief Returns the entry of the array given a flat index.
+     * @param index A flat index.
+     * @return The value of the internal array at the index.
+     */
+    Type& operator[](size_t index)
+    {
+        assert(index < m_count);
+        return m_y[index];
+    }
+
+    /**
+    * @brief Gets the first index of a group in the flat array.
+    * @tparam group The group for which the index should be returned.
+    * @return The index of the first entry of group in the flat array.
+    */
+    template <size_t Group>
+    size_t get_first_index_of_group() const
+    {
+        static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid.");
+        if constexpr (Group == 0) {
+            return 0;
+        }
+        else {
+            return get_first_index_of_group<Group - 1>() + type_at_index_t<Group - 1, LctStatesGroups>::Count;
+        }
+    }
+    /**
+     * @brief Returns an Eigen copy of the vector of populations. 
+     * This can be used as initial conditions for the ODE solver.
+     * @return Eigen::VectorXd of populations.
+     */
+    inline Vector<FP> get_compartments() const
+    {
+        return m_y.array().template cast<FP>();
+    }
+
+    /**
+     * @brief Returns the total population of a group.
+     * @tparam group The group for which the total population should be calculated.
+     * @return Total population of the group.
+     */
+    template <size_t Group>
+    FP get_group_total() const
+    {
+        return m_y.array()
+            .template cast<FP>()
+            .segment(get_first_index_of_group<Group>(), type_at_index_t<Group, LctStatesGroups>::Count)
+            .sum();
+    }
+
+    /**
+     * @brief Returns the total population of all compartments and groups.
+     * @return Total population.
+     */
+    FP get_total() const
+    {
+        return m_y.array().template cast<FP>().sum();
+    }
+
+    /**
+     * @brief Checks whether all compartments have non-negative values. 
+     * This function can be used to prevent slightly negative function values in compartment sizes that are produced 
+     * due to rounding errors if, e.g., population sizes were computed in a complex way.
+     *
+     * Attention: This function should be used with care. It can not and will not set model parameters and 
+     *            compartments to meaningful values. In most cases it is preferable to use check_constraints,
+     *            and correct values manually before proceeding with the simulation.
+     *            The main usage for apply_constraints is in automated tests using random values for initialization.
+     *
+     * @return Returns true if one (or more) constraint(s) were corrected, otherwise false.
+    */
+    bool apply_constraints()
+    {
+        bool corrected = false;
+        for (int i = 0; i < m_y.array().size(); i++) {
+            if (m_y.array()[i] < 0) {
+                log_warning("Constraint check: Compartment size {:d} changed from {:.4f} to {:d}", i, m_y.array()[i],
+                            0);
+                m_y.array()[i] = 0.;
+                corrected      = true;
+            }
+        }
+        return corrected;
+    }
+
+    /**
+     * @brief Checks whether all compartments have non-negative values and logs an error if constraint is not satisfied.
+     * @return Returns true if one or more constraints are not satisfied, false otherwise.
+     */
+    bool check_constraints() const
+    {
+        if ((m_y.array() < 0.).any()) {
+            log_error("Constraint check: At least one compartment size is smaller {}.", 0);
+            return true;
+        }
+        return false;
+    }
+
+private:
+    /**
+     * @brief Sets recursively the total number of (sub-)compartments over all groups.
+     * The number also corresponds to the size of the internal vector.
+     */
+    template <size_t Group = 0>
+    void set_count()
+    {
+        if constexpr (Group == 0) {
+            m_count = 0;
+        }
+        if constexpr (Group < num_groups) {
+            m_count += type_at_index_t<Group, LctStatesGroups>::Count;
+            set_count<Group + 1>();
+        }
+    }
+
+    size_t m_count; //< Number of groups stored.
+    InternalArrayType m_y{}; //< An array containing the number of people in the groups.
+};
+
+} // namespace mio
+
+#endif // MIO_EPI_LCT_POPULATIONS_H

cpp/models/lct_secir/CMakeLists.txt

@@ -11,4 +11,4 @@ target_include_directories(lct_secir PUBLIC
     $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
     $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
 )
-target_compile_options(lct_secir PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+target_compile_options(lct_secir PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})