diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 90cd4121b7..3d1cead066 100644 --- a/cpp/CMakeLists.txt +++ b/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) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index 980164d324..01f5c3b8e1 100644 --- a/cpp/examples/CMakeLists.txt +++ b/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}) diff --git a/cpp/examples/d_abm.cpp b/cpp/examples/d_abm.cpp index 0d5eea9863..28f2fb1e5d 100644 --- a/cpp/examples/d_abm.cpp +++ b/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 * diff --git a/cpp/examples/graph_abm.cpp b/cpp/examples/graph_abm.cpp index 34e89485c9..03b8ecdb57 100644 --- a/cpp/examples/graph_abm.cpp +++ b/cpp/examples/graph_abm.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2020-2024 MEmilio +* Copyright (C) 2020-2025 MEmilio * * Authors: Julia Bicker * diff --git a/cpp/examples/lct_secir_2_diseases.cpp b/cpp/examples/lct_secir_2_diseases.cpp new file mode 100644 index 0000000000..284582e71e --- /dev/null +++ b/cpp/examples/lct_secir_2_diseases.cpp @@ -0,0 +1,162 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Annika Jungklaus, Lena Ploetzke +* +* Contact: Martin J. Kuehn +* +* 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 + +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. + // The number of subcompartments can be chosen for most of the compartments: + 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; + // The compartment for Susceptible people and all compartments for Dead and Recovered people must have exactly one subcompartment: + constexpr size_t NumSusceptible = 1, NumDead_a = 1, NumDead_b = 1, NumRecovered_a = 1, NumRecovered_b = 1, + NumRecovered_ab = 1; + using InfState = mio::lsecir2d::InfectionState; + using LctState = + mio::LctInfectionState; + using Model = mio::lsecir2d::Model; + Model model; + + ScalarType tmax = 10; + + // Set Parameters. + model.parameters.get()[0] = 3.; + model.parameters.get()[0] = 3.; + model.parameters.get()[0] = 3.; + model.parameters.get()[0] = 3.; + model.parameters.get()[0] = 3.; + model.parameters.get()[0] = 3.; + model.parameters.get()[0] = 3.; + model.parameters.get()[0] = 3.; + model.parameters.get()[0] = 3.; + model.parameters.get()[0] = 3.; + + model.parameters.get()[0] = 0.1; + model.parameters.get()[0] = 0.1; + + mio::ContactMatrixGroup& contact_matrix = model.parameters.get(); + 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()[0] = 0.7; + model.parameters.get()[0] = 0.7; + model.parameters.get()[0] = 0.25; + model.parameters.get()[0] = 0.25; + model.parameters.get()[0] = 0.09; + model.parameters.get()[0] = 0.09; + model.parameters.get()[0] = 0.2; + model.parameters.get()[0] = 0.2; + model.parameters.get()[0] = 0.25; + model.parameters.get()[0] = 0.25; + model.parameters.get()[0] = 0.3; + model.parameters.get()[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> 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()) || + (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()) || + (initial_populations[(size_t)InfState::Dead_a].size() != + LctState::get_num_subcompartments()) || + (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())) { + 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 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 result = mio::simulate(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 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); +} diff --git a/cpp/examples/smm.cpp b/cpp/examples/smm.cpp index 1e7f5d285e..f3f216a6e5 100644 --- a/cpp/examples/smm.cpp +++ b/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 * diff --git a/cpp/memilio/epidemiology/lct_infection_state.h b/cpp/memilio/epidemiology/lct_infection_state.h index fbccebb4db..4badc7747c 100644 --- a/cpp/memilio/epidemiology/lct_infection_state.h +++ b/cpp/memilio/epidemiology/lct_infection_state.h @@ -96,34 +96,17 @@ class LctInfectionState Eigen::VectorX 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(), get_num_subcompartments()) - .sum(); - compartments[(Eigen::Index)InfectionState::InfectedNoSymptoms] = - subcompartments - .segment(get_first_index(), - get_num_subcompartments()) - .sum(); - compartments[(Eigen::Index)InfectionState::InfectedSymptoms] = - subcompartments - .segment(get_first_index(), - get_num_subcompartments()) - .sum(); - compartments[(Eigen::Index)InfectionState::InfectedSevere] = - subcompartments - .segment(get_first_index(), - get_num_subcompartments()) - .sum(); - compartments[(Eigen::Index)InfectionState::InfectedCritical] = - subcompartments - .segment(get_first_index(), - get_num_subcompartments()) - .sum(); - compartments[(Eigen::Index)InfectionState::Recovered] = - subcompartments[get_first_index()]; - compartments[(Eigen::Index)InfectionState::Dead] = subcompartments[get_first_index()]; + for (int i = 0; i < (Eigen::Index)InfectionState::Count; i++) { + InfectionState State = static_cast(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; } diff --git a/cpp/models/lct_secir_2_diseases/CMakeLists.txt b/cpp/models/lct_secir_2_diseases/CMakeLists.txt new file mode 100644 index 0000000000..ab0fcda2e2 --- /dev/null +++ b/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 + $ + $ +) +target_compile_options(lct_secir_2_diseases PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/models/lct_secir_2_diseases/README.md b/cpp/models/lct_secir_2_diseases/README.md new file mode 100644 index 0000000000..4e2e15c9e6 --- /dev/null +++ b/cpp/models/lct_secir_2_diseases/README.md @@ -0,0 +1,83 @@ +# LCT SECIR TWO DISEASES model + +This model describes infection with two independent (e.g. no co-infection) diseases a and b, based on the Linear Chain Trick. + +The Linear Chain Trick (LCT) provides the option to use Erlang-distributed stay times in the compartments through the use of subcompartments. +The normal ODE models have (possibly unrealistic) exponentially distributed stay times. +The LCT model can still be described by an ordinary differential equation system. + +For the concept of LCT models with one disease see +- Lena Plötzke, "Der Linear Chain Trick in der epidemiologischen Modellierung als Kompromiss zwischen gewöhnlichen und Integro-Differentialgleichungen", 2023. (https://elib.dlr.de/200381/, German only) +- P. J. Hurtado und A. S. Kirosingh, "Generalizations of the ‘Linear Chain Trick’: incorporating more flexible dwell time distributions into mean field ODE models“, 2019. (https://doi.org/10.1007/s00285-019-01412-w) + +For each infection there are the infection states Exposed, InfectedNoSymptoms, InfectedSymptoms, InfectedSevere, InfectedCritical, Recovered and Dead +The infectious compartments are InfectedNoSymptoms and InfectedSymptoms, all other compartments are considered to be not infectious + +There are two possibilities for a susceptible individual (since we assume no co-infection): +1. Get infected with disease a, then (if not dead) get infected with disease b +2. Get infected with disease b, then (if not dead) get infected with disease a + +This leads to the following compartments: +- `Susceptible` ($S$), may get infected with a or b (-> Exposed_1a or Exposed_1b) at any time +Infection 1a: +- `Exposed_1a` ($E_{1a}$), becomes InfectedNoSymptoms_1a after some time +- `InfectedNoSymptoms_1a` ($I_{NS, 1a}_$), becomes InfectedSymptoms_1a or Recovered_a after some time +- `InfectedSymptoms_1a` ($I_{Sy, 1a}$), becomes InfectedSevere_1a or Recovered_a after some time +- `InfectedSevere_1a` ($I_{Sev, 1a}$), becomes InfectedCritical_1a or Recovered_a after some time +- `InfectedCritical_1a` ($I_{Cr, 1a}$), becomes Recovered_a or Dead_a after some time +- `Recovered_a` ($R_a$), immune to a, may get infected with b (-> Exposed_2b) at any time +- `Dead_a` ($D_a$), absorbing state +Infection 1b: +- `Exposed_1b` ($E_{1b}$), becomes InfectedNoSymptoms_1b after some time +- `InfectedNoSymptoms_1b` ($I_{NS, 1b}$), becomes InfectedSymptoms_1b or Recovered_1b after some time +- `InfectedSymptoms_1b` ($I_{Sy, 1b}$), becomes InfectedSevere_1b or Recovered_b after some time +- `InfectedSevere_1b` ($I_{Sev, 1b}$), becomes InfectedCritical_1b or Recovered_b after some time +- `InfectedCritical_1b` ($I_{Cr, 1b}$), becomes Recovered_b or Dead_b after some time +- `Recovered_b` ($R_b$), immune to b, may get infected with b (-> Exposed_2a) at any time +- `Dead_b` ($D_b$), absorbing state +Infection 2a: +- `Exposed_2a` ($E_{2a}$), becomes InfectedNoSympotoms_2a after some time +- `InfectedNoSymptoms_2a` ($I_{NS, 2a}$), becomes InfectedSymptoms_2a or Recovered_ab after some time +- `InfectedSymptoms_2a` ($I_{Sy, 2a}$), becomes InfectedSevere_2a or Recovered_ab after some time +- `InfectedSevere_2a` ($I_{Sev, 2a}$), becomes InfectedCritical_2a or Recovered_ab after some time +- `InfectedCritical_2a` ($I_{Cr, 2a}$), becomes Recovered_ab or Dead_a after some time +- `Recovered_ab` ($R_{ab}$), absorbing state, immune to a and b +Infection 2b: +- `Exposed_2b` ($E_{2b}$), becomes infected after some time +- `InfectedNoSymptoms_2b` ($I_{NS,2b}$), becomes InfectedSymptoms_2b or Recovered_ab after some time +- `InfectedSymptoms_2b` ($I_{Sy,2b}$), becomes InfectedSevere_2b or Recovered_ab after some time +- `InfectedSevere_2b` ($I_{Sev, 2b}$), becomes InfectedCritical_2b or Recovered_ab after some time +- `InfectedCritical_2b` ($I_{Cr, 2b}$), becomes Recovered_ab or Dead_b after some time + + +It is possible to include subcompartments for the five compartments Exposed, InfectedNoSymptoms, InfectedSymptoms, InfectedSevere and InfectedCritical. +You can divide the population according to different groups, e.g. AgeGroups or gender and choose parameters according to groups. + +The parameters depend on the disease (a or b), the number of subcompartments depends on the individual compartment (since it can be set independently). + + +| Mathematical variable | C++ variable name | Description | +|---------------------------- | -------------------------------------- | ------------------------------------------------------------------------------------------- | +| $\phi$ | `ContactPatterns` | Average number of contacts of a person per day. | +| $\rho_i$ | `TransmissionProbabilityOnContact_i` | Transmission risk for people located in the susceptible compartments, i in {a,b}. | +| $\xi_{I_{NS},i}$ | `RelativeTransmissionNoSymptoms_i` | Proportion of nonsymptomatically infected people who are not isolated, i in {a,b}. | +| $\xi_{I_{Sy},i}$ | `RiskOfInfectionFromSymptomatic_i` | Proportion of infected people with symptoms who are not isolated, i in {a,b}. | +| $n_{E,i}$ | Defined in `LctStates` | Number of subcompartments of the Exposed compartment, i in {1a, 2a, 1b,2b}. | +| $n_{NS,i}$ | Defined in `LctStates` | Number of subcompartments of the InfectedNoSymptoms compartment, i in {1a, 2a, 1b,2b}. | +| $n_{Sy,i}$ | Defined in `LctStates` | Number of subcompartments of the InfectedSymptoms compartment, i in {1a, 2a, 1b,2b}. | +| $n_{Sev,i}$ | Defined in `LctStates` | Number of subcompartments of the InfectedSevere compartment, i in {1a, 2a, 1b,2b}. | +| $n_{Cr,i}$ | Defined in `LctStates` | Number of subcompartments of the InfectedCritical compartment, i in {1a, 2a, 1b,2b}. | +| $T_{E,i}$ | `TimeExposed` | Average time in days an individual stays in the Exposed compartment, i in {a,b}. | +| $T_{I_{NS},i}$ | `TimeInfectedNoSymptoms` | Average time in days an individual stays in the InfectedNoSymptoms compartment, i in {a,b}. | +| $T_{I_{Sy},i}$ | `TimeInfectedSymptoms` | Average time in days an individual stays in the InfectedSymptoms compartmen, i in {a,b}t. | +| $T_{I_{Sev},i}$ | `TimeInfectedSevere` | Average time in days an individual stays in the InfectedSevere compartment, i in {a,b}. | +| $T_{I_{Cr},i}$ | `TimeInfectedCritical` | Average time in days an individual stays in the InfectedCritical compartment, i in {a,b}. | +| $\mu_{I_{NS},i}^{R}$ | `RecoveredPerInfectedNoSymptoms` | Probability of transition from compartment InfectedNoSymptoms to Recovered, i in {a,b}. | +| $\mu_{I_{Sy},i}^{I_{Sev}}$ | `SeverePerInfectedSymptoms` | Probability of transition from compartment InfectedSymptoms to InfectedSevere, i in {a,b}. | +| $\mu_{I_{Sev},i}^{I_{Cr}}$ | `CriticalPerSevere` | Probability of transition from compartment InfectedSevere to InfectedCritical, i in {a,b}. | +| $\mu_{I_{Cr},i}^{D}$ | `DeathsPerCritical` | Probability of dying when in compartment InfectedCritical, i in {a,b}. | + + +## Examples + +A simple example can be found at [LCT2D minimal example](../../examples/lct_secir_2_diseases.cpp). diff --git a/cpp/models/lct_secir_2_diseases/infection_state.h b/cpp/models/lct_secir_2_diseases/infection_state.h new file mode 100644 index 0000000000..92ef39c442 --- /dev/null +++ b/cpp/models/lct_secir_2_diseases/infection_state.h @@ -0,0 +1,126 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Annika Jungklaus, Lena Ploetzke +* +* Contact: Martin J. Kuehn +* +* 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 LCT_SECIR_2_DISEASES_INFECTIONSTATE_H +#define LCT_SECIR_2_DISEASES_INFECTIONSTATE_H + +namespace mio +{ +namespace lsecir2d +{ + +/** + * @brief The InfectionState enum describes the basic + * categories for the infection state of persons. + */ +enum class InfectionState +{ + Susceptible = 0, + // Notation: State_[Infection number][disease] + // first infection with disease a + Exposed_1a = 1, + InfectedNoSymptoms_1a = 2, + InfectedSymptoms_1a = 3, + InfectedSevere_1a = 4, + InfectedCritical_1a = 5, + Recovered_a = 6, + Dead_a = 7, + // second infection with disease a + Exposed_2a = 8, + InfectedNoSymptoms_2a = 9, + InfectedSymptoms_2a = 10, + InfectedSevere_2a = 11, + InfectedCritical_2a = 12, + // first infection with disease b + Exposed_1b = 13, + InfectedNoSymptoms_1b = 14, + InfectedSymptoms_1b = 15, + InfectedSevere_1b = 16, + InfectedCritical_1b = 17, + Recovered_b = 18, + Dead_b = 19, + // second infection with disease b + Exposed_2b = 20, + InfectedNoSymptoms_2b = 21, + InfectedSymptoms_2b = 22, + InfectedSevere_2b = 23, + InfectedCritical_2b = 24, + // Recovered from both diseases + Recovered_ab = 25, + Count = 26 +}; + +/** + * @brief The InfectionTransition enum describes the possible + * transitions of the infectious state of persons. + */ +enum class InfectionTransition +{ + // first infection with a + SusceptibleToExposed_1a = 0, + Exposed_1aToInfectedNoSymptoms_1a = 1, + InfectedNoSymptoms_1aToInfectedSymptoms_1a = 2, + InfectedNoSymptoms_1aToRecovered_a = 3, + InfectedSymptoms_1aToInfectedSevere_1a = 4, + InfectedSymptoms_1aToRecovered_a = 5, + InfectedSevere_1aToInfectedCritical_1a = 6, + InfectedSevere_1aToRecovered_a = 7, + InfectedCritical_1aToDead_a = 8, + InfectedCritical_1aToRecovered_a = 9, + // second infection with a + Recovered_bToExposed_2a = 10, + Exposed_2aToInfectedNoSymptoms_2a = 11, + InfectedNoSymptoms_2aToInfectedSymptoms_2a = 12, + InfectedNoSymptoms_2aToRecovered_ab = 13, + InfectedSymptoms_2aToInfectedSevere_2a = 14, + InfectedSymptoms_2aToRecovered_ab = 15, + InfectedSevere_2aToInfectedCritical_2a = 16, + InfectedSevere_2aToRecovered_ab = 17, + InfectedCritical_2aToDead_a = 18, + InfectedCritical_2aToRecovered_ab = 19, + // first infection with b + SusceptibleToExposed_1b = 20, + Exposed_1bToInfectedNoSymptoms_1b = 21, + InfectedNoSymptoms_1bToInfectedSymptoms_1b = 22, + InfectedNoSymptoms_1bToRecovered_b = 23, + InfectedSymptoms_1bToInfectedSevere_1b = 24, + InfectedSymptoms_1bToRecovered_b = 25, + InfectedSevere_1bToInfectedCritical_1b = 26, + InfectedSevere_1bToRecovered_b = 27, + InfectedCritical_1bToDead_b = 28, + InfectedCritical_1bToRecovered_b = 29, + // second infection with b + Recovered_aToExposed_2b = 30, + Exposed_2bToInfectedNoSymptoms_2b = 31, + InfectedNoSymptoms_2bToInfectedSymptoms_2b = 32, + InfectedNoSymptoms_2bToRecovered_ab = 33, + InfectedSymptoms_2bToInfectedSevere_2b = 34, + InfectedSymptoms_2bToRecovered_ab = 35, + InfectedSevere_2bToInfectedCritical_2b = 36, + InfectedSevere_2bToRecovered_ab = 37, + InfectedCritical_2bToDead_b = 38, + InfectedCritical_2bToRecovered_ab = 39, + Count = 40 +}; + +} // namespace lsecir2d +} // namespace mio + +#endif // LCT_SECIR_2_DISEASES_INFECTIONSTATE_H \ No newline at end of file diff --git a/cpp/models/lct_secir_2_diseases/model.cpp b/cpp/models/lct_secir_2_diseases/model.cpp new file mode 100644 index 0000000000..77a3b42578 --- /dev/null +++ b/cpp/models/lct_secir_2_diseases/model.cpp @@ -0,0 +1,29 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Annika Jungklaus, Lena Ploetzke +* +* Contact: Martin J. Kuehn +* +* 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" + +namespace mio +{ +namespace lsecir2d +{ + +} // namespace lsecir2d +} // namespace mio \ No newline at end of file diff --git a/cpp/models/lct_secir_2_diseases/model.h b/cpp/models/lct_secir_2_diseases/model.h new file mode 100644 index 0000000000..4ebed88b47 --- /dev/null +++ b/cpp/models/lct_secir_2_diseases/model.h @@ -0,0 +1,741 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Annika Jungkalus, Lena Ploetzke +* +* Contact: Martin J. Kuehn +* +* 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 LCT_SECIR_2_DISEASE_MODEL_H +#define LCT_SECIR_2_DISEASE_MODEL_H + +#include "lct_secir_2_diseases/parameters.h" +#include "lct_secir_2_diseases/infection_state.h" +#include "memilio/compartments/compartmental_model.h" +#include "memilio/epidemiology/lct_populations.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "memilio/utils/logging.h" +#include "memilio/utils/type_list.h" +#include "memilio/utils/metaprogramming.h" + +namespace mio +{ +namespace lsecir2d +{ +/** + * @brief Class that defines an LCT-SECIR-2-DISEASE model. + * + * @tparam LctStates The LCT2D model can work with any number of LctStates, where each LctState corresponds to a group, + * e.g. one AgeGroup. The purpose of the LctStates is to define the number of subcompartments for each InfectionState. + * If you do not want to divide the population into groups, just use one LctState. + * If you want to divide the population according to more than one category, e.g. sex and age, + * you have to specify one LctState for each pair of groups, e.g. for (female, A00-A04), (female, A05-A14) etc. + * This is because the number of subcompartments can be different for each group. + * Therefore, the number of LctStates also determines the number of groups. + */ +template +class Model + : public CompartmentalModel, Parameters> +{ +public: + using LctStatesGroups = TypeList; + using Base = CompartmentalModel, Parameters>; + using typename Base::ParameterSet; + using typename Base::Populations; + static size_t constexpr num_groups = sizeof...(LctStates); + static_assert(num_groups >= 1, "The number of LctStates provided should be at least one."); + + /// @brief Default constructor. + Model() + : Base(Populations(), ParameterSet(num_groups)) + { + } + + /** @brief Constructor using Populations and ParameterSet. + * @param[in] pop An instance of the Populations class. + * @param[in] params Parameters used to be used in the simulation. + */ + Model(const Populations& pop, const ParameterSet& params) + : Base(pop, params) + { + } + + /** + * @brief Evaluates the right-hand-side f of the ODE dydt = f(y, t). + * + * The LCT-SECIR-2-DISEASE model is defined through ordinary differential equations of the form dydt = f(y, t). + * y is a vector containing the number of individuals for each (sub-) compartment. + * This function evaluates the right-hand-side f of the ODE and can be used in an ODE solver. + * @param[in] pop The current state of the population in the geographic unit we are considering. + * @param[in] y The current state of the model (or a subpopulation) as a flat array. + * @param[in] t The current time. + * @param[out] dydt A reference to the calculated output. + */ + void get_derivatives(Eigen::Ref> pop, + Eigen::Ref> y, ScalarType t, + Eigen::Ref> dydt) const override + { + // Vectors are sorted such that we first have all InfectionState%s for AgeGroup 0, + // afterwards all for AgeGroup 1 and so on. + dydt.setZero(); + get_derivatives_impl(pop, y, t, dydt); + } + + /** + * @brief Cumulates a simulation result with subcompartments to produce a result that divides the population only + * into the infection states defined in InfectionState. + * + * If the model is used for simulation, we will get a result in form of a TimeSeries with infection states divided + * in subcompartments. + * The function calculates a TimeSeries without subcompartments from another TimeSeries with subcompartments. + * This is done by summing up the numbers in the subcompartments. + * @param[in] subcompartments_ts Result of a simulation with the model. + * @return Result of the simulation divided in infection states without subcompartments. + * Returns TimeSeries with values -1 if calculation is not possible. + */ + TimeSeries calculate_compartments(const TimeSeries& subcompartments_ts) const + { + Eigen::Index count_InfStates = (Eigen::Index)InfectionState::Count; + Eigen::Index num_compartments = count_InfStates * num_groups; + TimeSeries compartments_ts(num_compartments); + if (!(this->populations.get_num_compartments() == (size_t)subcompartments_ts.get_num_elements())) { + log_error("Result does not match InfectionState of the Model."); + Eigen::VectorX wrong_size = Eigen::VectorX::Constant(num_compartments, -1); + compartments_ts.add_time_point(-1, wrong_size); + return compartments_ts; + } + Eigen::VectorX compartments(num_compartments); + for (Eigen::Index timepoint = 0; timepoint < subcompartments_ts.get_num_time_points(); ++timepoint) { + compress_vector(subcompartments_ts[timepoint], compartments); + compartments_ts.add_time_point(subcompartments_ts.get_time(timepoint), compartments); + } + + return compartments_ts; + } + + /** + * @brief Checks that the model satisfies all constraints (e.g. parameter or population constraints). + * @return Returns true if one or more constraints are not satisfied, false otherwise. + */ + bool check_constraints() const + { + + return (check_constraints_impl() || this->parameters.check_constraints() || + this->populations.check_constraints()); + } + +private: + /** + * @brief Converts a vector with subcompartments in a vector without subcompartments + * by summing up subcompartment values. + * This is done recursively for each group which corresponds to a slice of the vector. + * + * @tparam group The group specifying the slice of the vector being considered. + * @param[in] subcompartments The vector that should be converted. + * @param[out] compartments Reference to the vector where the output is stored. + */ + template + void compress_vector(const Eigen::VectorX& subcompartments, + Eigen::VectorX& compartments) const + { + static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid."); + using LctStateGroup = type_at_index_t; + + // Define first index of the group Group in a vector including all compartments without a resolution + // in subcompartments. + Eigen::Index count_InfStates = (Eigen::Index)InfectionState::Count; + Eigen::Index first_index_group_comps = Group * count_InfStates; + + // Use function from the LctState of the Group to calculate the vector without subcompartments + // using the corresponding vector with subcompartments. + compartments.segment(first_index_group_comps, count_InfStates) = + LctStateGroup::calculate_compartments(subcompartments.segment( + this->populations.template get_first_index_of_group(), LctStateGroup::Count)); + + // Function call for next group if applicable. + if constexpr (Group + 1 < num_groups) { + compress_vector(subcompartments, compartments); + } + } + + /** + * @brief Evaluates the right-hand-side f of the ODE dydt = f(y, t) recursively for each group. + * + * See also the function get_derivative. + * For each group, one slice of the output vector is calculated. + * @tparam group The group specifying the slice of the vector being considered. + * @param[in] pop The current state of the population in the geographic unit we are considering. + * @param[in] y The current state of the model (or a subpopulation) as a flat array. + * @param[in] t The current time. + * @param[out] dydt A reference to the calculated output. + */ + template + void get_derivatives_impl(Eigen::Ref> pop, + Eigen::Ref> y, ScalarType t, + Eigen::Ref> dydt) const + { + static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid."); + using LctStateGroup = type_at_index_t; + + size_t first_index_group = this->populations.template get_first_index_of_group(); + auto params = this->parameters; + ScalarType flow = 0; + + // Indices of first subcompartment of the InfectionState for the group in the vectors. + size_t Ei_1a_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t Ei_2a_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t Ei_1b_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t Ei_2b_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t INSi_1a_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t INSi_2a_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t INSi_1b_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t INSi_2b_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ISyi_1a_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ISyi_2a_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ISyi_1b_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ISyi_2b_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ISevi_1a_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ISevi_2a_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ISevi_1b_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ISevi_2b_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ICri_1a_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ICri_2a_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ICri_1b_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t ICri_2b_first_index = + first_index_group + LctStateGroup::template get_first_index(); + size_t Ri_a = first_index_group + LctStateGroup::template get_first_index(); + size_t Ri_b = first_index_group + LctStateGroup::template get_first_index(); + size_t Ri_ab = first_index_group + LctStateGroup::template get_first_index(); + size_t Di_a = first_index_group + LctStateGroup::template get_first_index(); + size_t Di_b = first_index_group + LctStateGroup::template get_first_index(); + + // Calculate derivative of the Susceptible compartment. + // outflow generated by disease a and disease b both + double part_a = 0.; // part of people getting infected with disease a + double part_b = 0.; // part of people getting infected with disease b + interact(pop, y, t, dydt, &part_a, &part_b, first_index_group, 2); + // split flow + double div_part_both = + ((part_a + part_b) < Limits::zero_tolerance()) ? 0.0 : 1.0 / (part_a + part_b); + + // Start with derivatives of First Infections (X_1a, X_1b) + // Calculate derivative of the Exposed_1x compartments, split flow from Susceptible. + // Exposed 1 a: + dydt[Ei_1a_first_index] = -dydt[first_index_group] * part_a * div_part_both; + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); subcomp++) { + // Variable flow stores the value of the flow from one subcompartment to the next one. + // Ei_1a_first_index + subcomp is always the index of a (sub-)compartment of Exposed and Ei_1a_first_index + // + subcomp + 1 can also be the index of the first (sub-)compartment of InfectedNoSymptoms. + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[Ei_1a_first_index + subcomp]; + // Subtract flow from dydt[Ei_1a_first_index + subcomp] and add to next subcompartment. + dydt[Ei_1a_first_index + subcomp] -= flow; + dydt[Ei_1a_first_index + subcomp + 1] = flow; + } + // Exposed 1 b: + dydt[Ei_1b_first_index] = -dydt[first_index_group] * part_b * div_part_both; + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); subcomp++) { + // Variable flow stores the value of the flow from one subcompartment to the next one. + // Ei_1a_first_index + subcomp is always the index of a (sub-)compartment of Exposed and Ei_1b_first_index + // + subcomp + 1 can also be the index of the first (sub-)compartment of InfectedNoSymptoms. + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[Ei_1b_first_index + subcomp]; + // Subtract flow from dydt[Ei_1b_first_index + subcomp] and add to next subcompartment. + dydt[Ei_1b_first_index + subcomp] -= flow; + dydt[Ei_1b_first_index + subcomp + 1] = flow; + } + + // Calculate derivative of the InfectedNoSymptoms_1a compartment + // flow from each subcompartment to the next + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = + (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[INSi_1a_first_index + subcomp]; + dydt[INSi_1a_first_index + subcomp] -= flow; + dydt[INSi_1a_first_index + subcomp + 1] = flow; + } + // Calculate derivative of the InfectedNoSymptoms_1b compartment + // flow from each subcompartment to the next + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = + (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[INSi_1b_first_index + subcomp]; + dydt[INSi_1b_first_index + subcomp] -= flow; + dydt[INSi_1b_first_index + subcomp + 1] = flow; + } + + // Calculate derivative of the InfectedSymptoms_1a compartment. + // Flow from last (sub-) compartment of C_1a must be split between + // the first subcompartment of InfectedSymptoms_1a and Recovered_a. + dydt[Ri_a] = dydt[ISyi_1a_first_index] * params.template get()[Group]; + dydt[ISyi_1a_first_index] = + dydt[ISyi_1a_first_index] * (1 - params.template get()[Group]); + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ISyi_1a_first_index + subcomp]; + dydt[ISyi_1a_first_index + subcomp] -= flow; + dydt[ISyi_1a_first_index + subcomp + 1] = flow; + } + // Calculate derivative of the InfectedSymptoms_1b compartment. + // Flow from last (sub-) compartment of C_1b must be split between + // the first subcompartment of InfectedSymptoms_1b and Recovered_b. + dydt[Ri_b] = dydt[ISyi_1b_first_index] * params.template get()[Group]; + dydt[ISyi_1b_first_index] = + dydt[ISyi_1b_first_index] * (1 - params.template get()[Group]); + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ISyi_1b_first_index + subcomp]; + dydt[ISyi_1b_first_index + subcomp] -= flow; + dydt[ISyi_1b_first_index + subcomp + 1] = flow; + } + + // Calculate derivative of the InfectedSevere_1a compartment. + // again split the flow from the last subcompartment of I_1a + dydt[Ri_a] += dydt[ISevi_1a_first_index] * (1 - params.template get()[Group]); + dydt[ISevi_1a_first_index] = + dydt[ISevi_1a_first_index] * params.template get()[Group]; + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ISevi_1a_first_index + subcomp]; + dydt[ISevi_1a_first_index + subcomp] -= flow; + dydt[ISevi_1a_first_index + subcomp + 1] = flow; + } + // Calculate derivative of the InfectedSevere_1b compartment. + // again split the flow from the last subcompartment of I_1b + dydt[Ri_b] += dydt[ISevi_1b_first_index] * (1 - params.template get()[Group]); + dydt[ISevi_1b_first_index] = + dydt[ISevi_1b_first_index] * params.template get()[Group]; + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ISevi_1b_first_index + subcomp]; + dydt[ISevi_1b_first_index + subcomp] -= flow; + dydt[ISevi_1b_first_index + subcomp + 1] = flow; + } + + // Calculate derivative of the InfectedCritical compartment. + // again split flow from last subcompartment of H_1a between R_a and U_1a + dydt[Ri_a] += dydt[ICri_1a_first_index] * (1 - params.template get()[Group]); + dydt[ICri_1a_first_index] = dydt[ICri_1a_first_index] * params.template get()[Group]; + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments() - 1; + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ICri_1a_first_index + subcomp]; + dydt[ICri_1a_first_index + subcomp] -= flow; + dydt[ICri_1a_first_index + subcomp + 1] = flow; + } + // Calculate derivative of the InfectedCritical compartment. + // again split flow from last subcompartment of H_1b between R_b and U_1b + dydt[Ri_b] += dydt[ICri_1b_first_index] * (1 - params.template get()[Group]); + dydt[ICri_1b_first_index] = dydt[ICri_1b_first_index] * params.template get()[Group]; + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments() - 1; + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ICri_1b_first_index + subcomp]; + dydt[ICri_1b_first_index + subcomp] -= flow; + dydt[ICri_1b_first_index + subcomp + 1] = flow; + } + + // Flow from InfectedCritical has to be divided between Recovered and Dead. + // Must be calculated separately in order not to overwrite the already calculated values ​​for Recovered. + // for InefectedCritical_1a: + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * + y[ICri_1a_first_index + + LctStateGroup::template get_num_subcompartments() - 1]; + dydt[ICri_1a_first_index + + LctStateGroup::template get_num_subcompartments() - 1] -= flow; + dydt[Ri_a] = dydt[Ri_a] + (1 - params.template get()[Group]) * flow; + dydt[Di_a] = dydt[Di_a] + params.template get()[Group] * flow; + // Flow from InfectedCritical_1b has to be divided between Recovered_b and Dead_b. + // Must be calculated separately in order not to overwrite the already calculated values ​​for Recovered. + // Outflow from InfectedCritical_1b: + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * + y[ICri_1b_first_index + + LctStateGroup::template get_num_subcompartments() - 1]; + dydt[ICri_1b_first_index + + LctStateGroup::template get_num_subcompartments() - 1] -= flow; + dydt[Ri_b] = dydt[Ri_b] + (1 - params.template get()[Group]) * flow; + dydt[Di_b] = dydt[Di_b] + params.template get()[Group] * flow; + + // Second Infection (X_2a, X_2b) + // outflow from Recovered, people getting infected for the second time + double temp_Ra = dydt[Ri_a]; + interact<0, 0>(pop, y, t, dydt, &part_a, &part_b, Ri_a, + 1); // outflow from R_a is only affected by b, no age groups + double temp_Rb = dydt[Ri_b]; + interact<0, 0>(pop, y, t, dydt, &part_a, &part_b, Ri_b, + 0); // outflow from R_b is only affected by a, no age groups + + // Calculate derivative of the Exposed_2i compartments + // Exposed 2 a + dydt[Ei_2a_first_index] = -(dydt[Ri_b] - temp_Rb); + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); subcomp++) { + // Variable flow stores the value of the flow from one subcompartment to the next one. + // Ei_2a_first_index + subcomp is always the index of a (sub-)compartment of Exposed and Ei_2a_first_index + // + subcomp + 1 can also be the index of the first (sub-)compartment of InfectedNoSymptoms. + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[Ei_2a_first_index + subcomp]; + // Subtract flow from dydt[Ei_2a_first_index + subcomp] and add to next subcompartment. + dydt[Ei_2a_first_index + subcomp] -= flow; + dydt[Ei_2a_first_index + subcomp + 1] = flow; + } + // Exposed 2 b + dydt[Ei_2b_first_index] = -(dydt[Ri_a] - temp_Ra); + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); subcomp++) { + // Variable flow stores the value of the flow from one subcompartment to the next one. + // Ei_2b_first_index + subcomp is always the index of a (sub-)compartment of Exposed and Ei_2b_first_index + // + subcomp + 1 can also be the index of the first (sub-)compartment of InfectedNoSymptoms. + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[Ei_2b_first_index + subcomp]; + // Subtract flow from dydt[Ei_2b_first_index + subcomp] and add to next subcompartment. + dydt[Ei_2b_first_index + subcomp] -= flow; + dydt[Ei_2b_first_index + subcomp + 1] = flow; + } + + // Calculate derivative of the InfectedNoSymptoms_2a compartment. + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = + (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[INSi_2a_first_index + subcomp]; + dydt[INSi_2a_first_index + subcomp] -= flow; + dydt[INSi_2a_first_index + subcomp + 1] = flow; + } + // Calculate derivative of the InfectedNoSymptoms_2b compartment. + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = + (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[INSi_2b_first_index + subcomp]; + dydt[INSi_2b_first_index + subcomp] -= flow; + dydt[INSi_2b_first_index + subcomp + 1] = flow; + } + + // Calculate derivative of the InfectedSymptoms_2a compartment. + // Flow from last (sub-) compartment of C_2a must be split between + // the first subcompartment of InfectedSymptoms_2a and Recovered_ab. + dydt[Ri_ab] += dydt[ISyi_2a_first_index] * params.template get()[Group]; + dydt[ISyi_2a_first_index] = + dydt[ISyi_2a_first_index] * (1 - params.template get()[Group]); + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ISyi_2a_first_index + subcomp]; + dydt[ISyi_2a_first_index + subcomp] -= flow; + dydt[ISyi_2a_first_index + subcomp + 1] = flow; + } + // Calculate derivative of the InfectedSymptoms_2b compartment. + // Flow from last (sub-) compartment of C_2b must be split between + // the first subcompartment of InfectedSymptoms_2b and Recovered_ab. + dydt[Ri_ab] += dydt[ISyi_2b_first_index] * params.template get()[Group]; + dydt[ISyi_2b_first_index] = + dydt[ISyi_2b_first_index] * (1 - params.template get()[Group]); + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ISyi_2b_first_index + subcomp]; + dydt[ISyi_2b_first_index + subcomp] -= flow; + dydt[ISyi_2b_first_index + subcomp + 1] = flow; + } + + // Calculate derivative of the InfectedSevere_2a compartment. + // again split the flow from the last subcompartment of I_2a + dydt[Ri_ab] += dydt[ISevi_2a_first_index] * (1 - params.template get()[Group]); + dydt[ISevi_2a_first_index] = + dydt[ISevi_2a_first_index] * params.template get()[Group]; + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ISevi_2a_first_index + subcomp]; + dydt[ISevi_2a_first_index + subcomp] -= flow; + dydt[ISevi_2a_first_index + subcomp + 1] = flow; + } + // Calculate derivative of the InfectedSevere compartment. + // again split the flow from the last subcompartment of I_2b + dydt[Ri_ab] += dydt[ISevi_2b_first_index] * (1 - params.template get()[Group]); + dydt[ISevi_2b_first_index] = + dydt[ISevi_2b_first_index] * params.template get()[Group]; + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments(); + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ISevi_2b_first_index + subcomp]; + dydt[ISevi_2b_first_index + subcomp] -= flow; + dydt[ISevi_2b_first_index + subcomp + 1] = flow; + } + + // Calculate derivative of the InfectedCritical compartment. + // again split flow from last subcompartment of H_2a between R_ab and U_2a + dydt[Ri_ab] += dydt[ICri_2a_first_index] * (1 - params.template get()[Group]); + dydt[ICri_2a_first_index] = dydt[ICri_2a_first_index] * params.template get()[Group]; + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments() - 1; + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ICri_2a_first_index + subcomp]; + dydt[ICri_2a_first_index + subcomp] -= flow; + dydt[ICri_2a_first_index + subcomp + 1] = flow; + } + // Calculate derivative of the InfectedCritical compartment. + // again split flow from last subcompartment of H_1b between R_b and U_1b + dydt[Ri_ab] += dydt[ICri_2b_first_index] * (1 - params.template get()[Group]); + dydt[ICri_2b_first_index] = dydt[ICri_2b_first_index] * params.template get()[Group]; + for (size_t subcomp = 0; + subcomp < LctStateGroup::template get_num_subcompartments() - 1; + subcomp++) { + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * y[ICri_2b_first_index + subcomp]; + dydt[ICri_2b_first_index + subcomp] -= flow; + dydt[ICri_2b_first_index + subcomp + 1] = flow; + } + + // Last flow from InfectedCritical has to be divided between Recovered and Dead. + // Must be calculated separately in order not to overwrite the already calculated values ​​for Recovered. + // for 2a + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * + y[ICri_2a_first_index + + LctStateGroup::template get_num_subcompartments() - 1]; + dydt[ICri_2a_first_index + + LctStateGroup::template get_num_subcompartments() - 1] -= flow; + dydt[Ri_ab] = dydt[Ri_ab] + (1 - params.template get()[Group]) * flow; + dydt[Di_a] = dydt[Di_a] + params.template get()[Group] * flow; + // for 2b + flow = (ScalarType)LctStateGroup::template get_num_subcompartments() * + (1 / params.template get()[Group]) * + y[ICri_2b_first_index + + LctStateGroup::template get_num_subcompartments() - 1]; + dydt[ICri_2b_first_index + + LctStateGroup::template get_num_subcompartments() - 1] -= flow; + dydt[Ri_ab] = dydt[Ri_ab] + (1 - params.template get()[Group]) * flow; + dydt[Di_b] = dydt[Di_b] + params.template get()[Group] * flow; + + // Function call for next group if applicable. + if constexpr (Group + 1 < num_groups) { + get_derivatives_impl(pop, y, t, dydt); + } + } + + /** + * @brief Calculates flows that are caused by people becoming infected (outflow from compartment S, Ra or Rb) for Group1. + * + * This is done recursively by calculating the interaction terms with each group. + * @tparam Group1 The group for which the derivative of the compartment should be calculated. + * @tparam Group2 The group that Group1 interacts with. + * @param[in] pop The current state of the population in the geographic unit we are considering. + * @param[in] y The current state of the model (or a subpopulation) as a flat array. + * @param[in] t The current time. + * @param[out] dydt A reference to the calculated output. + * The flow from Susceptible needs to be split into 2 parts for Exposed_1a and Exposed_1b: + * @param[out] part_a Reference to amount of flow caused by people infected with disease a. + * @param[out] part_b Reference to amount of flow caused by people infected with disease b. + * @param[in] compartment_index Index of the compartment for that the outflow will be calculated. + * @param[in] which_disease Index for which infected people cause the outflow (0 = a, 1 = b, 2 = a and b). + */ + template + void interact(Eigen::Ref> pop, Eigen::Ref> y, + ScalarType t, Eigen::Ref> dydt, double* part_a, double* part_b, + size_t compartment_index, int which_disease) const + { + static_assert((Group1 < num_groups) && (Group1 >= 0) && (Group2 < num_groups) && (Group2 >= 0), + "The template parameters Group1 & Group2 should be valid."); + using LctStateGroup2 = type_at_index_t; + ScalarType infectedNoSymptoms_2_a = 0; + ScalarType infectedSymptoms_2_a = 0; + ScalarType infectedNoSymptoms_2_b = 0; + ScalarType infectedSymptoms_2_b = 0; + auto params = this->parameters; + + size_t first_index_group2 = this->populations.template get_first_index_of_group(); + + // Calculate sum of all subcompartments for InfectedNoSymptoms for disease a of Group2. + infectedNoSymptoms_2_a = + pop.segment(first_index_group2 + + LctStateGroup2::template get_first_index(), + LctStateGroup2::template get_num_subcompartments()) + .sum() + + pop.segment(first_index_group2 + + LctStateGroup2::template get_first_index(), + LctStateGroup2::template get_num_subcompartments()) + .sum(); + // Calculate sum of all subcompartments for InfectedSymptoms for disease a of Group2. + infectedSymptoms_2_a = + pop.segment(first_index_group2 + + LctStateGroup2::template get_first_index(), + LctStateGroup2::template get_num_subcompartments()) + .sum() + + pop.segment(first_index_group2 + + LctStateGroup2::template get_first_index(), + LctStateGroup2::template get_num_subcompartments()) + .sum(); + // Calculate sum of all subcompartments for InfectedNoSymptoms for disease b of Group2. + infectedNoSymptoms_2_b = + pop.segment(first_index_group2 + + LctStateGroup2::template get_first_index(), + LctStateGroup2::template get_num_subcompartments()) + .sum() + + pop.segment(first_index_group2 + + LctStateGroup2::template get_first_index(), + LctStateGroup2::template get_num_subcompartments()) + .sum(); + // Calculate sum of all subcompartments for InfectedSymptoms for disease b of Group2. + infectedSymptoms_2_b = + pop.segment(first_index_group2 + + LctStateGroup2::template get_first_index(), + LctStateGroup2::template get_num_subcompartments()) + .sum() + + pop.segment(first_index_group2 + + LctStateGroup2::template get_first_index(), + LctStateGroup2::template get_num_subcompartments()) + .sum(); + // Size of the Subpopulation Group2 without dead people. + double N_2 = pop.segment(first_index_group2, LctStateGroup2::Count).sum(); // sum over all compartments + N_2 = N_2 - pop.segment(LctStateGroup2::template get_first_index(), 1).sum() - + pop.segment(LctStateGroup2::template get_first_index(), 1) + .sum(); // all people minus dead people + const double divN_2 = (N_2 < Limits::zero_tolerance()) ? 0.0 : 1.0 / N_2; + ScalarType season_val = 1 + params.template get() * + sin(3.141592653589793 * ((params.template get() + t) / 182.5 + 0.5)); + + if (which_disease == 0) { // disease a + dydt[compartment_index] += + -y[compartment_index] * divN_2 * season_val * + params.template get().get_cont_freq_mat().get_matrix_at(t)( + static_cast(Group1), static_cast(Group2)) * + params.template get()[Group1] * + (params.template get()[Group2] * infectedNoSymptoms_2_a + + params.template get()[Group2] * infectedSymptoms_2_a); + } + else if (which_disease == 1) { // disease b + dydt[compartment_index] += + -y[compartment_index] * divN_2 * season_val * + params.template get().get_cont_freq_mat().get_matrix_at(t)( + static_cast(Group1), static_cast(Group2)) * + params.template get()[Group1] * + (params.template get()[Group2] * infectedNoSymptoms_2_b + + params.template get()[Group2] * infectedSymptoms_2_b); + } + else if (which_disease == 2) { // both diseases + dydt[compartment_index] += + -y[compartment_index] * divN_2 * season_val * + params.template get().get_cont_freq_mat().get_matrix_at(t)( + static_cast(Group1), static_cast(Group2)) * + (params.template get()[Group1] * + (params.template get()[Group2] * infectedNoSymptoms_2_a + + params.template get()[Group2] * infectedSymptoms_2_a) + + params.template get()[Group1] * + (params.template get()[Group2] * infectedNoSymptoms_2_b + + params.template get()[Group2] * infectedSymptoms_2_b)); + } + // To split the outflow from S between E_1a and E_1b: + *part_a += params.template get()[Group1] * + (params.template get()[Group2] * infectedNoSymptoms_2_a + + params.template get()[Group2] * infectedSymptoms_2_a); + *part_b += params.template get()[Group1] * + (params.template get()[Group2] * infectedNoSymptoms_2_b + + params.template get()[Group2] * infectedSymptoms_2_b); + + if constexpr (Group2 + 1 < num_groups) { + interact(pop, y, t, dydt, part_a, part_b, compartment_index, which_disease); + } + } + + /** + * @brief Checks whether LctState of a group satisfies all constraints. + * Recursively, it checks that all groups satisfy the constraints. + * + * @tparam group The group for which the constraints should be checked. + * @return Returns true if one or more constraints are not satisfied, false otherwise. + */ + template + bool check_constraints_impl() const + { + static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid."); + using LctStateGroup = type_at_index_t; + + if (LctStateGroup::template get_num_subcompartments() != 1) { + log_warning("Constraint check: The number of subcompartments for Susceptibles of group {} should be one!", + Group); + return true; + } + if (LctStateGroup::template get_num_subcompartments() != 1 || + LctStateGroup::template get_num_subcompartments() != 1 || + LctStateGroup::template get_num_subcompartments() != 1) { + log_warning("Constraint check: The number of subcompartments for Recovered of group {} should be one!", + Group); + return true; + } + if (LctStateGroup::template get_num_subcompartments() != 1 || + LctStateGroup::template get_num_subcompartments() != 1) { + log_warning("Constraint check: The number of subcompartments for Dead of group {} should be one!", Group); + return true; + } + + if constexpr (Group == num_groups - 1) { + return false; + } + else { + return check_constraints_impl(); + } + } +}; + +} // namespace lsecir2d +} // namespace mio +#endif // LCT_SECIR_2_DISEASE_MODEL_H diff --git a/cpp/models/lct_secir_2_diseases/parameters.h b/cpp/models/lct_secir_2_diseases/parameters.h new file mode 100644 index 0000000000..fa44ddafbf --- /dev/null +++ b/cpp/models/lct_secir_2_diseases/parameters.h @@ -0,0 +1,665 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Annika Jungklaus, Lena Ploetzke +* +* Contact: Martin J. Kuehn +* +* 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 LCT_SECIR_2_DISEASES_PARAMS_H +#define LCT_SECIR_2_DISEASES_PARAMS_H + +#include "memilio/config.h" +#include "memilio/utils/parameter_set.h" +#include "memilio/utils/logging.h" +#include "memilio/utils/uncertain_value.h" +#include "memilio/epidemiology/uncertain_matrix.h" + +namespace mio +{ +namespace lsecir2d +{ + +/********************************************************* +* Define Parameters of the LCT-SECIHURD-2-DISEASES model * +**********************************************************/ + +/** + * @brief Average time spent in the Exposed compartment for disease a in day unit. + */ +struct TimeExposed_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TimeExposed_a"; + } +}; + +/** + * @brief Average time spent in the TimeInfectedNoSymptoms before developing + * symptoms or recover for disease a in day unit. + */ +struct TimeInfectedNoSymptoms_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TimeInfectedNoSymptoms_a"; + } +}; + +/** + * @brief Average time spent in the TimeInfectedSymptoms before going to hospital + * or recover for disease a in day unit. + */ +struct TimeInfectedSymptoms_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TimeInfectedNoSymptoms_a"; + } +}; + +/** + * @brief Average time being in the Hospital before treated by ICU or recover for disease a in day unit. + */ +struct TimeInfectedSevere_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TimeInfectedSevere_a"; + } +}; + +/** + * @brief Average time treated by ICU before dead or recover for disease a in day unit. + */ +struct TimeInfectedCritical_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TimeInfectedCritical_a"; + } +}; + +/** + * @brief Probability of getting infected from a contact for disease a. + */ +struct TransmissionProbabilityOnContact_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TransmissionProbabilityOnContact_a"; + } +}; + +/** + * @brief Average time spent in the Exposed compartment for disease b in day unit. + */ +struct TimeExposed_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TimeExposed_b"; + } +}; + +/** + * @brief Average time spent in the TimeInfectedNoSymptoms before developing + * symptoms or recover for disease b in day unit. + */ +struct TimeInfectedNoSymptoms_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TimeInfectedNoSymptoms_b"; + } +}; + +/** + * @brief Average time spent in the TimeInfectedSymptoms before going to hospital + * or recover for disease b in day unit. + */ +struct TimeInfectedSymptoms_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TimeInfectedNoSymptoms_b"; + } +}; + +/** + * @brief Average time being in the Hospital before treated by ICU or recover for disease b in day unit. + */ +struct TimeInfectedSevere_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TimeInfectedSevere_b"; + } +}; + +/** + * @brief Average time treated by ICU before dead or recover for disease b in day unit. + */ +struct TimeInfectedCritical_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TimeInfectedCritical_b"; + } +}; + +/** + * @brief Probability of getting infected from a contact for disease b. + */ +struct TransmissionProbabilityOnContact_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "TransmissionProbabilityOnContact_b"; + } +}; + +/** + * @brief The contact patterns within the society are modelled using an UncertainContactMatrix. + */ +struct ContactPatterns { + using Type = UncertainContactMatrix; + static Type get_default(size_t size = 1) // no age groups + { + mio::ContactMatrixGroup contact_matrix(1, (Eigen::Index)size); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant((Eigen::Index)size, (Eigen::Index)size, 10.)); + return Type(contact_matrix); + } + static std::string name() + { + return "ContactPatterns"; + } +}; + +/** + * @brief The relative InfectedNoSymptoms infectability for disease a. + */ +struct RelativeTransmissionNoSymptoms_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "RelativeTransmissionNoSymptoms_a"; + } +}; + +/** + * @brief The risk of infection from symptomatic cases for disease a. + */ +struct RiskOfInfectionFromSymptomatic_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "RiskOfInfectionFromSymptomatic_a"; + } +}; + +/** + * @brief The relative InfectedNoSymptoms infectability for disease b. + */ +struct RelativeTransmissionNoSymptoms_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "RelativeTransmissionNoSymptoms_b"; + } +}; + +/** + * @brief The risk of infection from symptomatic cases for disease b. + */ +struct RiskOfInfectionFromSymptomatic_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 1.); + } + static std::string name() + { + return "RiskOfInfectionFromSymptomatic_b"; + } +}; + +/** + * @brief The percentage of asymptomatic cases for disease a. + */ +struct RecoveredPerInfectedNoSymptoms_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 0.5); + } + static std::string name() + { + return "RecoveredPerInfectedNoSymptoms_a"; + } +}; + +/** + * @brief The percentage of hospitalized patients per infected patients for disease a. + */ +struct SeverePerInfectedSymptoms_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 0.5); + } + static std::string name() + { + return "SeverePerInfectedSymptoms_a"; + } +}; + +/** + * @brief The percentage of ICU patients per hospitalized patients for disease a. + */ +struct CriticalPerSevere_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 0.5); + } + static std::string name() + { + return "CriticalPerSevere_a"; + } +}; + +/** + * @brief The percentage of dead patients per ICU patients for disease a. + */ +struct DeathsPerCritical_a { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 0.1); + } + static std::string name() + { + return "DeathsPerCritical_a"; + } +}; + +/** + * @brief The percentage of asymptomatic cases for disease b. + */ +struct RecoveredPerInfectedNoSymptoms_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 0.5); + } + static std::string name() + { + return "RecoveredPerInfectedNoSymptoms_b"; + } +}; + +/** + * @brief The percentage of hospitalized patients per infected patients for disease b. + */ +struct SeverePerInfectedSymptoms_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 0.5); + } + static std::string name() + { + return "SeverePerInfectedSymptoms_b"; + } +}; + +/** + * @brief The percentage of ICU patients per hospitalized patients for disease b. + */ +struct CriticalPerSevere_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 0.5); + } + static std::string name() + { + return "CriticalPerSevere_b"; + } +}; + +/** + * @brief The percentage of dead patients per ICU patients for disease b. + */ +struct DeathsPerCritical_b { + using Type = Eigen::VectorX>; + static Type get_default(size_t size = 1) + { + return Type::Constant(size, 1, 0.1); + } + static std::string name() + { + return "DeathsPerCritical_b"; + } +}; + +/** + * @brief The start day in the LCT-SECIR-2-DISEASES model. + * The start day defines in which season the simulation is started. + * If the start day is 180 and simulation takes place from t0=0 to + * tmax=100 the days 180 to 280 of the year are simulated. + */ +struct StartDay { + using Type = ScalarType; + static Type get_default(size_t) + { + return 0.; + } + static std::string name() + { + return "StartDay"; + } +}; + +/** + * @brief The seasonality in the LCT-SECIR-2-DISEASES model. + * The seasonality is given as (1+k*sin()) where the sine + * curve is below one in summer and above one in winter. + */ +struct Seasonality { + using Type = ScalarType; + static Type get_default(size_t) + { + return 0.; + } + static std::string name() + { + return "Seasonality"; + } +}; + +using ParametersBase = + ParameterSet; + +/** + * @brief Parameters of an LCT-SECIR-2-DISEASES model. + */ +class Parameters : public ParametersBase +{ +public: + /** + * @brief Constructor. + * @param num_groups The number of groups considered in the LCT2D model. + */ + Parameters(size_t num_groups) + : ParametersBase(num_groups) + , m_num_groups{num_groups} + { + } + + size_t get_num_groups() const + { + return m_num_groups; + } + + /** + * @brief Checks whether all parameters satisfy their corresponding constraints and throws errors, if they do not. + * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false. + */ + bool check_constraints() const + { + for (size_t i = 0; i < m_num_groups; ++i) { + if (this->get() < 0.0 || this->get() > 0.5) { + log_warning("Constraint check: Parameter Seasonality should lie between {:0.4f} and {:.4f}", 0.0, 0.5); + return true; + } + + if (this->get()[i] < 1.0) { + log_error("Constraint check: Parameter TimeExposed_a is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get()[i] < 1.0) { + log_error("Constraint check: Parameter TimeExposed_b is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get()[i] < 1.0) { + log_error("Constraint check: Parameter TimeInfectedNoSymptoms_a is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get()[i] < 1.0) { + log_error("Constraint check: Parameter TimeInfectedNoSymptoms_b is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get()[i] < 1.0) { + log_error("Constraint check: Parameter TimeInfectedSymptoms_a is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get()[i] < 1.0) { + log_error("Constraint check: Parameter TimeInfectedSymptoms_b is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get()[i] < 1.0) { + log_error("Constraint check: Parameter TimeInfectedSevere_a is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get()[i] < 1.0) { + log_error("Constraint check: Parameter TimeInfectedSevere_b is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get()[i] < 1.0) { + log_error("Constraint check: Parameter TimeInfectedCritical_a is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get()[i] < 1.0) { + log_error("Constraint check: Parameter TimeInfectedCritical_b is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_error("Constraint check: Parameter TransmissionProbabilityOnContact_a smaller {:d} or larger {:d}", + 0, 1); + return true; + } + + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_error("Constraint check: Parameter TransmissionProbabilityOnContact_b smaller {:d} or larger {:d}", + 0, 1); + return true; + } + + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_error("Constraint check: Parameter RelativeTransmissionNoSymptoms_a smaller {:d} or larger {:d}", 0, + 1); + return true; + } + + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_error("Constraint check: Parameter RelativeTransmissionNoSymptoms_b smaller {:d} or larger {:d}", 0, + 1); + return true; + } + + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_error("Constraint check: Parameter RiskOfInfectionFromSymptomatic_a smaller {:d} or larger {:d}", + 0, 1); + return true; + } + + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_error("Constraint check: Parameter RiskOfInfectionFromSymptomatic_b smaller {:d} or larger {:d}", + 0, 1); + return true; + } + + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_error("Constraint check: Parameter RecoveredPerInfectedNoSymptoms_a smaller {:d} or larger {:d}", 0, + 1); + return true; + } + + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_error("Constraint check: Parameter RecoveredPerInfectedNoSymptoms_b smaller {:d} or larger {:d}", 0, + 1); + return true; + } + + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_error("Constraint check: Parameter SeverePerInfectedSymptoms_a smaller {:d} or larger {:d}", 0, 1); + return true; + } + + if (this->get()[i] < 0.0 || + this->get()[i] > 1.0) { + log_error("Constraint check: Parameter SeverePerInfectedSymptoms_b smaller {:d} or larger {:d}", 0, 1); + return true; + } + + if (this->get()[i] < 0.0 || this->get()[i] > 1.0) { + log_error("Constraint check: Parameter CriticalPerSevere_a smaller {:d} or larger {:d}", 0, 1); + return true; + } + + if (this->get()[i] < 0.0 || this->get()[i] > 1.0) { + log_error("Constraint check: Parameter CriticalPerSevere_b smaller {:d} or larger {:d}", 0, 1); + return true; + } + + if (this->get()[i] < 0.0 || this->get()[i] > 1.0) { + log_error("Constraint check: Parameter DeathsPerCritical_a smaller {:d} or larger {:d}", 0, 1); + return true; + } + + if (this->get()[i] < 0.0 || this->get()[i] > 1.0) { + log_error("Constraint check: Parameter DeathsPerCritical_b smaller {:d} or larger {:d}", 0, 1); + return true; + } + } + + return false; + } + +private: + Parameters(ParametersBase&& base) + : ParametersBase(std::move(base)) + , m_num_groups(this->template get().get_cont_freq_mat().get_num_groups()) + { + } + + size_t m_num_groups; + +public: + /** + * deserialize an object of this class. + * @see mio::deserialize + */ + template + static IOResult deserialize(IOContext& io) + { + BOOST_OUTCOME_TRY(auto&& base, ParametersBase::deserialize(io)); + return success(Parameters(std::move(base))); + } +}; + +} // namespace lsecir2d +} // namespace mio + +#endif // LCT_SECIR_2_DISEASES_PARAMS_H diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 126e17ff2c..2016620721 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -70,6 +70,7 @@ set(TESTSOURCES test_ide_secir_ageres.cpp test_state_age_function.cpp test_lct_secir.cpp + test_lct_secir_2_diseases.cpp test_lct_initializer_flows.cpp test_glct_secir.cpp test_ad.cpp diff --git a/cpp/tests/test_lct_secir_2_diseases.cpp b/cpp/tests/test_lct_secir_2_diseases.cpp new file mode 100644 index 0000000000..ed5d8b8f38 --- /dev/null +++ b/cpp/tests/test_lct_secir_2_diseases.cpp @@ -0,0 +1,956 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Annika Jungklaus, Lena Ploetzke +* +* Contact: Martin J. Kuehn +* +* 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 "lct_secir_2_diseases/parameters.h" +#include "lct_secir/model.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "memilio/utils/logging.h" +#include "memilio/epidemiology/contact_matrix.h" +#include "memilio/data/analyze_result.h" +#include "memilio/compartments/simulation.h" +#include "load_test_data.h" +#include +#include +#include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp" + +// Test confirms that default construction of an LCT2D model works. +TEST(TestLCTSecir2d, simulateDefault) +{ + using InfState = mio::lsecir2d::InfectionState; + using LctState = + mio::LctInfectionState; + using Model = mio::lsecir2d::Model; + ScalarType t0 = 0; + ScalarType tmax = 1; + ScalarType dt = 0.1; + + Eigen::VectorX init = Eigen::VectorX::Constant((Eigen::Index)InfState::Count, 0); + init[0] = 200; + init[3] = 50; // people infected with disease a + init[5] = 30; + init[15] = 50; // people infected with disease b + init[17] = 30; + + Model model; + for (size_t i = 0; i < LctState::Count; i++) { + model.populations[i] = init[i]; + } + + mio::TimeSeries result = mio::simulate(t0, tmax, dt, model); + + EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); + ScalarType sum_pop = init.sum(); + for (Eigen::Index i = 0; i < result.get_num_time_points(); i++) { + EXPECT_NEAR(sum_pop, result[i].sum(), 1e-5); // check that total pop. is constant + } +} + +/* Tests comparing the result for an LCT SECIR 2 DISEASE model (with transmission prob. 0 for one disease and 2 age groups) + with the result of the equivalent LCT SECIR model. */ +// 1. Test: First infection for disease a (transmission prob. of disease b = 0) +TEST(TestLCTSecir2d, compareWithLCTSecir1) +{ + using InfState2d = mio::lsecir2d::InfectionState; + using LctState2d = mio::LctInfectionState; + using Model_2d = mio::lsecir2d::Model; + + using InfState_lct = mio::lsecir::InfectionState; + using LctState_lct = mio::LctInfectionState; + using Model_lct = mio::lsecir::Model; + + ScalarType t0 = 0; + ScalarType tmax = 5; + ScalarType dt = 0.1; + + // Initialization vector for lct2d model. + Eigen::VectorX init_lct2d = Eigen::VectorX::Constant((Eigen::Index)InfState2d::Count, 0); + init_lct2d[0] = 200; // lct and lct2d use different infection states + init_lct2d[3] = 50; // make sure initial pop. is in the same compartments for lct and lct2d + init_lct2d[5] = 30; + + // Define LCT2D model. + Model_2d model_lct2d; + //Set initial values + for (size_t i = 0; i < LctState2d::Count; i++) { + model_lct2d.populations[i] = init_lct2d[i]; + } + + // Set Parameters. + model_lct2d.parameters.get()[0] = 3.2; + model_lct2d.parameters.get()[0] = 2; + model_lct2d.parameters.get()[0] = 5.8; + model_lct2d.parameters.get()[0] = 9.5; + model_lct2d.parameters.get()[0] = 7.1; + model_lct2d.parameters.get()[0] = 3.2; + model_lct2d.parameters.get()[0] = 2; + model_lct2d.parameters.get()[0] = 5.8; + model_lct2d.parameters.get()[0] = 9.5; + model_lct2d.parameters.get()[0] = 7.1; + + model_lct2d.parameters.get()[0] = 0.05; + model_lct2d.parameters.get()[0] = 0.; + + mio::ContactMatrixGroup& contact_matrix_lct2d = model_lct2d.parameters.get(); + contact_matrix_lct2d[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(2, 2, 10)); + contact_matrix_lct2d[0].add_damping(0.7, mio::SimulationTime(2.)); + + model_lct2d.parameters.get()[0] = 0.7; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.09; + model_lct2d.parameters.get()[0] = 0.2; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.3; + model_lct2d.parameters.get()[0] = 0.7; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.09; + model_lct2d.parameters.get()[0] = 0.2; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.3; + model_lct2d.parameters.get() = 50; + model_lct2d.parameters.get() = 0.1; + + // Initialization vector for LCT model. + Eigen::VectorX init_lct = Eigen::VectorX::Constant((Eigen::Index)InfState_lct::Count, 0); + init_lct[0] = 200; + init_lct[3] = 50; + init_lct[5] = 30; + + // Define LCT model. + Model_lct model_lct; + //Set initial values + for (size_t i = 0; i < LctState_lct::Count; i++) { + model_lct.populations[i] = init_lct[i]; + } + + // Set Parameters. + model_lct.parameters.get()[0] = 3.2; + model_lct.parameters.get()[0] = 2; + model_lct.parameters.get()[0] = 5.8; + model_lct.parameters.get()[0] = 9.5; + model_lct.parameters.get()[0] = 7.1; + + model_lct.parameters.get()[0] = 0.05; + + mio::ContactMatrixGroup& contact_matrix_lct = model_lct.parameters.get(); + contact_matrix_lct[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(2, 2, 10)); + contact_matrix_lct[0].add_damping(0.7, mio::SimulationTime(2.)); + + model_lct.parameters.get()[0] = 0.7; + model_lct.parameters.get()[0] = 0.25; + model_lct.parameters.get() = 50; + model_lct.parameters.get() = 0.1; + model_lct.parameters.get()[0] = 0.09; + model_lct.parameters.get()[0] = 0.2; + model_lct.parameters.get()[0] = 0.25; + model_lct.parameters.get()[0] = 0.3; + + // Simulate + mio::TimeSeries result_lct2d = mio::simulate(t0, tmax, dt, model_lct2d); + + mio::TimeSeries result_lct = mio::simulate(t0, tmax, dt, model_lct); + + // Simulation results should be equal. + // Compare LCT with infection 1a in LCT2D + ASSERT_EQ(result_lct.get_num_time_points(), result_lct2d.get_num_time_points()); + for (int i = 0; i < 4; ++i) { + EXPECT_NEAR(result_lct.get_time(i), result_lct2d.get_time(i), 1e-5); + + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::Susceptible], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Susceptible], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::Exposed], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Exposed_1a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::InfectedNoSymptoms], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedNoSymptoms_1a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::InfectedSymptoms], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedSymptoms_1a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::InfectedCritical], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedCritical_1a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::InfectedSevere], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedSevere_1a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::Recovered], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Recovered_a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::Dead], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Dead_a], 1e-5); + } +} + +// 2. Test: First infection with disease b (transmission prob. of disease a = 0) +TEST(TestLCTSecir2d, compareWithLCTSecir2) +{ + using InfState2d = mio::lsecir2d::InfectionState; + using LctState2d = mio::LctInfectionState; + using Model_2d = mio::lsecir2d::Model; + ScalarType t0 = 0; + ScalarType tmax = 5; + ScalarType dt = 0.1; + + // Initialization vector for lct2d model. + Eigen::VectorX init_lct2d = Eigen::VectorX::Constant((Eigen::Index)InfState2d::Count, 0); + init_lct2d[0] = 200; // lct and lct2d use different infection states + init_lct2d[15] = 50; // make sure initial pop. is in the same compartments for lct and lct2d + init_lct2d[17] = 30; + + // Define LCT2D model. + Model_2d model_lct2d; + //Set initial values + for (size_t i = 0; i < LctState2d::Count; i++) { + model_lct2d.populations[i] = init_lct2d[i]; + } + + // Set Parameters. + model_lct2d.parameters.get()[0] = 3.2; + model_lct2d.parameters.get()[0] = 2; + model_lct2d.parameters.get()[0] = 5.8; + model_lct2d.parameters.get()[0] = 9.5; + model_lct2d.parameters.get()[0] = 7.1; + model_lct2d.parameters.get()[0] = 3.2; + model_lct2d.parameters.get()[0] = 2; + model_lct2d.parameters.get()[0] = 5.8; + model_lct2d.parameters.get()[0] = 9.5; + model_lct2d.parameters.get()[0] = 7.1; + + model_lct2d.parameters.get()[0] = 0.; + model_lct2d.parameters.get()[0] = 0.05; + + mio::ContactMatrixGroup& contact_matrix_lct2d = model_lct2d.parameters.get(); + contact_matrix_lct2d[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(2, 2, 10)); + contact_matrix_lct2d[0].add_damping(0.7, mio::SimulationTime(2.)); + + model_lct2d.parameters.get()[0] = 0.7; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.09; + model_lct2d.parameters.get()[0] = 0.2; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.3; + model_lct2d.parameters.get()[0] = 0.7; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.09; + model_lct2d.parameters.get()[0] = 0.2; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.3; + model_lct2d.parameters.get() = 50; + model_lct2d.parameters.get() = 0.1; + + using InfState = mio::lsecir::InfectionState; + using LctState = mio::LctInfectionState; + using Model = mio::lsecir::Model; + + // Initialization vector for LCT model. + Eigen::VectorX init_lct = Eigen::VectorX::Constant((Eigen::Index)InfState::Count, 0); + init_lct[0] = 200; + init_lct[3] = 50; + init_lct[5] = 30; + + // Define LCT model. + Model model_lct; + //Set initial values + for (size_t i = 0; i < LctState::Count; i++) { + model_lct.populations[i] = init_lct[i]; + } + + // Set Parameters. + model_lct.parameters.get()[0] = 3.2; + model_lct.parameters.get()[0] = 2; + model_lct.parameters.get()[0] = 5.8; + model_lct.parameters.get()[0] = 9.5; + model_lct.parameters.get()[0] = 7.1; + + model_lct.parameters.get()[0] = 0.05; + + mio::ContactMatrixGroup& contact_matrix_lct = model_lct.parameters.get(); + contact_matrix_lct[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(2, 2, 10)); + contact_matrix_lct[0].add_damping(0.7, mio::SimulationTime(2.)); + + model_lct.parameters.get()[0] = 0.7; + model_lct.parameters.get()[0] = 0.25; + model_lct.parameters.get() = 50; + model_lct.parameters.get() = 0.1; + model_lct.parameters.get()[0] = 0.09; + model_lct.parameters.get()[0] = 0.2; + model_lct.parameters.get()[0] = 0.25; + model_lct.parameters.get()[0] = 0.3; + + // Simulate + mio::TimeSeries result_lct2d = mio::simulate(t0, tmax, dt, model_lct2d); + + mio::TimeSeries result_lct = mio::simulate(t0, tmax, dt, model_lct); + + // Simulation results should be equal. + // Compare LCT with Infection 1b in LCT2D + ASSERT_EQ(result_lct.get_num_time_points(), result_lct2d.get_num_time_points()); + for (int i = 0; i < 4; ++i) { + EXPECT_NEAR(result_lct.get_time(i), result_lct2d.get_time(i), 1e-5); + + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::Susceptible], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Susceptible], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::Exposed], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Exposed_1b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::InfectedNoSymptoms], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedNoSymptoms_1b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::InfectedSymptoms], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedSymptoms_1b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::InfectedCritical], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedCritical_1b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::InfectedSevere], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedSevere_1b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::Recovered], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Recovered_b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::Dead], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Dead_b], 1e-5); + } +} + +// 3. Test: Second infection with disease a (transmission prob. of disease b = 0, start in Recovered_b)*/ +TEST(TestLCTSecir2d, compareWithLCTSecir3) +{ + using InfState_2d = mio::lsecir2d::InfectionState; + using LctState_2d = mio::LctInfectionState; + using Model_2d = mio::lsecir2d::Model; + + using InfState_lct = mio::lsecir::InfectionState; + using LctState_lct = mio::LctInfectionState; + using Model_lct = mio::lsecir::Model; + + ScalarType t0 = 0; + ScalarType tmax = 5; + ScalarType dt = 0.1; + + // Initialization vector for lct2d model. + Eigen::VectorX init_lct2d = Eigen::VectorX::Constant((Eigen::Index)InfState_2d::Count, 0); + init_lct2d[18] = 200; // lct and lct2d use different infection states + init_lct2d[10] = 50; // make sure initial pop. is in the same compartments for lct and lct2d + init_lct2d[12] = 30; + + // Define LCT2D model. + Model_2d model_lct2d; + //Set initial values + for (size_t i = 0; i < LctState_2d::Count; i++) { + model_lct2d.populations[i] = init_lct2d[i]; + } + + // Set Parameters. + model_lct2d.parameters.get()[0] = 3.2; + model_lct2d.parameters.get()[0] = 2; + model_lct2d.parameters.get()[0] = 5.8; + model_lct2d.parameters.get()[0] = 9.5; + model_lct2d.parameters.get()[0] = 7.1; + model_lct2d.parameters.get()[0] = 3.2; + model_lct2d.parameters.get()[0] = 2; + model_lct2d.parameters.get()[0] = 5.8; + model_lct2d.parameters.get()[0] = 9.5; + model_lct2d.parameters.get()[0] = 7.1; + + model_lct2d.parameters.get()[0] = 0.05; + model_lct2d.parameters.get()[0] = 0.; + + mio::ContactMatrixGroup& contact_matrix_lct2d = model_lct2d.parameters.get(); + contact_matrix_lct2d[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(2, 2, 10)); + contact_matrix_lct2d[0].add_damping(0.7, mio::SimulationTime(2.)); + + model_lct2d.parameters.get()[0] = 0.7; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.09; + model_lct2d.parameters.get()[0] = 0.2; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.3; + model_lct2d.parameters.get()[0] = 0.7; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.09; + model_lct2d.parameters.get()[0] = 0.2; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.3; + model_lct2d.parameters.get() = 50; + model_lct2d.parameters.get() = 0.1; + + // Initialization vector for LCT model. + Eigen::VectorX init_lct = Eigen::VectorX::Constant((Eigen::Index)InfState_lct::Count, 0); + init_lct[0] = 200; + init_lct[3] = 50; + init_lct[5] = 30; + + // Define LCT model. + Model_lct model_lct; + //Set initial values + for (size_t i = 0; i < LctState_lct::Count; i++) { + model_lct.populations[i] = init_lct[i]; + } + + // Set Parameters. + model_lct.parameters.get()[0] = 3.2; + model_lct.parameters.get()[0] = 2; + model_lct.parameters.get()[0] = 5.8; + model_lct.parameters.get()[0] = 9.5; + model_lct.parameters.get()[0] = 7.1; + + model_lct.parameters.get()[0] = 0.05; + + mio::ContactMatrixGroup& contact_matrix_lct = model_lct.parameters.get(); + contact_matrix_lct[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(2, 2, 10)); + contact_matrix_lct[0].add_damping(0.7, mio::SimulationTime(2.)); + + model_lct.parameters.get()[0] = 0.7; + model_lct.parameters.get()[0] = 0.25; + model_lct.parameters.get() = 50; + model_lct.parameters.get() = 0.1; + model_lct.parameters.get()[0] = 0.09; + model_lct.parameters.get()[0] = 0.2; + model_lct.parameters.get()[0] = 0.25; + model_lct.parameters.get()[0] = 0.3; + + // Simulate + mio::TimeSeries result_lct2d = mio::simulate(t0, tmax, dt, model_lct2d); + + mio::TimeSeries result_lct = mio::simulate(t0, tmax, dt, model_lct); + + // Simulation results should be equal. + // Compare LCT with Infection 2a in LCT2D + ASSERT_EQ(result_lct.get_num_time_points(), result_lct2d.get_num_time_points()); + for (int i = 0; i < 4; ++i) { + EXPECT_NEAR(result_lct.get_time(i), result_lct2d.get_time(i), 1e-5); + + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::Susceptible], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Recovered_b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::Exposed], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Exposed_2a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::InfectedNoSymptoms], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedNoSymptoms_2a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::InfectedSymptoms], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedSymptoms_2a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::InfectedCritical], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedCritical_2a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::InfectedSevere], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedSevere_2a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::Recovered], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Recovered_ab], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState_lct::Dead], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Dead_a], 1e-5); + } +} + +// 4. Test: Second infection with disease b (transmission prob. of disease a = 0, start in Recovered_a)*/ +TEST(TestLCTSecir2d, compareWithLCTSecir4) +{ + using InfState2d = mio::lsecir2d::InfectionState; + using LctState2d = mio::LctInfectionState; + using Model_2d = mio::lsecir2d::Model; + ScalarType t0 = 0; + ScalarType tmax = 5; + ScalarType dt = 0.1; + + // Initialization vector for lct2d model. + Eigen::VectorX init_lct2d = Eigen::VectorX::Constant((Eigen::Index)InfState2d::Count, 0); + init_lct2d[6] = 200; // lct and lct2d use different infection states + init_lct2d[22] = 50; // make sure initial pop. is in the same compartments for lct and lct2d + init_lct2d[24] = 30; + + // Define LCT2D model. + Model_2d model_lct2d; + //Set initial values + for (size_t i = 0; i < LctState2d::Count; i++) { + model_lct2d.populations[i] = init_lct2d[i]; + } + + // Set Parameters. + model_lct2d.parameters.get()[0] = 3.2; + model_lct2d.parameters.get()[0] = 2; + model_lct2d.parameters.get()[0] = 5.8; + model_lct2d.parameters.get()[0] = 9.5; + model_lct2d.parameters.get()[0] = 7.1; + model_lct2d.parameters.get()[0] = 3.2; + model_lct2d.parameters.get()[0] = 2; + model_lct2d.parameters.get()[0] = 5.8; + model_lct2d.parameters.get()[0] = 9.5; + model_lct2d.parameters.get()[0] = 7.1; + + model_lct2d.parameters.get()[0] = 0.; + model_lct2d.parameters.get()[0] = 0.05; + + mio::ContactMatrixGroup& contact_matrix_lct2d = model_lct2d.parameters.get(); + contact_matrix_lct2d[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(2, 2, 10)); + contact_matrix_lct2d[0].add_damping(0.7, mio::SimulationTime(2.)); + + model_lct2d.parameters.get()[0] = 0.7; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.09; + model_lct2d.parameters.get()[0] = 0.2; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.3; + model_lct2d.parameters.get()[0] = 0.7; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.09; + model_lct2d.parameters.get()[0] = 0.2; + model_lct2d.parameters.get()[0] = 0.25; + model_lct2d.parameters.get()[0] = 0.3; + model_lct2d.parameters.get() = 50; + model_lct2d.parameters.get() = 0.1; + + using InfState = mio::lsecir::InfectionState; + using LctState = mio::LctInfectionState; + using Model = mio::lsecir::Model; + + // Initialization vector for LCT model. + Eigen::VectorX init_lct = Eigen::VectorX::Constant((Eigen::Index)InfState::Count, 0); + init_lct[0] = 200; + init_lct[3] = 50; + init_lct[5] = 30; + + // Define LCT model. + Model model_lct; + //Set initial values + for (size_t i = 0; i < LctState::Count; i++) { + model_lct.populations[i] = init_lct[i]; + } + + // Set Parameters. + model_lct.parameters.get()[0] = 3.2; + model_lct.parameters.get()[0] = 2; + model_lct.parameters.get()[0] = 5.8; + model_lct.parameters.get()[0] = 9.5; + model_lct.parameters.get()[0] = 7.1; + + model_lct.parameters.get()[0] = 0.05; + + mio::ContactMatrixGroup& contact_matrix_lct = model_lct.parameters.get(); + contact_matrix_lct[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(2, 2, 10)); + contact_matrix_lct[0].add_damping(0.7, mio::SimulationTime(2.)); + + model_lct.parameters.get()[0] = 0.7; + model_lct.parameters.get()[0] = 0.25; + model_lct.parameters.get() = 50; + model_lct.parameters.get() = 0.1; + model_lct.parameters.get()[0] = 0.09; + model_lct.parameters.get()[0] = 0.2; + model_lct.parameters.get()[0] = 0.25; + model_lct.parameters.get()[0] = 0.3; + + // Simulate + mio::TimeSeries result_lct2d = mio::simulate(t0, tmax, dt, model_lct2d); + + mio::TimeSeries result_lct = mio::simulate(t0, tmax, dt, model_lct); + + // Simulation results should be equal. + // Compare LCT with Infection 2b in LCT2D + ASSERT_EQ(result_lct.get_num_time_points(), result_lct2d.get_num_time_points()); + for (int i = 0; i < 4; ++i) { + EXPECT_NEAR(result_lct.get_time(i), result_lct2d.get_time(i), 1e-5); + + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::Susceptible], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Recovered_a], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::Exposed], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Exposed_2b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::InfectedNoSymptoms], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedNoSymptoms_2b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::InfectedSymptoms], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedSymptoms_2b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::InfectedCritical], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedCritical_2b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::InfectedSevere], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::InfectedSevere_2b], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::Recovered], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Recovered_ab], 1e-5); + EXPECT_NEAR(result_lct[i][(Eigen::Index)InfState::Dead], + result_lct2d[i][(Eigen::Index)mio::lsecir2d::InfectionState::Dead_b], 1e-5); + } +} + +// Run the model with more than one subcompartment for states E,I,C,H,U +// and calculate the TimeSeries with no subcompartments from the result +TEST(TestLCTSecir2d, testSubcompartments) +{ + using InfState = mio::lsecir2d::InfectionState; + using LctState = + mio::LctInfectionState; + using Model = mio::lsecir2d::Model; + ScalarType t0 = 0; + ScalarType tmax = 1; + ScalarType dt = 0.1; + + // Initial population, split into subcompartments + std::vector> init = { + {200}, {0, 0}, {30, 10, 10}, {0, 0, 0}, {10, 10, 10}, {0, 0, 0}, {0}, {0}, {0, 0}, + {30, 10, 10}, {0, 0, 0}, {10, 10, 10}, {0, 0, 0}, {0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, + {0}, {0}, {0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0}}; + + Model model; + + // Transfer the initial values in initial_populations to the model. + std::vector flat_init; + for (auto&& vec : init) { + flat_init.insert(flat_init.end(), vec.begin(), vec.end()); + } + for (size_t i = 0; i < LctState::Count; i++) { + model.populations[i] = flat_init[i]; + } + + mio::TimeSeries result = mio::simulate(t0, tmax, dt, model); + mio::TimeSeries population_no_subcompartments = model.calculate_compartments(result); + auto interpolated_results = mio::interpolate_simulation_result(population_no_subcompartments); + + EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); +} + +// Model setup to compare result with a previous output. +class ModelTestLCTSecir2d : public testing::Test +{ +public: + using InfState = mio::lsecir2d::InfectionState; + using LctState = + mio::LctInfectionState; + using Model = mio::lsecir2d::Model; + +protected: + virtual void SetUp() + { + // Define initial distribution of the population in the subcompartments. + std::vector> initial_populations = {{200}, {0}, {0}, {30}, {0}, {0}, {0}, {0}, {0}, + {0}, {0}, {0}, {0}, {0}, {0}, {30}, {0}, {0}, + {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}}; + model = new Model(); + // Transfer the initial values in initial_populations to the model. + std::vector 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]; + } + + // Set parameters. + model->parameters.get()[0] = 3.2; + model->parameters.get()[0] = 2; + model->parameters.get()[0] = 5.8; + model->parameters.get()[0] = 9.5; + model->parameters.get()[0] = 7.1; + model->parameters.get()[0] = 0.05; + model->parameters.get()[0] = 3.2; + model->parameters.get()[0] = 2; + model->parameters.get()[0] = 5.8; + model->parameters.get()[0] = 9.5; + model->parameters.get()[0] = 7.1; + model->parameters.get()[0] = 0.05; + + mio::ContactMatrixGroup& contact_matrix = model->parameters.get(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10)); + contact_matrix[0].add_damping(0.7, mio::SimulationTime(2.)); + + model->parameters.get()[0] = 0.7; + model->parameters.get()[0] = 0.25; + model->parameters.get()[0] = 0.09; + model->parameters.get()[0] = 0.2; + model->parameters.get()[0] = 0.25; + model->parameters.get()[0] = 0.3; + model->parameters.get()[0] = 0.7; + model->parameters.get()[0] = 0.25; + model->parameters.get()[0] = 0.09; + model->parameters.get()[0] = 0.2; + model->parameters.get()[0] = 0.25; + model->parameters.get()[0] = 0.3; + } + + virtual void TearDown() + { + delete model; + } + +public: + Model* model = nullptr; +}; + +// Test calculate_compartments with a TimeSeries that has an incorrect number of elements. +TEST_F(ModelTestLCTSecir2d, testCalculatePopWrongSize) +{ + // Deactivate temporarily log output because an error is expected. + mio::set_log_level(mio::LogLevel::off); + // TimeSeries has to have LctState::Count elements. + size_t wrong_size = LctState::Count - 2; + // Define TimeSeries with wrong_size elements. + mio::TimeSeries wrong_num_elements(wrong_size); + Eigen::VectorX vec_wrong_size = Eigen::VectorX::Ones(wrong_size); + wrong_num_elements.add_time_point(-10, vec_wrong_size); + wrong_num_elements.add_time_point(-9, vec_wrong_size); + // Call the calculate_compartments function with the TimeSeries with a wrong number of elements. + mio::TimeSeries population = model->calculate_compartments(wrong_num_elements); + // A TimeSeries of the right size with values -1 is expected. + ASSERT_EQ(1, population.get_num_time_points()); + for (int i = 0; i < population.get_num_elements(); i++) { + EXPECT_EQ(-1, population.get_last_value()[i]); + } + // Reactive log output. + mio::set_log_level(mio::LogLevel::warn); +} + +//Check constraints of Parameters class. +TEST(TestLCTSecir2d, testConstraintsParameters) +{ + // Deactivate temporarily log output for next tests. + mio::set_log_level(mio::LogLevel::off); + + // Check for exceptions of parameters. + mio::lsecir2d::Parameters parameters_lct2d(1); + parameters_lct2d.get()[0] = 0; + parameters_lct2d.get()[0] = 3.1; + parameters_lct2d.get()[0] = 6.1; + parameters_lct2d.get()[0] = 11.1; + parameters_lct2d.get()[0] = 17.1; + parameters_lct2d.get()[0] = 0.01; + parameters_lct2d.get()[0] = 3.1; + parameters_lct2d.get()[0] = 3.1; + parameters_lct2d.get()[0] = 6.1; + parameters_lct2d.get()[0] = 11.1; + parameters_lct2d.get()[0] = 17.1; + parameters_lct2d.get()[0] = 0.01; + mio::ContactMatrixGroup& contact_matrix = parameters_lct2d.get(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10)); + + parameters_lct2d.get()[0] = 1; + parameters_lct2d.get()[0] = 1; + parameters_lct2d.get()[0] = 0.1; + parameters_lct2d.get()[0] = 0.1; + parameters_lct2d.get()[0] = 0.1; + parameters_lct2d.get()[0] = 0.1; + parameters_lct2d.get()[0] = 1; + parameters_lct2d.get()[0] = 1; + parameters_lct2d.get()[0] = 0.1; + parameters_lct2d.get()[0] = 0.1; + parameters_lct2d.get()[0] = 0.1; + parameters_lct2d.get()[0] = 0.1; + + // Check improper TimeExposed. + bool constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 3.1; + + parameters_lct2d.get()[0] = 0.1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 3.1; + + // Check TimeInfectedNoSymptoms. + parameters_lct2d.get()[0] = 0.1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 3.1; + + parameters_lct2d.get()[0] = 0.1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 3.1; + + // Check TimeInfectedSymptoms. + parameters_lct2d.get()[0] = -0.1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 6.1; + + parameters_lct2d.get()[0] = -0.1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 6.1; + + // Check TimeInfectedSevere. + parameters_lct2d.get()[0] = 0.5; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 11.1; + + parameters_lct2d.get()[0] = 0.5; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 11.1; + + // Check TimeInfectedCritical. + parameters_lct2d.get()[0] = 0.; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 17.1; + + parameters_lct2d.get()[0] = 0.; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 17.1; + + // Check TransmissionProbabilityOnContact. + parameters_lct2d.get()[0] = -1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 0.01; + + parameters_lct2d.get()[0] = -1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 0.01; + + // Check RelativeTransmissionNoSymptoms. + parameters_lct2d.get()[0] = 1.5; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 1; + + parameters_lct2d.get()[0] = 1.5; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 1; + + // Check RiskOfInfectionFromSymptomatic. + parameters_lct2d.get()[0] = 1.5; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 1; + + parameters_lct2d.get()[0] = 1.5; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 1; + + // Check RecoveredPerInfectedNoSymptoms. + parameters_lct2d.get()[0] = 1.5; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 0.1; + + parameters_lct2d.get()[0] = 1.5; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 0.1; + + // Check SeverePerInfectedSymptoms. + parameters_lct2d.get()[0] = -1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 0.1; + + parameters_lct2d.get()[0] = -1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 0.1; + + // Check CriticalPerSevere. + parameters_lct2d.get()[0] = -1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 0.1; + + parameters_lct2d.get()[0] = -1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 0.1; + + // Check DeathsPerCritical. + parameters_lct2d.get()[0] = -1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 0.1; + + parameters_lct2d.get()[0] = -1; + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.get()[0] = 0.1; + + // Check Seasonality. + parameters_lct2d.set(1); + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct2d.set(0.1); + + // Check with correct parameters. + constraint_check = parameters_lct2d.check_constraints(); + EXPECT_FALSE(constraint_check); + + // Reactive log output. + mio::set_log_level(mio::LogLevel::warn); +} + +// Check constraints of the Model setup. +TEST(TestLCTSecir2d, testConstraintsModel) +{ + // Deactivate temporarily log output for next tests. + mio::set_log_level(mio::LogLevel::off); + + using InfState = mio::lsecir2d::InfectionState; + + // Check for improper number of subcompartments for Susceptible. + using LctStatewrongSusceptibles = + mio::LctInfectionState; + using ModelwrongSusceptibles = mio::lsecir2d::Model; + ModelwrongSusceptibles modelwrongSusceptibles; + bool constraint_check = modelwrongSusceptibles.check_constraints(); + EXPECT_TRUE(constraint_check); + + // Check for improper number of subcompartments for Recovered. + using LctStatewrongRecovered_a = + mio::LctInfectionState; + using ModelwrongRecovered_a = mio::lsecir2d::Model; + ModelwrongRecovered_a modelwrongRecovered_a; + constraint_check = modelwrongRecovered_a.check_constraints(); + EXPECT_TRUE(constraint_check); + + using LctStatewrongRecovered_b = + mio::LctInfectionState; + using ModelwrongRecovered_b = mio::lsecir2d::Model; + ModelwrongRecovered_b modelwrongRecovered_b; + constraint_check = modelwrongRecovered_b.check_constraints(); + EXPECT_TRUE(constraint_check); + + using LctStatewrongRecovered_ab = + mio::LctInfectionState; + using ModelwrongRecovered_ab = mio::lsecir2d::Model; + ModelwrongRecovered_ab modelwrongRecovered_ab; + constraint_check = modelwrongRecovered_ab.check_constraints(); + EXPECT_TRUE(constraint_check); + + // Check for improper number of subcompartments for Dead. + using LctStatewrongDead_a = + mio::LctInfectionState; + using ModelwrongDead_a = mio::lsecir2d::Model; + ModelwrongDead_a modelwrongDead_a; + constraint_check = modelwrongDead_a.check_constraints(); + EXPECT_TRUE(constraint_check); + + using LctStatewrongDead_b = + mio::LctInfectionState; + using ModelwrongDead_b = mio::lsecir2d::Model; + ModelwrongDead_b modelwrongDead_b; + constraint_check = modelwrongDead_b.check_constraints(); + EXPECT_TRUE(constraint_check); + + // Check with a negative number in the initial population distribution. + using LctStatevalid = + mio::LctInfectionState; + using Model = mio::lsecir2d::Model; + Model model; + model.populations[0] = -1000; + constraint_check = model.check_constraints(); + EXPECT_TRUE(constraint_check); + + // Reactive log output. + mio::set_log_level(mio::LogLevel::warn); + + // Check for valid Setup. + model.populations[0] = 1000; + constraint_check = model.check_constraints(); + EXPECT_FALSE(constraint_check); +} diff --git a/docs/source/cpp/lct.rst b/docs/source/cpp/lct.rst index 56cf4a8a98..807cae7b93 100644 --- a/docs/source/cpp/lct.rst +++ b/docs/source/cpp/lct.rst @@ -114,3 +114,4 @@ List of models :titlesonly: models/lsecir + models/lsecir2d diff --git a/docs/source/cpp/models/lsecir.rst b/docs/source/cpp/models/lsecir.rst index b818ffbef7..c7f3deefb7 100644 --- a/docs/source/cpp/models/lsecir.rst +++ b/docs/source/cpp/models/lsecir.rst @@ -1,7 +1,7 @@ LCT-based SECIR-type model ========================== -The LCT-SECIR module models and simulates an epidemic using an ODE-based approach while making use of the Linear Chain +The LCT-SECIR module models and simulates an epidemic using an ODE-based approach while making use of the Linear Chain Trick. This provides the option of Erlang distributed stay times in the compartments through the use of subcompartments. The model is particularly suited for pathogens with pre- or asymptomatic infection states and when severe or critical states are possible. The model assumes perfect immunity after recovery and is thus only suited for epidemic use cases. diff --git a/docs/source/cpp/models/lsecir2d.rst b/docs/source/cpp/models/lsecir2d.rst new file mode 100644 index 0000000000..3dbee34a62 --- /dev/null +++ b/docs/source/cpp/models/lsecir2d.rst @@ -0,0 +1,466 @@ +2 Diseases in LCT-based SECIR-type model +========================== + +| The LCT-SECIR-2-DISEASES model is an extension of the :doc:`model with one disease `. +| The model is ODE-based and uses the Linear Chain Trick to allow for more general Erlang distributed stay times in each compartment instead of just exponentially distributed stay times induced by basic ODE-based models. +| With the SECIR structure the model is particularly suited for pathogens with pre- or asymptomatic infection states and when severe or critical states are possible. +| For the two diseases :math:`a` and :math:`b` the model assumes no co-infection, a certain independence in the sense that prior infection with one disease does not affect the infection with the other disease (e.g. probability to get infected, time spend in each state, chances of recovery etc.), and perfect immunity after recovery for both diseases. + + +There are two possibilities for a susceptible individual (since we assume no co-infection): + 1. Get infected with disease a, then (if not dead) get infected with disease b + 2. Get infected with disease b, then (if not dead) get infected with disease a + +Each infection is simulated using a LCT-SECIR model. The states are labeled according to infection (first or second) and disease (:math:`a` or :math:`b`), +so the full model is given by the combination of infections :math:`1a`, :math:`2a`, :math:`1b`, and :math:`2b`. + +Below is a visualization of the infection states split into LCT-states and transitions without a stratification according to socisdemographic groups. + +.. image:: "" + :alt: tikz_lct-2d + +With infection states for :math:`i \in \{1,2\}, x \in \{a,b\}`: + +.. list-table:: + :header-rows: 1 + :widths: 20 20 60 + + * - Letter + - SECIR-State + - In the code + * - S + - `Susceptible` + - Susceptible + * - E + - `Exposed` + - Exposed_ix + * - C + - `Carrier` + - InfectedNoSymptoms_ix + * - I + - `Infected` + - InfectedSymptoms_ix + * - H + - `Hospitalized` + - InfectedSevere_ix + * - U + - `in Intensive Care Unit` + - InfectedCritical_ix + * - R + - `Recovered` + - Recovered_x, Recoverd_ab + * - D + - `Dead` + - Dead_x + +The compartments :math:`C, I` are infectious (red), the other compartments are considered to be not infectious (blue), +due to extensive isolation of the hospitalized individuals. + + +Infection States +---------------- + +The model contains the following list of **InfectionState**\s: + +.. code-block:: RST + + `Susceptible` + + `Exposed_1a` + `InfectedNoSymptoms_1a` + `InfectedSymptoms_1a` + `InfectedSevere_1a` + `InfectedCritical_1a` + + `Recovered_a` + `Dead_a` + + `Exposed_2a` + `InfectedNoSymptoms_2a` + `InfectedSymptoms_2a` + `InfectedSevere_2a` + `InfectedCritical_2a` + + `Exposed_1b` + `InfectedNoSymptoms_1b` + `InfectedSymptoms_1b` + `InfectedSevere_1b` + `InfectedCritical_1b` + + `Recovered_b` + `Dead_b` + + `Exposed_2b` + `InfectedNoSymptoms_2b` + `InfectedSymptoms_2b` + `InfectedSevere_2b` + `InfectedCritical_2b` + + `Recovered_ab` + +| It is possible to include subcompartments for the states `E`, `C`, `I`, `H`, and `U`, so compartments +| Exposed_1a, Exposed_2a, Exposed_1b, Exposed_2b, +| InfectedNoSymptoms_1a, InfectedNoSymptoms_2a, InfectedNoSymptoms_1b, InfectedNoSymptoms_2b, +| InfectedSymptoms_1a , InfectedSymptoms_2a, InfectedSymptoms_1b , InfectedSymptoms_2b, +| InfectedSevere_1a, InfectedSevere_2a , InfectedSevere_1b , InfectedSevere_2b, +| InfectedCritical_1a, InfectedCritical_2a, InfectedCritical_1b, and InfectedCritical_2b. + +The number of subcompartments can be set individually for each compartment. + + +Sociodemographic Stratification +------------------------------- + +In the LCT-SECIR-2D model, the population can be stratified by one sociodemographic dimension. +This dimension is denoted **Group**. It can be used for age groups as well as for other interpretations. +Different age groups can have different numbers of subcompartments. + +Parameters +---------- + +The model implements the following parameters: + +.. list-table:: + :header-rows: 1 + :widths: 20 20 60 + + * - Mathematical variable + - C++ variable name + - Description + * - :math:`\phi` + - ``ContactPatterns`` + - Average number of contacts for person per day, for multiple age groups this is a matrix. + * - :math:`\rho_a` + - ``TransmissionProbabilityOnContact_a`` + - Transmission risk for people located in the susceptible compartments for disease :math:`a`. + * - :math:`\rho_b` + - ``TransmissionProbabilityOnContact_b`` + - Transmission risk for people located in the susceptible compartments for disease :math:`b`. + * - :math:`\xi_{I_{NS}, a}` + - ``RelativeTransmissionNoSymptoms_a`` + - Proportion of nonsymptomatically infected people who are not isolated for disease :math:`a`. + * - :math:`\xi_{I_{NS}, b}` + - ``RelativeTransmissionNoSymptoms_b`` + - Proportion of nonsymptomatically infected people who are not isolated for disease :math:`b`. + * - :math:`\xi_{I_{Sy}, a}` + - ``RiskOfInfectionFromSymptomatic_a`` + - Proportion of infected people with symptoms who are not isolated for disease :math:`a`. + * - :math:`\xi_{I_{Sy}, b}` + - ``RiskOfInfectionFromSymptomatic_b`` + - Proportion of infected people with symptoms who are not isolated for disease :math:`b`. + * - :math:`N` + - ``m_N0`` + - Total population. + * - :math:`n_{E,1a}` + - Defined in ``LctStates`` + - Number of subcompartments of the Exposed_1a compartment. + * - :math:`n_{E,2a}` + - Defined in ``LctStates`` + - Number of subcompartments of the Exposed_2a compartment. + * - :math:`n_{E,1b}` + - Defined in ``LctStates`` + - Number of subcompartments of the Exposed_1b compartment. + * - :math:`n_{E,2b}` + - Defined in ``LctStates`` + - Number of subcompartments of the Exposed_2b compartment. + * - :math:`n_{NS,1a}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedNoSymptoms_1a compartment. + * - :math:`n_{NS,2a}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedNoSymptoms_2a compartment. + * - :math:`n_{NS,1b}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedNoSymptoms_1b compartment. + * - :math:`n_{NS,2b}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedNoSymptoms_2b compartment. + * - :math:`n_{Sy,1a}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedSymptoms_1a compartment. + * - :math:`n_{Sy,2a}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedSymptoms_2a compartment. + * - :math:`n_{Sy,1b}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedSymptoms_1b compartment. + * - :math:`n_{Sy,2b}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedSymptoms_2b compartment. + * - :math:`n_{Sev,1a}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedSevere_1a compartment. + * - :math:`n_{Sev,2a}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedSevere_2a compartment. + * - :math:`n_{Sev,1b}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedSevere_1b compartment. + * - :math:`n_{Sev,2b}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedSevere_2b compartment. + * - :math:`n_{Cr,1a}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedCritical_1a compartment. + * - :math:`n_{Cr,2a}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedCritical_2a compartment. + * - :math:`n_{Cr,1b}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedCritical_1b compartment. + * - :math:`n_{Cr,2b}` + - Defined in ``LctStates`` + - Number of subcompartments of the InfectedCritical_2b compartment. + * - :math:`T_{E,a}` + - ``TimeExposed_a`` + - Average time in days an individual stays in the Exposed_1a or Exposed_2a compartment. + * - :math:`T_{E,b}` + - ``TimeExposed_b`` + - Average time in days an individual stays in the Exposed_1b or Exposed_2b compartment. + * - :math:`T_{I_{NS},a}` + - ``TimeInfectedNoSymptoms_a`` + - Average time in days an individual stays in the InfectedNoSymptomsa_1a or InfectedNoSymptoms_2a compartment. + * - :math:`T_{I_{NS},b}` + - ``TimeInfectedNoSymptoms_b`` + - Average time in days an individual stays in the InfectedNoSymptomsa_1b or InfectedNoSymptoms_2b compartment. + * - :math:`T_{I_{Sy},a}` + - ``TimeInfectedSymptoms_a`` + - Average time in days an individual stays in the InfectedSymptoms_1a or InfectedSymptoms_2a compartment. + * - :math:`T_{I_{Sy},b}` + - ``TimeInfectedSymptoms_b`` + - Average time in days an individual stays in the InfectedSymptoms_1b or InfectedSymptoms_2b compartment. + * - :math:`T_{I_{Sev},a}` + - ``TimeInfectedSevere_a`` + - Average time in days an individual stays in the InfectedSevere_1a or InfectedSevere_2a compartment. + * - :math:`T_{I_{Sev},b}` + - ``TimeInfectedSevere_b`` + - Average time in days an individual stays in the InfectedSevere_1b or InfectedSevere_2b compartment. + * - :math:`T_{I_{Cr},a}` + - ``TimeInfectedCritical_a`` + - Average time in days an individual stays in the InfectedCritical_1a or InfectedCritical_2a compartment. + * - :math:`T_{I_{Cr},b}` + - ``TimeInfectedCritical_b`` + - Average time in days an individual stays in the InfectedCritical_1b or InfectedCritical_2b compartment. + * - :math:`\mu_{I_{NS},a}^{R}` + - ``RecoveredPerInfectedNoSymptoms_a`` + - Probability of transition from compartment InfectedNoSymptoms_1a to Recovered_a or from InfectedNoSymptoms_2a to Recovered_ab. + * - :math:`\mu_{I_{NS},b}^{R}` + - ``RecoveredPerInfectedNoSymptoms_b`` + - Probability of transition from compartment InfectedNoSymptoms_1b to Recovered_b or from InfectedNoSymptoms_2b to Recovered_ab. + * - :math:`\mu_{I_{Sy},a}^{I_{Sev}}` + - ``SeverePerInfectedSymptoms_a`` + - Probability of transition from compartment InfectedSymptoms_1a to InfectedSevere_1a or from InfectedSymptoms_2a to InfectedSevere_2a. + * - :math:`\mu_{I_{Sy},b}^{I_{Sev}}` + - ``SeverePerInfectedSymptoms_b`` + - Probability of transition from compartment InfectedSymptoms_1b to InfectedSevere_1b or from InfectedSymptoms_2b to InfectedSevere_2b. + * - :math:`\mu_{I_{Sev},a}^{I_{Cr}}` + - ``CriticalPerSevere_a`` + - Probability of transition from compartment InfectedSevere_1a to InfectedCritical_1a or from InfectedSevere_2a to InfectedCritical_2a. + * - :math:`\mu_{I_{Sev},b}^{I_{Cr}}` + - ``CriticalPerSevere_b`` + - Probability of transition from compartment InfectedSevere_1b to InfectedCritical_1b or from InfectedSevere_2b to InfectedCritical_2b. + * - :math:`\mu_{I_{Cr},a}^{D}` + - ``DeathsPerCritical_a`` + - Probability of dying when in compartment InfectedCritical_1a or InfectedCritical_2a. + * - :math:`\mu_{I_{Cr},b}^{D}` + - ``DeathsPerCritical_b`` + - Probability of dying when in compartment InfectedCritical_1b or InfectedCritical_2b. + + +Initial conditions +------------------ + +To initialize the model, we start by defining the number of subcompartments for every compartment and constructing the model with it. + +.. code-block:: cpp + + constexpr size_t NumExposed_1a = 1, NumInfectedNoSymptoms_1a = 1, NumInfectedSymptoms_1a = 1, + NumInfectedSevere_1a = 1, NumInfectedCritical_1a = 1, NumExposed_2a = 1, + NumInfectedNoSymptoms_2a = 1, NumInfectedSymptoms_2a = 1, NumInfectedSevere_2a = 1, + NumInfectedCritical_2a = 1, NumExposed_1b = 1, NumInfectedNoSymptoms_1b = 1, + NumInfectedSymptoms_1b = 1, NumInfectedSevere_1b = 1, NumInfectedCritical_1b = 1, + NumExposed_2b = 1, NumInfectedNoSymptoms_2b = 1, NumInfectedSymptoms_2b = 1, + NumInfectedSevere_2b = 1, 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; + Model model; + +For the simulation, we need initial values for all (sub)compartments. If we do not set the initial values manually, these are set to :math:`0` by default. + +We start with constructing a vector ``initial_populations`` that we will pass on to the model. It contains vectors for each compartment, +that contains a vector with initial values for the respective subcompartments. + +.. code-block:: cpp + + std::vector> 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}}; + + +We assert that the vector has the correct size by checking that the number of `InfectionState`\s and the numbers of subcompartments are correct. + +.. code-block:: cpp + + 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()) || + (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()) || + (initial_populations[(size_t)InfState::Dead_a].size() != + LctState::get_num_subcompartments()) || + (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())) { + mio::log_error("The length of at least one vector in initial_populations does not match the related number of " + "subcompartments."); + return 1; + } + +The initial populations in the model are set via: + +.. code-block:: cpp + + std::vector 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]; + } + +.. _Nonpharmaceutical Interventions: +Nonpharmaceutical Interventions +------------------------------- + +In the LCT-SECIR-2D model, nonpharmaceutical interventions (NPIs) are implemented through dampings in the contact matrix. +These dampings reduce the contact rates between different groups to simulate interventions. + +Basic dampings can be added to the contact matrix as follows: + +.. code-block:: cpp + + // Create a contact matrix with constant contact rate 10 (one age group). + mio::ContactMatrixGroup& contact_matrix = model.parameters.get(); + 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.)); + +For age-resolved models, you can apply different dampings to different groups: + +.. code-block:: cpp + + // Create a contact matrix with constant contact rate 10 between all age groups + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_agegroups, num_agegroups, 10)); + + // Add a damping that reduces contacts within the same age group by 70% starting at day 5. + contact_matrix.add_damping(Eigen::VectorX::Constant(num_agegroups, 0.7).asDiagonal(), + mio::SimulationTime(5.)); + +Simulation +---------- + +We can simulate the model from :math:`t_0` to :math:`t_{\max}` with initial step size :math:`dt` as follows: + +.. code-block:: cpp + + ScalarType t0 = 0; + ScalarType tmax = 10; + ScalarType dt = 0.5; + mio::TimeSeries result = mio::simulate(t0, tmax, dt, model); + + +Output +------ + +The simulation result is stratefied by subcompartments. The function ``calculate_compartments()`` aggregates the subcompartments by `InfectionState`\s. + +.. code-block:: cpp + + mio::TimeSeries population_no_subcompartments = model.calculate_compartments(result); + +You can access the data in the `mio::TimeSeries` object as follows: + +.. code-block:: cpp + + // Get the number of time points. + auto num_points = static_cast(result.get_num_time_points()); + + // Access data at a specific time point. + Eigen::VectorX value_at_time_i = result.get_value(i); + ScalarType time_i = result.get_time(i); + + // Access the last time point. + Eigen::VectorX last_value = result.get_last_value(); + ScalarType last_time = result.get_last_time(); + + +You can print the simulation results as a formatted table: + +.. code-block:: cpp + + // Print results to console with default formatting. + result.print_table(); + + // Print with custom column labels and width, and custom number of decimals. + 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); + +Additionally, you can export the results to a CSV file: + +.. code-block:: cpp + + // Export results to CSV with default settings. + result.export_csv("simulation_results.csv"); + + +Visualization +------------- + +To visualize the results of a simulation, you can use the Python package :doc:`m-plot <../../python/m-plot>` and its documentation. + + +Examples +-------- + +An example can be found at: + +- `examples/lct_secir_2_diseases.cpp `_ + + +Overview of the ``lsecir2d`` namespace: +----------------------------------------- + +.. doxygennamespace:: mio::lsecir2d