1294 add LCT SECIR model with two disease strains

cpp/CMakeLists.txt

@@ -179,6 +179,7 @@ if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/ode_secirts)
     add_subdirectory(models/ode_secirvvs)
     add_subdirectory(models/lct_secir)
+    add_subdirectory(models/lct_secir_2_diseases)
     add_subdirectory(models/glct_secir)
     add_subdirectory(models/ide_secir)
     add_subdirectory(models/ide_seir)

cpp/examples/CMakeLists.txt

@@ -132,6 +132,10 @@ add_executable(lct_secir_example lct_secir.cpp)
 target_link_libraries(lct_secir_example PRIVATE memilio lct_secir)
 target_compile_options(lct_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(lct_secir_2_diseases_example lct_secir_2_diseases.cpp)
+target_link_libraries(lct_secir_2_diseases_example PRIVATE memilio lct_secir_2_diseases)
+target_compile_options(lct_secir_2_diseases_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
 add_executable(glct_secir_example glct_secir.cpp)
 target_link_libraries(glct_secir_example PRIVATE memilio glct_secir)
 target_compile_options(glct_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

cpp/examples/d_abm.cpp

@@ -1,5 +1,5 @@
 /* 
-* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
 *
 * Authors: Julia Bicker, René Schmieding
 *

cpp/examples/graph_abm.cpp

@@ -1,5 +1,5 @@
 /* 
-* Copyright (C) 2020-2024 MEmilio
+* Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Julia Bicker
 *

cpp/examples/lct_secir.cpp

@@ -37,31 +37,31 @@ 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;
+    constexpr size_t NumExposed = 1, NumInfectedNoSymptoms = 1, NumInfectedSymptoms = 1, NumInfectedSevere = 1,
+                     NumInfectedCritical = 1;
     using InfState                       = mio::lsecir::InfectionState;
     using LctState = mio::LctInfectionState<InfState, 1, NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms,
                                             NumInfectedSevere, NumInfectedCritical, 1, 1>;
-    using Model    = mio::lsecir::Model<LctState>;
+    using Model    = mio::lsecir::Model<LctState, LctState>;
     Model model;
 
     // Variable defines whether the class Initializer is used to define an initial vector from flows or whether a manually
     // defined initial vector is used to initialize the LCT model.
     bool use_initializer_flows = false;
 
-    ScalarType tmax = 10;
+    ScalarType tmax = 5;
 
     // Set Parameters.
-    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::TimeExposed>()[0]            = 3.;
+    model.parameters.get<mio::lsecir::TimeInfectedNoSymptoms>()[0] = 3.;
+    model.parameters.get<mio::lsecir::TimeInfectedSymptoms>()[0]   = 3.;
+    model.parameters.get<mio::lsecir::TimeInfectedSevere>()[0]     = 3.;
+    model.parameters.get<mio::lsecir::TimeInfectedCritical>()[0]   = 3.;
 
-    model.parameters.get<mio::lsecir::TransmissionProbabilityOnContact>()[0] = 0.05;
+    model.parameters.get<mio::lsecir::TransmissionProbabilityOnContact>()[0] = 0.1;
 
     mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::lsecir::ContactPatterns>();
-    contact_matrix[0]                       = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10));
+    contact_matrix[0]                       = mio::ContactMatrix(Eigen::MatrixXd::Constant(2, 2, 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.));
 
@@ -117,8 +117,7 @@ int main()
         // This method of defining the initial values using a vector of vectors is not necessary, but should remind you
         // how the entries of the initial value vector relate to the defined template parameters of the model or the number
         // of subcompartments. It is also possible to define the initial values directly.
-        std::vector<std::vector<ScalarType>> initial_populations = {{750}, {30, 20},          {20, 10, 10}, {50},
-                                                                    {50},  {10, 10, 5, 3, 2}, {20},         {10}};
+        std::vector<std::vector<ScalarType>> initial_populations = {{9000}, {1000}, {0}, {0}, {0}, {0}, {0}, {0}};
 
         // Assert that initial_populations has the right shape.
         if (initial_populations.size() != (size_t)InfState::Count) {
@@ -154,10 +153,10 @@ int main()
     }
 
     // Perform a simulation.
-    mio::TimeSeries<ScalarType> result = mio::simulate<ScalarType, Model>(0, tmax, 0.5, model);
+    mio::TimeSeries<ScalarType> result = mio::simulate<ScalarType, Model>(0, tmax, 0.1, model);
     // The simulation result is divided by subcompartments.
     // We call the function calculate_compartments 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 "}, 12, 4);
+    interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 12, 3);
 }

cpp/examples/lct_secir_2_diseases.cpp

@@ -0,0 +1,156 @@
+/* 
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Annika Jungklaus, 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.
+*/
+
+#include "lct_secir_2_diseases/model.h"
+#include "lct_secir_2_diseases/infection_state.h"
+#include "memilio/config.h"
+#include "memilio/utils/time_series.h"
+#include "memilio/epidemiology/lct_infection_state.h"
+#include "memilio/utils/logging.h"
+#include "memilio/compartments/simulation.h"
+#include "memilio/data/analyze_result.h"
+#include <vector>
+
+int main()
+{
+    // Simple example to demonstrate how to run a simulation using an LCT-SECIR-2-DISEASE 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_1a = 2, NumInfectedNoSymptoms_1a = 3, NumInfectedSymptoms_1a = 3,
+                     NumInfectedSevere_1a = 3, NumInfectedCritical_1a = 2, NumExposed_2a = 1,
+                     NumInfectedNoSymptoms_2a = 2, NumInfectedSymptoms_2a = 2, NumInfectedSevere_2a = 2,
+                     NumInfectedCritical_2a = 1, NumExposed_1b = 2, NumInfectedNoSymptoms_1b = 3,
+                     NumInfectedSymptoms_1b = 3, NumInfectedSevere_1b = 3, NumInfectedCritical_1b = 2,
+                     NumExposed_2b = 1, NumInfectedNoSymptoms_2b = 2, NumInfectedSymptoms_2b = 2,
+                     NumInfectedSevere_2b = 2, NumInfectedCritical_2b = 1;
+    using InfState = mio::lsecir2d::InfectionState;
+    using LctState = mio::LctInfectionState<
+        InfState, 1, NumExposed_1a, NumInfectedNoSymptoms_1a, NumInfectedSymptoms_1a, NumInfectedSevere_1a,
+        NumInfectedCritical_1a, 1, 1, NumExposed_2a, NumInfectedNoSymptoms_2a, NumInfectedSymptoms_2a,
+        NumInfectedSevere_2a, NumInfectedCritical_2a, NumExposed_1b, NumInfectedNoSymptoms_1b, NumInfectedSymptoms_1b,
+        NumInfectedSevere_1b, NumInfectedCritical_1b, 1, 1, NumExposed_2b, NumInfectedNoSymptoms_2b,
+        NumInfectedSymptoms_2b, NumInfectedSevere_2b, NumInfectedCritical_2b, 1>;
+    using Model = mio::lsecir2d::Model<LctState>;
+    Model model;
+
+    ScalarType tmax = 10;
+
+    // Set Parameters.
+    model.parameters.get<mio::lsecir2d::TimeExposed_a>()[0]            = 3.;
+    model.parameters.get<mio::lsecir2d::TimeExposed_b>()[0]            = 3.;
+    model.parameters.get<mio::lsecir2d::TimeInfectedNoSymptoms_a>()[0] = 3.;
+    model.parameters.get<mio::lsecir2d::TimeInfectedNoSymptoms_b>()[0] = 3.;
+    model.parameters.get<mio::lsecir2d::TimeInfectedSymptoms_a>()[0]   = 3.;
+    model.parameters.get<mio::lsecir2d::TimeInfectedSymptoms_b>()[0]   = 3.;
+    model.parameters.get<mio::lsecir2d::TimeInfectedSevere_a>()[0]     = 3.;
+    model.parameters.get<mio::lsecir2d::TimeInfectedSevere_b>()[0]     = 3.;
+    model.parameters.get<mio::lsecir2d::TimeInfectedCritical_a>()[0]   = 3.;
+    model.parameters.get<mio::lsecir2d::TimeInfectedCritical_b>()[0]   = 3.;
+
+    model.parameters.get<mio::lsecir2d::TransmissionProbabilityOnContact_a>()[0] = 0.1;
+    model.parameters.get<mio::lsecir2d::TransmissionProbabilityOnContact_b>()[0] = 0.1;
+
+    mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::lsecir2d::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::lsecir2d::RelativeTransmissionNoSymptoms_a>()[0] = 0.7;
+    model.parameters.get<mio::lsecir2d::RelativeTransmissionNoSymptoms_b>()[0] = 0.7;
+    model.parameters.get<mio::lsecir2d::RiskOfInfectionFromSymptomatic_a>()[0] = 0.25;
+    model.parameters.get<mio::lsecir2d::RiskOfInfectionFromSymptomatic_b>()[0] = 0.25;
+    model.parameters.get<mio::lsecir2d::RecoveredPerInfectedNoSymptoms_a>()[0] = 0.09;
+    model.parameters.get<mio::lsecir2d::RecoveredPerInfectedNoSymptoms_b>()[0] = 0.09;
+    model.parameters.get<mio::lsecir2d::SeverePerInfectedSymptoms_a>()[0]      = 0.2;
+    model.parameters.get<mio::lsecir2d::SeverePerInfectedSymptoms_b>()[0]      = 0.2;
+    model.parameters.get<mio::lsecir2d::CriticalPerSevere_a>()[0]              = 0.25;
+    model.parameters.get<mio::lsecir2d::CriticalPerSevere_b>()[0]              = 0.25;
+    model.parameters.get<mio::lsecir2d::DeathsPerCritical_a>()[0]              = 0.3;
+    model.parameters.get<mio::lsecir2d::DeathsPerCritical_b>()[0]              = 0.3;
+
+    // Simple example how to initialize model without flows.
+    // Define the initial values with the distribution of the population into subcompartments.
+    // This method of defining the initial values using a vector of vectors is not necessary, but should remind you
+    // how the entries of the initial value vector relate to the defined template parameters of the model or the number
+    // of subcompartments. It is also possible to define the initial values directly.
+    std::vector<std::vector<ScalarType>> initial_populations = {
+        {200},  {0, 0},  {30, 10, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0},     {0},       {0},       {0},
+        {0, 0}, {10, 0}, {0, 0},      {0},       {10, 0},   {30, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0},
+        {0},    {0},     {100},       {0, 0},    {0, 0},    {0, 0},     {0},       {0}};
+
+    // Assert that initial_populations has the right shape.
+    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)InfState::Susceptible].size() !=
+         LctState::get_num_subcompartments<InfState::Susceptible>()) ||
+        (initial_populations[(size_t)InfState::Exposed_1a].size() != NumExposed_1a) ||
+        (initial_populations[(size_t)InfState::InfectedNoSymptoms_1a].size() != NumInfectedNoSymptoms_1a) ||
+        (initial_populations[(size_t)InfState::InfectedSymptoms_1a].size() != NumInfectedSymptoms_1a) ||
+        (initial_populations[(size_t)InfState::InfectedSevere_1a].size() != NumInfectedSevere_1a) ||
+        (initial_populations[(size_t)InfState::InfectedCritical_1a].size() != NumInfectedCritical_1a) ||
+        (initial_populations[(size_t)InfState::Exposed_2a].size() != NumExposed_2a) ||
+        (initial_populations[(size_t)InfState::InfectedNoSymptoms_2a].size() != NumInfectedNoSymptoms_2a) ||
+        (initial_populations[(size_t)InfState::InfectedSymptoms_2a].size() != NumInfectedSymptoms_2a) ||
+        (initial_populations[(size_t)InfState::InfectedSevere_2a].size() != NumInfectedSevere_2a) ||
+        (initial_populations[(size_t)InfState::InfectedCritical_2a].size() != NumInfectedCritical_2a) ||
+        (initial_populations[(size_t)InfState::Recovered_a].size() !=
+         LctState::get_num_subcompartments<InfState::Recovered_a>()) ||
+        (initial_populations[(size_t)InfState::Dead_a].size() !=
+         LctState::get_num_subcompartments<InfState::Dead_a>()) ||
+        (initial_populations[(size_t)InfState::Exposed_1b].size() != NumExposed_1b) ||
+        (initial_populations[(size_t)InfState::InfectedNoSymptoms_1b].size() != NumInfectedNoSymptoms_1b) ||
+        (initial_populations[(size_t)InfState::InfectedSymptoms_1b].size() != NumInfectedSymptoms_1b) ||
+        (initial_populations[(size_t)InfState::InfectedSevere_1b].size() != NumInfectedSevere_1b) ||
+        (initial_populations[(size_t)InfState::InfectedCritical_1b].size() != NumInfectedCritical_1b) ||
+        (initial_populations[(size_t)InfState::Exposed_2b].size() != NumExposed_2b) ||
+        (initial_populations[(size_t)InfState::InfectedNoSymptoms_2b].size() != NumInfectedNoSymptoms_2b) ||
+        (initial_populations[(size_t)InfState::InfectedSymptoms_2b].size() != NumInfectedSymptoms_2b) ||
+        (initial_populations[(size_t)InfState::InfectedSevere_2b].size() != NumInfectedSevere_2b) ||
+        (initial_populations[(size_t)InfState::InfectedCritical_2b].size() != NumInfectedCritical_2b) ||
+        (initial_populations[(size_t)InfState::Recovered_ab].size() !=
+         LctState::get_num_subcompartments<InfState::Recovered_ab>())) {
+        mio::log_error("The length of at least one vector in initial_populations does not match the related number of "
+                       "subcompartments.");
+        return 1;
+    }
+    // Transfer the initial values in initial_populations to the model.
+    std::vector<ScalarType> flat_initial_populations;
+    for (auto&& vec : initial_populations) {
+        flat_initial_populations.insert(flat_initial_populations.end(), vec.begin(), vec.end());
+    }
+    for (size_t i = 0; i < LctState::Count; i++) {
+        model.populations[i] = flat_initial_populations[i];
+    }
+
+    // Perform a simulation.
+    mio::TimeSeries<ScalarType> result = mio::simulate<ScalarType, Model>(0, tmax, 0.5, model);
+    // The simulation result is divided by subcompartments.
+    // We call the function calculate_compartments 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",   "   E1a", "   C1a", "   I1a", "   H1a", "   U1a", "   Ra",
+                                      "   Da",  "   E2a", "   C2a", "   I2a", "   H2a", "   U2a", "   E1b",
+                                      "   C1b", "   I1b", "   H1b", "   U1b", "   Rb",  "   Db",  "   E2b",
+                                      "   C2b", "   I2b", "   H2b", "   U2b", "   Rab"},
+                                     6, 2);
+}

cpp/examples/smm.cpp

@@ -1,5 +1,5 @@
 /* 
-* Copyright (C) 2020-2024 German Aerospace Center (DLR-SC)
+* Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
 *
 * Authors: Julia Bicker, René Schmieding
 *

cpp/memilio/epidemiology/lct_infection_state.h

@@ -96,34 +96,17 @@ class LctInfectionState
 
         Eigen::VectorX<ScalarType> compartments((Eigen::Index)InfectionState::Count);
         // Use segment of the vector subcompartments of each InfectionState and sum up the values of subcompartments.
-        compartments[(Eigen::Index)InfectionState::Susceptible] = subcompartments[0];
-        compartments[(Eigen::Index)InfectionState::Exposed] =
-            subcompartments
-                .segment(get_first_index<InfectionState::Exposed>(), get_num_subcompartments<InfectionState::Exposed>())
-                .sum();
-        compartments[(Eigen::Index)InfectionState::InfectedNoSymptoms] =
-            subcompartments
-                .segment(get_first_index<InfectionState::InfectedNoSymptoms>(),
-                         get_num_subcompartments<InfectionState::InfectedNoSymptoms>())
-                .sum();
-        compartments[(Eigen::Index)InfectionState::InfectedSymptoms] =
-            subcompartments
-                .segment(get_first_index<InfectionState::InfectedSymptoms>(),
-                         get_num_subcompartments<InfectionState::InfectedSymptoms>())
-                .sum();
-        compartments[(Eigen::Index)InfectionState::InfectedSevere] =
-            subcompartments
-                .segment(get_first_index<InfectionState::InfectedSevere>(),
-                         get_num_subcompartments<InfectionState::InfectedSevere>())
-                .sum();
-        compartments[(Eigen::Index)InfectionState::InfectedCritical] =
-            subcompartments
-                .segment(get_first_index<InfectionState::InfectedCritical>(),
-                         get_num_subcompartments<InfectionState::InfectedCritical>())
-                .sum();
-        compartments[(Eigen::Index)InfectionState::Recovered] =
-            subcompartments[get_first_index<InfectionState::Recovered>()];
-        compartments[(Eigen::Index)InfectionState::Dead] = subcompartments[get_first_index<InfectionState::Dead>()];
+        for (int i = 0; i < (Eigen::Index)InfectionState::Count; i++) {
+            InfectionState State = static_cast<InfectionState>(i);
+            // first index of first subcompartment:
+            size_t index = 0;
+            for (size_t j = 0; j < (size_t)(State); j++) {
+                index = index + m_subcompartment_numbers[j];
+            }
+            // number of subcompartments:
+            size_t num_subcomp = m_subcompartment_numbers[(size_t)State];
+            compartments[i]    = subcompartments.segment(index, num_subcomp).sum();
+        }
 
         return compartments;
     }

cpp/models/lct_secir_2_diseases/CMakeLists.txt

@@ -0,0 +1,12 @@
+add_library(lct_secir_2_diseases
+    infection_state.h
+    model.h
+    model.cpp
+    parameters.h
+)
+target_link_libraries(lct_secir_2_diseases PUBLIC memilio)
+target_include_directories(lct_secir_2_diseases PUBLIC
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+target_compile_options(lct_secir_2_diseases PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})