diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 5b192a4f5e..33be503ef9 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -146,6 +146,7 @@ if(MEMILIO_BUILD_MODELS) add_subdirectory(models/ode_secir) add_subdirectory(models/ode_secirvvs) add_subdirectory(models/lct_secir) + add_subdirectory(models/glct_secir) add_subdirectory(models/ide_secir) add_subdirectory(models/ide_seir) add_subdirectory(models/ode_seir) diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index b73975d79f..4cffc21467 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -115,6 +115,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(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}) + add_executable(history_example history.cpp) target_link_libraries(history_example PRIVATE memilio) target_compile_options(history_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) diff --git a/cpp/examples/glct_secir.cpp b/cpp/examples/glct_secir.cpp new file mode 100644 index 0000000000..ba1718f15c --- /dev/null +++ b/cpp/examples/glct_secir.cpp @@ -0,0 +1,211 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: 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 "glct_secir/model.h" +#include "glct_secir/parameters.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/math/eigen.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 GLCT SECIR model. + // This example recreates the lct_secir.cpp example using the GLCT model. + // This means, that we use the corresponding initial numbers, parameters and numbers of subcompartments that are + // equivalent to the choices in the LCT example. + // Parameters, initial values and the number of subcompartments are not meant to represent a realistic scenario. + // We need to double the choices of the number of subcompartments for some compartments, + // as we define different strains for different transition possibilities. For both strains, the same number of + // subcompartments is chosen. The transition probabilities are defined in the StartingProbabilities. + constexpr size_t NumExposed = 2, NumInfectedNoSymptoms = 6, NumInfectedSymptoms = 2, NumInfectedSevere = 2, + NumInfectedCritical = 10; + using Model = mio::glsecir::Model; + using LctState = Model::LctState; + using InfectionState = LctState::InfectionState; + + Model model; + + const ScalarType tmax = 10; + const ScalarType t0 = 0; + const ScalarType dt_init = 10; ///< Initially used step size for adaptive method. + + // Define some epidemiological parameters needed throughput the model definition and initialization. + const ScalarType timeExposed = 3.2; + const ScalarType timeInfectedNoSymptoms = 2.; + const ScalarType timeInfectedSymptoms = 5.8; + const ScalarType timeInfectedSevere = 9.5; + const ScalarType timeInfectedCritical = 7.1; + const ScalarType recoveredPerInfectedNoSymptoms = 0.09; + const ScalarType severePerInfectedSymptoms = 0.2; + const ScalarType criticalPerSevere = 0.25; + const ScalarType deathsPerCritical = 0.3; + + // 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. + // Split the initial population defined in the lct_secir.cpp example according to the transition probabilities in + // the two strains per compartment. + std::vector> initial_populations = { + {750}, + {30, 20}, + {20 * (1 - recoveredPerInfectedNoSymptoms), 10 * (1 - recoveredPerInfectedNoSymptoms), + 10 * (1 - recoveredPerInfectedNoSymptoms), 20 * recoveredPerInfectedNoSymptoms, + 10 * recoveredPerInfectedNoSymptoms, 10 * recoveredPerInfectedNoSymptoms}, + {50 * severePerInfectedSymptoms, 50 * (1 - severePerInfectedSymptoms)}, + {50 * criticalPerSevere, 50 * (1 - criticalPerSevere)}, + {10 * deathsPerCritical, 10 * deathsPerCritical, 5 * deathsPerCritical, 3 * deathsPerCritical, + 2 * deathsPerCritical, 10 * (1 - deathsPerCritical), 10 * (1 - deathsPerCritical), 5 * (1 - deathsPerCritical), + 3 * (1 - deathsPerCritical), 2 * (1 - deathsPerCritical)}, + {20}, + {10}}; + + // Assert that initial_populations has the right shape. + if (initial_populations.size() != (size_t)InfectionState::Count) { + mio::log_error("The number of vectors in initial_populations does not match the number of InfectionStates."); + return 1; + } + if ((initial_populations[(size_t)InfectionState::Susceptible].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(size_t)InfectionState::Exposed].size() != NumExposed) || + (initial_populations[(size_t)InfectionState::InfectedNoSymptoms].size() != NumInfectedNoSymptoms) || + (initial_populations[(size_t)InfectionState::InfectedSymptoms].size() != NumInfectedSymptoms) || + (initial_populations[(size_t)InfectionState::InfectedSevere].size() != NumInfectedSevere) || + (initial_populations[(size_t)InfectionState::InfectedCritical].size() != NumInfectedCritical) || + (initial_populations[(size_t)InfectionState::Recovered].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(size_t)InfectionState::Dead].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[mio::Index(i)] = flat_initial_populations[i]; + } + + // Set Parameters according to the LCT model definitions, e.g. use Erlang-distributions. + // Exposed. + // The get_default of the StartingProbabilities returns the first unit vector of the defined size. + // It is necessary to set it although the default method is used to define the length of the vector. + model.parameters.get() = + mio::glsecir::StartingProbabilitiesExposed().get_default( + LctState::get_num_subcompartments()); + // The get_default function returns the TransitionMatrix that is required to have an Erlang-distributed + // stay time with an average of timeExposed. + model.parameters.get() = + mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms().get_default( + LctState::get_num_subcompartments(), timeExposed); + // This definition of the StartingProbability and the TransitionMatrix lead to an Erlang-distributed latent stage. + + // InfectedNoSymptoms. + // For InfectedNoSymptoms, two strains has to be defined, one for the Transition + // InfectedNoSymptomsToInfectedSymptoms and one for InfectedNoSymptomsToRecovered. + // The strains have a length of NumInfectedNoSymptoms/2. each. + // The transition probability is included in the StartingProbability vector. + mio::Vector StartingProbabilitiesInfectedNoSymptoms = + mio::Vector::Zero(LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedNoSymptoms[0] = 1 - recoveredPerInfectedNoSymptoms; + StartingProbabilitiesInfectedNoSymptoms[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = recoveredPerInfectedNoSymptoms; + model.parameters.get() = + StartingProbabilitiesInfectedNoSymptoms; + // Define equal TransitionMatrices for the strains. + // They follow the same Erlang-distribution such that we get the same result as with one strain in the LCT model. + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedNoSymptomsToInfectedSymptoms().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), + timeInfectedNoSymptoms); + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), + timeInfectedNoSymptoms); + // Do the same for all compartments. + // InfectedSymptoms. + mio::Vector StartingProbabilitiesInfectedSymptoms = + mio::Vector::Zero(LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedSymptoms[0] = severePerInfectedSymptoms; + StartingProbabilitiesInfectedSymptoms[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 1 - severePerInfectedSymptoms; + model.parameters.get() = StartingProbabilitiesInfectedSymptoms; + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedSymptomsToInfectedSevere().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedSymptoms); + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedSymptoms); + // InfectedSevere. + mio::Vector StartingProbabilitiesInfectedSevere = + mio::Vector::Zero(LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedSevere[0] = criticalPerSevere; + StartingProbabilitiesInfectedSevere[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 1 - criticalPerSevere; + model.parameters.get() = StartingProbabilitiesInfectedSevere; + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedSevere); + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedSevereToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedSevere); + // InfectedCritical. + mio::Vector StartingProbabilitiesInfectedCritical = + mio::Vector::Zero(LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedCritical[0] = deathsPerCritical; + StartingProbabilitiesInfectedCritical[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 1 - deathsPerCritical; + model.parameters.get() = StartingProbabilitiesInfectedCritical; + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedCriticalToDead().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedCritical); + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedCriticalToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), timeInfectedCritical); + + model.parameters.get() = 0.05; + + 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.7; + model.parameters.get() = 0.25; + + // Perform a simulation. + mio::TimeSeries result = mio::simulate(t0, tmax, dt_init, model); + // The simulation result is divided by subcompartments as in the LCT model. + // 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_result = mio::interpolate_simulation_result(population_no_subcompartments, 0.1); + interpolated_result.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 10, 4); +} diff --git a/cpp/examples/lct_secir.cpp b/cpp/examples/lct_secir.cpp index 9a6308a5c3..6678633171 100644 --- a/cpp/examples/lct_secir.cpp +++ b/cpp/examples/lct_secir.cpp @@ -47,9 +47,9 @@ int main() // Variable defines whether the class Initializer is used to define an initial vector from flows or whether a manually // defined initial vector is used to initialize the LCT model. - bool use_initializer_flows = true; + bool use_initializer_flows = false; - ScalarType tmax = 20; + ScalarType tmax = 10; // Set Parameters. model.parameters.get()[0] = 3.2; @@ -155,10 +155,9 @@ int main() // Perform a simulation. mio::TimeSeries result = mio::simulate(0, tmax, 0.5, model); - //result.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8); // The simulation result is divided by subcompartments. - // We call the function calculate_comparttments to get a result according to the InfectionStates. + // 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", "E", "C", "I", "H", "U", "R", "D "}, 16, 8); + interpolated_results.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 12, 4); } diff --git a/cpp/memilio/epidemiology/lct_infection_state.h b/cpp/memilio/epidemiology/lct_infection_state.h index 457210f8e5..bf80a80079 100644 --- a/cpp/memilio/epidemiology/lct_infection_state.h +++ b/cpp/memilio/epidemiology/lct_infection_state.h @@ -20,6 +20,9 @@ #ifndef MIO_EPI_LCT_INFECTION_STATE_H #define MIO_EPI_LCT_INFECTION_STATE_H +#include "memilio/config.h" +#include "memilio/math/eigen.h" + #include namespace mio @@ -78,6 +81,53 @@ class LctInfectionState return index; } + /** + * @brief Cumulates a vector with the number of individuals in each subcompartment (with subcompartments + * according to the LctInfectionState) to produce a Vector that divides the population only into the infection + * states defined in InfectionStates. + * + * @param[in] subcompartments Vector with number of individuals in each subcompartment. + * The size of the vector has to match the LctInfectionState. + * @return Vector with accumulated values for the InfectionStates. + */ + static Vector calculate_compartments(const Vector& subcompartments) + { + assert(subcompartments.rows() == Count); + + Vector 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()]; + + return compartments; + } + static constexpr size_t Count{(... + Ns)}; private: diff --git a/cpp/models/glct_secir/CMakeLists.txt b/cpp/models/glct_secir/CMakeLists.txt new file mode 100644 index 0000000000..4534e383ef --- /dev/null +++ b/cpp/models/glct_secir/CMakeLists.txt @@ -0,0 +1,12 @@ +add_library(glct_secir + infection_state.h + model.h + model.cpp + parameters.h +) +target_link_libraries(glct_secir PUBLIC memilio) +target_include_directories(glct_secir PUBLIC + $ + $ +) +target_compile_options(glct_secir PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) \ No newline at end of file diff --git a/cpp/models/glct_secir/README.md b/cpp/models/glct_secir/README.md new file mode 100755 index 0000000000..a557d7623d --- /dev/null +++ b/cpp/models/glct_secir/README.md @@ -0,0 +1,74 @@ +# GLCT SECIR model + +This model is based on the Generalized Linear Chain Trick (GLCT). + +The GLCT provides the option to use phase-type distributed stay times in the compartments through the use of subcompartments. The Generalized Linear Chain Trick is an extension of the Linear Chain Trick (as the name already suggests). Phase-type distributions are dense in the field of all positive-valued distributions. Therefore, for any positive-valued distribution, a phase-type distribution of arbitrary precision can be identified. +The normal ODE models have (possibly unrealistic) exponentially distributed stay times. +The GLCT model can still be described by an ordinary differential equation system. + +For the concept see +- P. J. Hurtado and 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) +- P. J. Hurtado and C. Richards, "Building mean field ODE models using the generalized linear chain trick & Markov chain theory", 2021. (https://doi.org/10.1080/17513758.2021.1912418) + +Here, the eight compartments +- `Susceptible` ($S$), may become exposed at any time +- `Exposed` ($E$), becomes infected after some time +- `InfectedNoSymptoms` ($I_{NS}$), becomes InfectedSymptoms or Recovered after some time +- `InfectedSymptoms` ($I_{Sy}$), becomes InfectedSevere or Recovered after some time +- `InfectedSevere` ($I_{Sev}$), becomes InfectedCritical or Recovered after some time +- `InfectedCritical` ($I_{Cr}$), becomes Recovered or Dead after some time +- `Recovered` ($R$) +- `Dead` ($D$) + +are used to simulate the spread of the disease. +It is possible to include phase-type distributed stay times for the five compartments Exposed, InfectedNoSymptoms, InfectedSymptoms, InfectedSevere and InfectedCritical. + +## Model equations +Below is an overview of the model variables and the model equations are stated. For a simpler description let $\mathcal{Z}=\{{E,I_{NS},I_{Sy},I_{Sev},I_{Cr}}\}$ be the set of the compartments that can be divided into subcompartments. + +| Mathematical variable | C++ variable name | Description | +|---------------------------- | --------------- | -------------------------------------------------------------------------------------------------- | +| $\phi$ | `ContactPatterns` | Average number of contacts of a person per day. | +| $\rho$ | `TransmissionProbabilityOnContact` | Transmission risk for people located in the Susceptible compartments. | +| $\xi_{I_{NS}}$ | `RelativeTransmissionNoSymptoms` | Proportion of nonsymptomatically infected people who are not isolated. | +| $\xi_{I_{Sy}}$ | `RiskOfInfectionFromSymptomatic` | Proportion of infected people with symptoms who are not isolated. | +| $n_{z}$ | `Num(...)` | Number of subcompartments of compartment $z\in\mathcal{Z}$. (...) refers to the C++-name of $z$ as stated above.| +| $\boldsymbol{\alpha_{z}}$ | `StartingProbabilities(...)` | Vector of size $n_{z}$ with the initial probability of starting in any of the subcompartments of compartment $z\in\mathcal{Z}$. The entries should sum up to $1$. | +| $\mathbf{A_{z}^{*}}$ | `TransitionMatrix(...z)To(...*)` | Matrix describing the transitions in between of the subcompartments of $z\in\mathcal{Z}$ that describes the transition to the compartment $*$. | + +![equations](https://github.com/SciCompMod/memilio/assets/70579874/e1da5e1d-e719-4c16-9f14-45374be7c353) + +Note that the bold notation $\mathbf{z}(t)$ for $z\in\mathcal{Z}$ stands for a vector. If several transitions are possible from a compartment, the vector is split in order to be able to select the stay times until the transitions in each case phase-type distributed. +For example, the order $\mathbf{I_{\text{NS}}}(t)=[\mathbf{I_{\text{NS}}^{\text{Sy}}}(t),\mathbf{I_{\text{NS}}^{\text{R}}}(t)]^{T}$ is used. Similar holds true for the other compartments $z\in\mathcal{Z}$. + +Implicitly, the matrices $\mathbf{A_{z}^{*}}$ for one $z\in\mathcal{Z}$ are a block of a matrix $\mathbf{A_{z}}$ corresponding to the whole vector $\mathbf{z}(t)$. As we have no transitions in between the strains defined for different transition probabilities, we would have many zeros in the matrix. The matrix can be defined as + +```math +\mathbf{A_{z}}= +\begin{bmatrix} +\mathbf{A_{z}^{*_1}} & \mathbf{0} \\ +\mathbf{0} & \mathbf{A_{z}^{*_2}} +\end{bmatrix}, +``` + +where ${\*}\_{1}$ is the compartment of the first transition, e.g. $I_{\text{Sy}}$ for $z=I_{\text{NS}}$ and $*_2$ the compartment of the second possible transition, e.g. $R$. +Therefore, we just store the non-zero blocks of the matrix. +Using these parameters, the phase-type distribution that defines the stay time in compartment $z\in\mathcal{Z}$ has the probability density function + +$f(x)=\boldsymbol{\alpha_z}^T\hspace{3mu}e^{x\hspace{0.1em}\mathbf{A_z}}\hspace{3mu}(-\mathbf{A_z}\hspace{3mu}\boldsymbol{\Bbb{1}})$ for $x\in\mathbb{R}^{+}$ + +and the cumulative distribution function + +$F(x)=1-\boldsymbol{\alpha_z}^T\hspace{3mu}e^{x\hspace{0.1em}\mathbf{A_z}}\hspace{3mu}\boldsymbol{\Bbb{1}},$ + +where $e^{x\hspace{0.1em}\mathbf{A_z}}=\sum_{j=0}^{\infty}\frac{(x\hspace{2.5mu}\mathbf{A_z})^j}{j!}$ is the matrix exponential and $\boldsymbol{\Bbb{1}}$ is the vector containing ones of the matching size. Therefore, via changing the vector $\boldsymbol{\alpha_z}$ and the matrices $\mathbf{A_{z}^{*}}$, one can choose the stay time distribution appropriately. + +It is important that the sizes of the vectors and matrices match each other and satisfy some other conditions that are checked before a simulation. + +The compartment structure with subcompartments is the same as in the LCT-SECIR model. An overview of the model architecture can be found in the [README of the LCT model](../lct_secir/README.md). +For the GLCT model, some additional transitions are possible and we have more arrows in the model architecture. Below is an example for the Exposed compartment. Note that some indices are omitted (e.g. $n$ instead of $n_E$) to keep the picture simple. + +![tikzGLCTSECIR](https://github.com/user-attachments/assets/fc075b7a-6cd2-4e70-bdd0-a2f4b9f2cf53) +## Examples + +A simple example can be found at [GLCT minimal example](../../examples/glct_secir.cpp). diff --git a/cpp/models/glct_secir/infection_state.h b/cpp/models/glct_secir/infection_state.h new file mode 100644 index 0000000000..cec3a15d69 --- /dev/null +++ b/cpp/models/glct_secir/infection_state.h @@ -0,0 +1,46 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: 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 MIO_GLCT_SECIR_INFECTIONSTATE_H +#define MIO_GLCT_SECIR_INFECTIONSTATE_H + +namespace mio +{ +namespace glsecir +{ + +/// @brief The InfectionState enum describes the basic categories for the infection state of persons. +enum class InfectionState +{ + Susceptible = 0, + Exposed = 1, + InfectedNoSymptoms = 2, + InfectedSymptoms = 3, + InfectedSevere = 4, + InfectedCritical = 5, + Recovered = 6, + Dead = 7, + Count = 8 +}; + +} // namespace glsecir +} // namespace mio + +#endif // MIO_GLCT_SECIR_INFECTIONSTATE_H \ No newline at end of file diff --git a/cpp/models/glct_secir/model.cpp b/cpp/models/glct_secir/model.cpp new file mode 100644 index 0000000000..4e2d3f531d --- /dev/null +++ b/cpp/models/glct_secir/model.cpp @@ -0,0 +1,29 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: 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 "glct_secir/model.h" + +namespace mio +{ +namespace glsecir +{ + +} // namespace glsecir +} // namespace mio \ No newline at end of file diff --git a/cpp/models/glct_secir/model.h b/cpp/models/glct_secir/model.h new file mode 100644 index 0000000000..6501fe5d53 --- /dev/null +++ b/cpp/models/glct_secir/model.h @@ -0,0 +1,364 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: 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 MIO_GLCT_SECIR_MODEL_H +#define MIO_GLCT_SECIR_MODEL_H + +#include "glct_secir/parameters.h" +#include "glct_secir/infection_state.h" +#include "memilio/epidemiology/lct_infection_state.h" +#include "memilio/compartments/compartmentalmodel.h" +#include "memilio/epidemiology/populations.h" +#include "memilio/config.h" +#include "memilio/utils/logging.h" +#include "memilio/utils/time_series.h" +#include "memilio/math/eigen.h" + +namespace mio +{ +namespace glsecir +{ + +/** + * @brief Class that defines an GLCT-SECIR model. + * + * @tparam NumExposed The number of subcompartments used for the Exposed compartment. + * @tparam NumInfectedNoSymptoms The number of subcompartments used for the InfectedNoSymptoms compartment. + * @tparam NumInfectedSymptoms The number of subcompartments used for the InfectedSymptoms compartment. + * @tparam NumInfectedSevere The number of subcompartments used for the InfectedSevere compartment. + * @tparam NumInfectedCritical The number of subcompartments used for the InfectedCritical compartment. + */ +template +class Model + : public CompartmentalModel< + ScalarType, + LctInfectionState, + mio::Populations>, + Parameters> +{ + +public: + using LctState = + LctInfectionState; ///< This class specifies the number of subcompartments. + using Base = CompartmentalModel, Parameters>; + using typename Base::ParameterSet; + using typename Base::Populations; + + /// @brief Default constructor. + Model() + : Base(Populations({Index(LctState::Count)}, 0.), ParameterSet()) + { + } + + /** + * @brief Checks that the model satisfies all constraints (e.g. parameter or population constraints), and + * logs an error if constraints are not satisfied. + * + * @return Returns true if one or more constraints are not satisfied, false otherwise. + */ + bool check_constraints() const + { + auto params = this->parameters; + + // --- Check that the dimensions are consistent. --- + if ((Eigen::Index)LctState::template get_num_subcompartments() != + params.template get().rows()) { + log_error("Constraint check: Dimension of the parameters does not match the number of subcompartments for " + "the Exposed " + "compartment."); + return true; + } + if ((Eigen::Index)LctState::template get_num_subcompartments() != + params.template get().rows()) { + log_error( + "Constraint check: Dimension of the parameters does not match the number of subcompartments for the " + "InfectedNoSymptoms compartment."); + return true; + } + if ((Eigen::Index)LctState::template get_num_subcompartments() != + params.template get().rows()) { + log_error( + "Constraint check: Dimension of the parameters does not match the number of subcompartments for the " + "InfectedSymptoms compartment."); + return true; + } + if ((Eigen::Index)LctState::template get_num_subcompartments() != + params.template get().rows()) { + log_error("Constraint check: Dimension of the parameters does not match the number of subcompartments for " + "the InfectedSevere " + "compartment."); + return true; + } + if ((Eigen::Index)LctState::template get_num_subcompartments() != + params.template get().rows()) { + log_error( + "Constraint check: Dimension of the parameters does not match the number of subcompartments for the " + "InfectedCritical compartment."); + return true; + } + + return (params.check_constraints() || this->populations.check_constraints()); + } + + /** + * @brief Evaluates the right-hand-side f of the GLCT dydt = f(y, t). + * + * The GLCT-SECIR model is defined through ordinary differential equations of the form dydt = f(y, t). + * y is a vector containing 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 + { + dydt.setZero(); + + auto params = this->parameters; + auto total_population = pop.sum() - pop[LctState::template get_first_index()]; + + // Calculate sum of all subcompartments for InfectedNoSymptoms. + ScalarType InfectedNoSymptoms_sum = + pop.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) + .sum(); + // Calculate sum of all subcompartments for InfectedSymptoms. + ScalarType InfectedSymptoms_sum = + pop.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) + .sum(); + + // --- Susceptibles. --- + ScalarType season_val = 1 + params.template get() * + sin(3.141592653589793 * ((params.template get() + t) / 182.5 + 0.5)); + dydt[0] = -y[0] / total_population * season_val * params.template get() * + params.template get().get_cont_freq_mat().get_matrix_at(t)(0, 0) * + (params.template get() * InfectedNoSymptoms_sum + + params.template get() * InfectedSymptoms_sum); + + // --- Exposed. --- + dydt.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) -= + dydt[0] * params.template get(); + dydt.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) += + params.template get().transpose() * + y.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()); + + // --- InfectedNoSymptoms. --- + // Flow from Exposed To InfectedNoSymptoms. + dydt.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) = + -(params.template get() * + Vector::Ones(LctState::template get_num_subcompartments())) + .transpose() * + y.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) * + params.template get(); + // Flow from InfectedNoSymptoms To InfectedSymptoms. + size_t dimensionInfectedNoSymptomsToInfectedSymptoms = + params.template get().rows(); + dydt.segment(LctState::template get_first_index(), + dimensionInfectedNoSymptomsToInfectedSymptoms) += + params.template get().transpose() * + y.segment(LctState::template get_first_index(), + dimensionInfectedNoSymptomsToInfectedSymptoms); + // Flow from InfectedNoSymptoms To Recovered. + size_t dimensionInfectedNoSymptomsToRecovered = + params.template get().rows(); + dydt.segment(LctState::template get_first_index() + + dimensionInfectedNoSymptomsToInfectedSymptoms, + dimensionInfectedNoSymptomsToRecovered) += + params.template get().transpose() * + y.segment(LctState::template get_first_index() + + dimensionInfectedNoSymptomsToInfectedSymptoms, + dimensionInfectedNoSymptomsToRecovered); + // Add flow directly to Recovered compartment. + dydt[LctState::template get_first_index()] += + -(params.template get() * + Vector::Ones(dimensionInfectedNoSymptomsToRecovered)) + .transpose() * + y.segment(LctState::template get_first_index() + + dimensionInfectedNoSymptomsToInfectedSymptoms, + dimensionInfectedNoSymptomsToRecovered); + + // --- InfectedSymptoms. --- + // Flow from InfectedNoSymptoms To InfectedSymptoms. + dydt.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) = + -params.template get() * + (params.template get() * + Vector::Ones(dimensionInfectedNoSymptomsToInfectedSymptoms)) + .transpose() * + y.segment(LctState::template get_first_index(), + dimensionInfectedNoSymptomsToInfectedSymptoms); + // Flow from InfectedSymptoms To InfectedSevere. + size_t dimensionInfectedSymptomsToInfectedSevere = + params.template get().rows(); + dydt.segment(LctState::template get_first_index(), + dimensionInfectedSymptomsToInfectedSevere) += + params.template get().transpose() * + y.segment(LctState::template get_first_index(), + dimensionInfectedSymptomsToInfectedSevere); + // Flow from InfectedSymptoms To Recovered. + size_t dimensionInfectedSymptomsToRecovered = + params.template get().rows(); + dydt.segment(LctState::template get_first_index() + + dimensionInfectedSymptomsToInfectedSevere, + dimensionInfectedSymptomsToRecovered) += + params.template get().transpose() * + y.segment(LctState::template get_first_index() + + dimensionInfectedSymptomsToInfectedSevere, + dimensionInfectedSymptomsToRecovered); + // Add flow directly to Recovered compartment. + dydt[LctState::template get_first_index()] += + -(params.template get() * + Vector::Ones(dimensionInfectedSymptomsToRecovered)) + .transpose() * + y.segment(LctState::template get_first_index() + + dimensionInfectedSymptomsToInfectedSevere, + dimensionInfectedSymptomsToRecovered); + + // --- InfectedSevere. --- + // Flow from InfectedSymptoms To InfectedSevere. + dydt.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) = + -params.template get() * + (params.template get() * + Vector::Ones(dimensionInfectedSymptomsToInfectedSevere)) + .transpose() * + y.segment(LctState::template get_first_index(), + dimensionInfectedSymptomsToInfectedSevere); + // Flow from InfectedSevere To InfectedCritical. + size_t dimensionInfectedSevereToInfectedCritical = + params.template get().rows(); + dydt.segment(LctState::template get_first_index(), + dimensionInfectedSevereToInfectedCritical) += + params.template get().transpose() * + y.segment(LctState::template get_first_index(), + dimensionInfectedSevereToInfectedCritical); + // Flow from InfectedSevere To Recovered. + size_t dimensionInfectedSevereToRecovered = + params.template get().rows(); + dydt.segment(LctState::template get_first_index() + + dimensionInfectedSevereToInfectedCritical, + dimensionInfectedSevereToRecovered) += + params.template get().transpose() * + y.segment(LctState::template get_first_index() + + dimensionInfectedSevereToInfectedCritical, + dimensionInfectedSevereToRecovered); + // Add flow directly to Recovered compartment. + dydt[LctState::template get_first_index()] += + -(params.template get() * + Vector::Ones(dimensionInfectedSevereToRecovered)) + .transpose() * + y.segment(LctState::template get_first_index() + + dimensionInfectedSevereToInfectedCritical, + dimensionInfectedSevereToRecovered); + + // --- InfectedCritical. --- + // Flow from InfectedSevere To InfectedCritical. + dydt.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) = + -params.template get() * + (params.template get() * + Vector::Ones(dimensionInfectedSevereToInfectedCritical)) + .transpose() * + y.segment(LctState::template get_first_index(), + dimensionInfectedSevereToInfectedCritical); + // Flow from InfectedCritical To Dead. + size_t dimensionInfectedCriticalToDead = params.template get().rows(); + dydt.segment(LctState::template get_first_index(), + dimensionInfectedCriticalToDead) += + params.template get().transpose() * + y.segment(LctState::template get_first_index(), + dimensionInfectedCriticalToDead); + // Flow from InfectedCritical To Recovered. + size_t dimensionInfectedCriticalToRecovered = + params.template get().rows(); + dydt.segment(LctState::template get_first_index() + + dimensionInfectedCriticalToDead, + dimensionInfectedCriticalToRecovered) += + params.template get().transpose() * + y.segment(LctState::template get_first_index() + + dimensionInfectedCriticalToDead, + dimensionInfectedCriticalToRecovered); + // Add flow directly to Recovered compartment. + dydt[LctState::template get_first_index()] += + -(params.template get() * + Vector::Ones(dimensionInfectedCriticalToRecovered)) + .transpose() * + y.segment(LctState::template get_first_index() + + dimensionInfectedCriticalToDead, + dimensionInfectedCriticalToRecovered); + + // --- Dead. --- + dydt[LctState::template get_first_index()] = + -(params.template get() * + Vector::Ones(dimensionInfectedCriticalToDead)) + .transpose() * + y.segment(LctState::template get_first_index(), + dimensionInfectedCriticalToDead); + } + + /** + * @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 corresponding 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 + { + TimeSeries compartments_ts((Eigen::Index)InfectionState::Count); + if (!(LctState::Count == subcompartments_ts.get_num_elements())) { + log_error("Result does not match InfectionStates of the model."); + // Return a TimeSeries with values -1. + Vector error_output = Vector::Constant((Eigen::Index)InfectionState::Count, -1); + compartments_ts.add_time_point(-1, error_output); + return compartments_ts; + } + Vector compartments((Eigen::Index)InfectionState::Count); + for (Eigen::Index i = 0; i < subcompartments_ts.get_num_time_points(); ++i) { + compartments_ts.add_time_point(subcompartments_ts.get_time(i), + LctState::calculate_compartments(subcompartments_ts[i])); + } + return compartments_ts; + } +}; + +} // namespace glsecir +} // namespace mio + +#endif // MIO_GLCT_SECIR_MODEL_H diff --git a/cpp/models/glct_secir/parameters.h b/cpp/models/glct_secir/parameters.h new file mode 100644 index 0000000000..e87b72b7ea --- /dev/null +++ b/cpp/models/glct_secir/parameters.h @@ -0,0 +1,681 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: 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 MIO_GLCT_SECIR_PARAMS_H +#define MIO_GLCT_SECIR_PARAMS_H + +#include "memilio/config.h" +#include "memilio/utils/parameter_set.h" +#include "memilio/math/eigen.h" +#include "memilio/math/floating_point.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/utils/logging.h" + +namespace mio +{ +namespace glsecir +{ + +/*********************************************** +* Define Parameters of the GLCT-SECIHURD model * +***********************************************/ + +/// @brief Vector with the probability to start in any of the subcompartments of the Exposed compartment. +struct StartingProbabilitiesExposed { + using Type = Vector; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in the Exposed compartment. + * @param[in] numExposed Number of subcompartments of the Exposed compartment. + */ + static Type get_default(size_t numExposed) + { + Vector def = Vector::Zero(numExposed); + def[0] = 1.; + return def; + } + static std::string name() + { + return "StartingProbabilitiesExposed"; + } +}; + +/// @brief Transition matrix of the Exposed compartment. +struct TransitionMatrixExposedToInfectedNoSymptoms { + using Type = Eigen::MatrixXd; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in the Exposed compartment. + * @param[in] numExposed Number of subcompartments of the Exposed compartment. + * @param[in] timeExposed Average time spent in Exposed compartment in day unit. + */ + static Type get_default(size_t numExposed, ScalarType timeExposed = 1.) + { + Eigen::MatrixXd def = + Vector::Constant(numExposed, -(ScalarType)numExposed / timeExposed).asDiagonal(); + def.diagonal(1).setConstant((ScalarType)numExposed / timeExposed); + return def; + } + static std::string name() + { + return "TransitionMatrixExposedToInfectedNoSymptoms"; + } +}; + +/// @brief Vector with the probability to start in any of the subcompartments of the InfectedNoSymptoms compartment. +struct StartingProbabilitiesInfectedNoSymptoms { + using Type = Vector; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedNoSymptoms compartment. + * @param[in] numInfectedNoSymptoms Number of subcompartments of the InfectedNoSymptoms compartment. + */ + static Type get_default(size_t numInfectedNoSymptoms) + { + Vector def = Vector::Zero(numInfectedNoSymptoms); + def[0] = 1.; + return def; + } + static std::string name() + { + return "StartingProbabilitiesInfectedNoSymptoms"; + } +}; + +/** + * @brief Transition matrix of the phase-type distribution describing the stay time in the InfectedNoSymptoms + * compartment before developing symptoms. + */ +struct TransitionMatrixInfectedNoSymptomsToInfectedSymptoms { + using Type = Eigen::MatrixXd; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedNoSymptoms compartment + * before developing symptoms. + * @param[in] dimension Number of rows/columns of the transition matrix. + * @param[in] time Average time spent in InfectedNoSymptoms before developing symptoms in day unit. + */ + static Type get_default(size_t dimension, ScalarType time = 1.) + { + Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + def.diagonal(1).setConstant((ScalarType)dimension / time); + return def; + } + static std::string name() + { + return "TransitionMatrixInfectedNoSymptomsToInfectedSymptoms"; + } +}; + +/** + * @brief Transition matrix of the phase-type distribution describing the stay time in the InfectedNoSymptoms + * compartment before recovery. + */ +struct TransitionMatrixInfectedNoSymptomsToRecovered { + using Type = Eigen::MatrixXd; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedNoSymptoms compartment + * before recovery. + * @param[in] dimension Number of rows/columns of the transition matrix. + * @param[in] time Average time spent in InfectedNoSymptoms before recovery in day unit. + */ + static Type get_default(size_t dimension, ScalarType time = 1.) + { + Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + def.diagonal(1).setConstant((ScalarType)dimension / time); + return def; + } + static std::string name() + { + return "TransitionMatrixInfectedNoSymptomsToInfectedSymptomsToRecovered"; + } +}; + +/// @brief Vector with the probability to start in any of the subcompartments of the InfectedSymptoms compartment. +struct StartingProbabilitiesInfectedSymptoms { + using Type = Vector; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedSymptoms compartment. + * @param[in] numInfectedSymptoms Number of subcompartments of the InfectedSymptoms compartment. + */ + static Type get_default(size_t numInfectedSymptoms) + { + Vector def = Vector::Zero(numInfectedSymptoms); + def[0] = 1.; + return def; + } + static std::string name() + { + return "StartingProbabilitiesInfectedSymptoms"; + } +}; + +/** + * @brief Transition matrix of the phase-type distribution describing the stay time in the InfectedNoSymptoms + * compartment before going to hospital. + */ +struct TransitionMatrixInfectedSymptomsToInfectedSevere { + using Type = Eigen::MatrixXd; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in the InfectedSymptoms compartment + * before going to hospital. + * @param[in] dimension Number of rows/columns of the transition matrix. + * @param[in] time Average time spent in InfectedSymptoms before going to hospital in day unit. + */ + static Type get_default(size_t dimension, ScalarType time = 1.) + { + Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + def.diagonal(1).setConstant((ScalarType)dimension / time); + return def; + } + static std::string name() + { + return "TransitionMatrixInfectedSymptomsToInfectedSevere"; + } +}; + +/** + * @brief Transition matrix of the phase-type distribution describing the stay time in the InfectedSymptoms + * compartment before recovery. + */ +struct TransitionMatrixInfectedSymptomsToRecovered { + using Type = Eigen::MatrixXd; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in the InfectedSymptoms compartment + * before recovery. + * @param[in] dimension Number of rows/columns of the transition matrix. + * @param[in] time Average time spent in InfectedSymptoms before recovery in day unit. + */ + static Type get_default(size_t dimension, ScalarType time = 1.) + { + Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + def.diagonal(1).setConstant((ScalarType)dimension / time); + return def; + } + static std::string name() + { + return "TransitionMatrixInfectedSymptomsToRecovered"; + } +}; + +/// @brief Vector with the probability to start in any of the subcompartments of the InfectedSevere compartment. +struct StartingProbabilitiesInfectedSevere { + using Type = Vector; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedSevere compartment. + * @param[in] numInfectedSevere Number of subcompartments of the InfectedSevere compartment. + */ + static Type get_default(size_t numInfectedSevere) + { + Vector def = Vector::Zero(numInfectedSevere); + def[0] = 1.; + return def; + } + static std::string name() + { + return "StartingProbabilitiesInfectedSevere"; + } +}; + +/** + * @brief Transition matrix of the phase-type distribution describing the stay time in the InfectedSevere + * compartment before treated by ICU. + */ +struct TransitionMatrixInfectedSevereToInfectedCritical { + using Type = Eigen::MatrixXd; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedSevere compartment + * before treated by ICU. + * @param[in] dimension Number of rows/columns of the transition matrix. + * @param[in] time Average time spent in InfectedSevere before treated by ICU in day unit. + */ + static Type get_default(size_t dimension, ScalarType time = 1.) + { + Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + def.diagonal(1).setConstant((ScalarType)dimension / time); + return def; + } + static std::string name() + { + return "TransitionMatrixInfectedSevereToInfectedCritical"; + } +}; + +/** + * @brief Transition matrix of the phase-type distribution describing the stay time in the InfectedSevere + * compartment before recovery. + */ +struct TransitionMatrixInfectedSevereToRecovered { + using Type = Eigen::MatrixXd; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedSevere compartment + * before recovery. + * @param[in] dimension Number of rows/columns of the transition matrix. + * @param[in] time Average time spent in InfectedSevere before recovery in day unit. + */ + static Type get_default(size_t dimension, ScalarType time = 1.) + { + Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + def.diagonal(1).setConstant((ScalarType)dimension / time); + return def; + } + static std::string name() + { + return "TransitionMatrixInfectedSevereToRecovered"; + } +}; + +/// @brief Vector with the probability to start in any of the subcompartments of the InfectedCritical compartment. +struct StartingProbabilitiesInfectedCritical { + using Type = Vector; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedCritical compartment. + * @param[in] numInfectedCritical Number of subcompartments of the InfectedCritical compartment. + */ + static Type get_default(size_t numInfectedCritical) + { + Vector def = Vector::Zero(numInfectedCritical); + def[0] = 1.; + return def; + } + static std::string name() + { + return "StartingProbabilitiesInfectedCritical"; + } +}; + +/** + * @brief Transition matrix of the phase-type distribution describing the stay time in the InfectedCritical + * compartment before death. + */ +struct TransitionMatrixInfectedCriticalToDead { + using Type = Eigen::MatrixXd; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedCritical compartment + * before death. + * @param[in] dimension Number of rows/columns of the transition matrix. + * @param[in] time Average time treated by ICU before dying in day unit. + */ + static Type get_default(size_t dimension, ScalarType time = 1.) + { + Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + def.diagonal(1).setConstant((ScalarType)dimension / time); + return def; + } + static std::string name() + { + return "TransitionMatrixInfectedCriticalToDead"; + } +}; + +/** + * @brief Transition matrix of the phase-type distribution describing the stay time in the InfectedCritical + * compartment before recovery. + */ +struct TransitionMatrixInfectedCriticalToRecovered { + using Type = Eigen::MatrixXd; + /** + * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedCritical compartment + * before recovery. + * @param[in] dimension Number of rows/columns of the transition matrix. + * @param[in] time Average time treated by ICU before recovery in day unit. + */ + static Type get_default(size_t dimension, ScalarType time = 1.) + { + Eigen::MatrixXd def = Vector::Constant(dimension, -(ScalarType)dimension / time).asDiagonal(); + def.diagonal(1).setConstant((ScalarType)dimension / time); + return def; + } + static std::string name() + { + return "TransitionMatrixInfectedCriticalToRecovered"; + } +}; + +/// @brief Probability of getting infected from a contact. +struct TransmissionProbabilityOnContact { + using Type = ScalarType; + static Type get_default() + { + return 1.0; + } + static std::string name() + { + return "TransmissionProbabilityOnContact"; + } +}; + +/// @brief The contact patterns within the society are modelled using an UncertainContactMatrix. +struct ContactPatterns { + using Type = UncertainContactMatrix; + + static Type get_default() + { + ContactMatrixGroup contact_matrix = ContactMatrixGroup(1, 1); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10.)); + return Type(contact_matrix); + } + static std::string name() + { + return "ContactPatterns"; + } +}; + +/// @brief The relative InfectedNoSymptoms infectability. +struct RelativeTransmissionNoSymptoms { + using Type = ScalarType; + static Type get_default() + { + return 0.5; + } + static std::string name() + { + return "RelativeTransmissionNoSymptoms"; + } +}; + +/// @brief The risk of infection from symptomatic cases in the GLCT-SECIR model. +struct RiskOfInfectionFromSymptomatic { + using Type = ScalarType; + static Type get_default() + { + return 0.5; + } + static std::string name() + { + return "RiskOfInfectionFromSymptomatic"; + } +}; + +/** + * @brief The start day in the GLCT-SECIR 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() + { + return 0.; + } + static std::string name() + { + return "StartDay"; + } +}; + +/** + * @brief The seasonality in the GLCT-SECIR 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() + { + return Type(0.); + } + static std::string name() + { + return "Seasonality"; + } +}; + +using ParametersBase = + ParameterSet; + +/// @brief Parameters of an GLCT-SECIR model. +class Parameters : public ParametersBase +{ +public: + /// @brief Default constructor. + Parameters() + : ParametersBase() + { + } + + /** + * @brief Checks that all parameters satisfy their corresponding constraints and logs an error + * if constraints are not satisfied. + * + * @return Returns true if one or more constraints are not satisfied, false otherwise. + */ + bool check_constraints() const + { + // --- Parameters affecting the transmission of the virus. --- + if (this->get() < 0.0 || + this->get() > 1.0) { + log_error("Constraint check: Parameter TransmissionProbabilityOnContact smaller {:d} or larger {:d}", 0, 1); + return true; + } + + if (this->get() < 0.0 || this->get() > 1.0) { + log_error("Constraint check: Parameter RelativeTransmissionNoSymptoms smaller {:d} or larger {:d}", 0, 1); + return true; + } + + if (this->get() < 0.0 || this->get() > 1.0) { + log_error("Constraint check: Parameter RiskOfInfectionFromSymptomatic smaller {:d} or larger {:d}", 0, 1); + return true; + } + + 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; + } + + // --- Parameters affecting the phase-type distributions. --- + // --- Check that the dimensions are consistent. --- + if ((this->get().cols() != + this->get().rows()) || + (this->get().cols() != + this->get().rows()) || + (this->get().cols() != + this->get().rows()) || + (this->get().cols() != + this->get().rows()) || + (this->get().cols() != + this->get().rows()) || + (this->get().cols() != + this->get().rows()) || + (this->get().cols() != + this->get().rows()) || + (this->get().cols() != + this->get().rows()) || + (this->get().cols() != + this->get().rows())) { + log_error("Constraint check: At least one of the matrices used for the TransitionMatrix parameters is not " + "quadratic."); + return true; + } + + if (this->get().rows() != + this->get().rows()) { + log_error("Constraint check: Dimensions of StartingProbabilitiesExposed and " + "TransitionMatrixExposedToInfectedNoSymptoms are not matching."); + return true; + } + + if (this->get().rows() != + this->get().rows() + + this->get().rows()) { + log_error("Constraint check: Dimensions of StartingProbabilitiesInfectedNoSymptoms and " + "TransitionMatrices of InfectedNoSymptoms compartment are not matching."); + return true; + } + + if (this->get().rows() != + this->get().rows() + + this->get().rows()) { + log_error("Constraint check: Dimensions of StartingProbabilitiesInfectedSymptoms and " + "TransitionMatrices of InfectedSymptoms compartment are not matching."); + return true; + } + + if (this->get().rows() != + this->get().rows() + + this->get().rows()) { + log_error("Constraint check: Dimensions of StartingProbabilitiesInfectedSevere and " + "TransitionMatrices of InfectedSevere compartment are not matching."); + return true; + } + + if (this->get().rows() != + this->get().rows() + + this->get().rows()) { + log_error("Constraint check: Dimensions of StartingProbabilitiesInfectedCritical and " + "TransitionMatrices of InfectedCritical compartment are not matching."); + return true; + } + + // --- Check constraints of the starting probability vectors. --- + if ((!floating_point_equal(1., this->get().sum())) || + (!floating_point_equal(1., this->get().sum())) || + (!floating_point_equal(1., this->get().sum())) || + (!floating_point_equal(1., this->get().sum())) || + (!floating_point_equal(1., this->get().sum()))) { + log_warning( + "Constraint check: At least one of the vectors for the starting probabilities does not sum to one."); + return true; + } + + if ((this->get().array() < -1e-10).any() || + (this->get().array() < -1e-10).any() || + (this->get().array() < -1e-10).any() || + (this->get().array() < -1e-10).any() || + (this->get().array() < -1e-10).any()) { + log_warning("Constraint check: At least one of the vectors for the starting probabilities has at least one " + "negative entry."); + return true; + } + + // --- Check that we have no flows back from one compartment to the previous one + // (only in between of the subcompartments). --- + if (((this->get() * + Vector::Ones(this->get().rows())) + .array() > 1e-10) + .any()) { + log_warning( + "Constraint check: The entries of TransitionMatrixExposedToInfectedNoSymptoms lead to a negative " + "flow ExposedToInfectedNoSymptoms."); + return true; + } + if (((this->get() * + Vector::Ones(this->get().rows())) + .array() > 1e-10) + .any()) { + log_warning("Constraint check: The entries of TransitionMatrixInfectedNoSymptomsToInfectedSymptoms lead to " + "a negative " + "flow InfectedNoSymptomsToInfectedSymptoms."); + return true; + } + if (((this->get() * + Vector::Ones(this->get().rows())) + .array() > 1e-10) + .any()) { + log_warning( + "Constraint check: The entries of TransitionMatrixInfectedNoSymptomsToRecovered lead to a negative " + "flow InfectedNoSymptomsToRecovered."); + return true; + } + if (((this->get() * + Vector::Ones(this->get().rows())) + .array() > 1e-10) + .any()) { + log_warning( + "Constraint check: The entries of TransitionMatrixInfectedSymptomsToInfectedSevere lead to a negative " + "flow InfectedSymptomsToInfectedSevere."); + return true; + } + if (((this->get() * + Vector::Ones(this->get().rows())) + .array() > 1e-10) + .any()) { + log_warning( + "Constraint check: The entries of TransitionMatrixInfectedSymptomsToRecovered lead to a negative " + "flow InfectedSymptomsToRecovered."); + return true; + } + if (((this->get() * + Vector::Ones(this->get().rows())) + .array() > 1e-10) + .any()) { + log_warning( + "Constraint check: The entries of TransitionMatrixInfectedSevereToInfectedCritical lead to a negative " + "flow InfectedSevereToInfectedCritical."); + return true; + } + if (((this->get() * + Vector::Ones(this->get().rows())) + .array() > 1e-10) + .any()) { + log_warning("Constraint check: The entries of TransitionMatrixInfectedSevereToRecovered lead to a negative " + "flow InfectedSevereToRecovered."); + return true; + } + if (((this->get() * + Vector::Ones(this->get().rows())) + .array() > 1e-10) + .any()) { + log_warning("Constraint check: The entries of TransitionMatrixInfectedCriticalToDead lead to a negative " + "flow InfectedCriticalToDead."); + return true; + } + if (((this->get() * + Vector::Ones(this->get().rows())) + .array() > 1e-10) + .any()) { + log_warning( + "Constraint check: The entries of TransitionMatrixInfectedCriticalToRecovered lead to a negative " + "flow InfectedCriticalToRecovered."); + return true; + } + + return false; + } + +private: + Parameters(ParametersBase&& base) + : ParametersBase(std::move(base)) + { + } + +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 glsecir +} // namespace mio + +#endif // MIO_GLCT_SECIR_PARAMS_H diff --git a/cpp/models/lct_secir/README.md b/cpp/models/lct_secir/README.md index 4d6c990094..ee5b3c024f 100644 --- a/cpp/models/lct_secir/README.md +++ b/cpp/models/lct_secir/README.md @@ -33,14 +33,14 @@ Below is an overview of the model architecture and its compartments without a st | $\phi$ | `ContactPatterns` | Average number of contacts of a person per day. | | $\rho$ | `TransmissionProbabilityOnContact` | Transmission risk for people located in the susceptible compartments. | | $\xi_{I_{NS}}$ | `RelativeTransmissionNoSymptoms` | Proportion of nonsymptomatically infected people who are not isolated. | -| $\xi_{I_{Sy}}$ | `RiskOfInfectionFromSymptomatic` | Proportion of infected people with symptomps who are not isolated. | +| $\xi_{I_{Sy}}$ | `RiskOfInfectionFromSymptomatic` | Proportion of infected people with symptoms who are not isolated. | | $N$ | `m_N0` | Total population. | | $D$ | `D` | Number of death people. | -| $n_E$ | `NumExposed` | Number of subcompartments of the Exposed compartment. | -| $n_{NS}$ | `NumInfectedNoSymptoms` | Number of subcompartments of the InfectedNoSymptoms compartment. | -| $n_{Sy}$ | `NumInfectedSymptoms` | Number of subcompartments of the InfectedSymptoms compartment. | -| $n_{Sev}$ |`NumInfectedSevere` | Number of subcompartments of the InfectedSevere compartment.| -| $n_{Cr}$ | `NumInfectedCritical` | Number of subcompartments of the InfectedCritical compartment. | +| $n_E$ | Defined in `LctStates` | Number of subcompartments of the Exposed compartment. | +| $n_{NS}$ | Defined in `LctStates` | Number of subcompartments of the InfectedNoSymptoms compartment. | +| $n_{Sy}$ | Defined in `LctStates` | Number of subcompartments of the InfectedSymptoms compartment. | +| $n_{Sev}$ | Defined in `LctStates` | Number of subcompartments of the InfectedSevere compartment.| +| $n_{Cr}$ | Defined in `LctStates` | Number of subcompartments of the InfectedCritical compartment. | | $T_E$ | `TimeExposed` | Average time in days an individual stays in the Exposed compartment. | | $T_{I_{NS}}$ | `TimeInfectedNoSymptoms` | Average time in days an individual stays in the InfectedNoSymptoms compartment. | | $T_{I_{Sy}}$ | `TimeInfectedSymptoms` | Average time in days an individual stays in the InfectedSymptoms compartment. | diff --git a/cpp/models/lct_secir/infection_state.h b/cpp/models/lct_secir/infection_state.h index 0cea2df28b..25b24da169 100644 --- a/cpp/models/lct_secir/infection_state.h +++ b/cpp/models/lct_secir/infection_state.h @@ -18,8 +18,8 @@ * limitations under the License. */ -#ifndef LCTSECIR_INFECTIONSTATE_H -#define LCTSECIR_INFECTIONSTATE_H +#ifndef LCT_SECIR_INFECTIONSTATE_H +#define LCT_SECIR_INFECTIONSTATE_H namespace mio { @@ -65,4 +65,4 @@ enum class InfectionTransition } // namespace lsecir } // namespace mio -#endif \ No newline at end of file +#endif // LCT_SECIR_INFECTIONSTATE_H \ No newline at end of file diff --git a/cpp/models/lct_secir/model.h b/cpp/models/lct_secir/model.h index aef4287796..c20e9df290 100644 --- a/cpp/models/lct_secir/model.h +++ b/cpp/models/lct_secir/model.h @@ -153,53 +153,17 @@ class Model static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid."); using LctStateGroup = type_at_index_t; - // Specify first indices of the vector slices considered. - // First index of the group Group in an vector including all subcompartments. - Eigen::Index first_index_group_subcomps = this->populations.template get_first_index_of_group(); - Eigen::Index count_InfStates = (Eigen::Index)InfectionState::Count; - // First index of the group Group in an vector including all compartments without a resolution + // 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 segment of the vector subcompartments of each InfectionState and sum up the values of subcompartments. - compartments[first_index_group_comps + (Eigen::Index)InfectionState::Susceptible] = - subcompartments[first_index_group_subcomps]; - compartments[first_index_group_comps + (Eigen::Index)InfectionState::Exposed] = - subcompartments - .segment(first_index_group_subcomps + - LctStateGroup::template get_first_index(), - LctStateGroup::template get_num_subcompartments()) - .sum(); - compartments[first_index_group_comps + (Eigen::Index)InfectionState::InfectedNoSymptoms] = - subcompartments - .segment(first_index_group_subcomps + - LctStateGroup::template get_first_index(), - LctStateGroup::template get_num_subcompartments()) - .sum(); - compartments[first_index_group_comps + (Eigen::Index)InfectionState::InfectedSymptoms] = - subcompartments - .segment(first_index_group_subcomps + - LctStateGroup::template get_first_index(), - LctStateGroup::template get_num_subcompartments()) - .sum(); - compartments[first_index_group_comps + (Eigen::Index)InfectionState::InfectedSevere] = - subcompartments - .segment(first_index_group_subcomps + - LctStateGroup::template get_first_index(), - LctStateGroup::template get_num_subcompartments()) - .sum(); - compartments[first_index_group_comps + (Eigen::Index)InfectionState::InfectedCritical] = - subcompartments - .segment(first_index_group_subcomps + - LctStateGroup::template get_first_index(), - LctStateGroup::template get_num_subcompartments()) - .sum(); - compartments[first_index_group_comps + (Eigen::Index)InfectionState::Recovered] = - subcompartments[first_index_group_subcomps + - LctStateGroup::template get_first_index()]; - compartments[first_index_group_comps + (Eigen::Index)InfectionState::Dead] = - subcompartments[first_index_group_subcomps + - LctStateGroup::template get_first_index()]; + // 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); diff --git a/cpp/models/lct_secir/parameters.h b/cpp/models/lct_secir/parameters.h index 21a32f8f37..099b6c22ab 100644 --- a/cpp/models/lct_secir/parameters.h +++ b/cpp/models/lct_secir/parameters.h @@ -44,7 +44,7 @@ struct TimeExposed { using Type = Vector>; static Type get_default(size_t size) { - return Type::Constant(size, 1, 2.); + return Type::Constant(size, 1, 1.); } static std::string name() { @@ -76,7 +76,7 @@ struct TimeInfectedSymptoms { using Type = Vector>; static Type get_default(size_t size) { - return Type::Constant(size, 1, 1.5); + return Type::Constant(size, 1, 1.); } static std::string name() { diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 4a7be57cb2..b806543d93 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -66,6 +66,7 @@ set(TESTSOURCES test_state_age_function.cpp test_lct_secir.cpp test_lct_initializer_flows.cpp + test_glct_secir.cpp test_ad.cpp abm_helpers.h abm_helpers.cpp @@ -97,7 +98,7 @@ endif() add_executable(memilio-test ${TESTSOURCES}) target_include_directories(memilio-test PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) -target_link_libraries(memilio-test PRIVATE memilio ode_secir ode_seir ode_secirvvs ode_seair ide_seir ide_secir lct_secir abm gtest_main AD::AD) +target_link_libraries(memilio-test PRIVATE memilio ode_secir ode_seir ode_secirvvs ode_seair ide_seir ide_secir lct_secir glct_secir abm gtest_main AD::AD) target_compile_options(memilio-test PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) # make unit tests find the test data files diff --git a/cpp/tests/test_glct_secir.cpp b/cpp/tests/test_glct_secir.cpp new file mode 100644 index 0000000000..6b7b9c4a60 --- /dev/null +++ b/cpp/tests/test_glct_secir.cpp @@ -0,0 +1,527 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: 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 "glct_secir/model.h" +#include "glct_secir/parameters.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "memilio/utils/logging.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/compartments/simulation.h" +#include "memilio/math/eigen.h" +#include "load_test_data.h" + +#include +#include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp" + +// Test if the function eval_right_hand_side() is working using a hand calculated result. +TEST(TestGLCTSecir, testEvalRightHandSide) +{ + // Define initial values, parameters and numbers of subcompartments according to the choices of the + // testEvalRightHandSide of the LCT testing suite. For more details, + // we refer to the example glct_secir.cpp. + using Model = mio::glsecir::Model<2, 6, 4, 4, 4>; + using LctState = Model::LctState; + using InfectionState = LctState::InfectionState; + + Model model; + + // Set parameters such that the stay times are Erlang-distributed as in the corresponding LCT model. + // Exposed. + // Default functions are used to set the parameters but the corresponding dimensions have to be set manually. + model.parameters.get() = + mio::glsecir::StartingProbabilitiesExposed().get_default( + LctState::get_num_subcompartments()); + model.parameters.get() = + mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms().get_default( + LctState::get_num_subcompartments(), 3.2); + // InfectedNoSymptoms. + mio::Vector StartingProbabilitiesInfectedNoSymptoms = mio::Vector::Zero( + (Eigen::Index)LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedNoSymptoms[0] = 1 - 0.09; + StartingProbabilitiesInfectedNoSymptoms[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 0.09; + model.parameters.get() = + StartingProbabilitiesInfectedNoSymptoms; + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedNoSymptomsToInfectedSymptoms().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 2.); + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 2.); + // InfectedSymptoms. + mio::Vector StartingProbabilitiesInfectedSymptoms = mio::Vector::Zero( + (Eigen::Index)LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedSymptoms[0] = 0.2; + StartingProbabilitiesInfectedSymptoms[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 1 - 0.2; + model.parameters.get() = StartingProbabilitiesInfectedSymptoms; + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedSymptomsToInfectedSevere().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 5.8); + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 5.8); + // InfectedSevere. + mio::Vector StartingProbabilitiesInfectedSevere = mio::Vector::Zero( + (Eigen::Index)LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedSevere[0] = 0.25; + StartingProbabilitiesInfectedSevere[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 1 - 0.25; + model.parameters.get() = StartingProbabilitiesInfectedSevere; + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 9.5); + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedSevereToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 9.5); + // InfectedCritical. + mio::Vector StartingProbabilitiesInfectedCritical = mio::Vector::Zero( + (Eigen::Index)LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedCritical[0] = 0.3; + StartingProbabilitiesInfectedCritical[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 1 - 0.3; + model.parameters.get() = StartingProbabilitiesInfectedCritical; + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedCriticalToDead().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 7.1); + model.parameters.get() = + mio::glsecir::TransitionMatrixInfectedCriticalToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 7.1); + + model.parameters.get() = 0.05; + + mio::ContactMatrixGroup& contact_matrix = model.parameters.get(); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10)); + + model.parameters.get() = 0.7; + model.parameters.get() = 0.25; + model.parameters.get() = 0.; + model.parameters.get() = 0; + + // Define initial population distribution in infection states, one entry per subcompartment. + std::vector> initial_populations = { + {750}, + {30, 20}, + {20 * StartingProbabilitiesInfectedNoSymptoms[0], 10 * StartingProbabilitiesInfectedNoSymptoms[0], + 10 * StartingProbabilitiesInfectedNoSymptoms[0], 20 * (1 - StartingProbabilitiesInfectedNoSymptoms[0]), + 10 * (1 - StartingProbabilitiesInfectedNoSymptoms[0]), 10 * (1 - StartingProbabilitiesInfectedNoSymptoms[0])}, + {30 * StartingProbabilitiesInfectedSymptoms[0], 20 * StartingProbabilitiesInfectedSymptoms[0], + 30 * (1 - StartingProbabilitiesInfectedSymptoms[0]), 20 * (1 - StartingProbabilitiesInfectedSymptoms[0])}, + {40 * StartingProbabilitiesInfectedSevere[0], 10 * StartingProbabilitiesInfectedSevere[0], + 40 * (1 - StartingProbabilitiesInfectedSevere[0]), 10 * (1 - StartingProbabilitiesInfectedSevere[0])}, + {10 * StartingProbabilitiesInfectedCritical[0], 20 * StartingProbabilitiesInfectedCritical[0], + 10 * (1 - StartingProbabilitiesInfectedCritical[0]), 20 * (1 - StartingProbabilitiesInfectedCritical[0])}, + {20}, + {10}}; + std::vector flat_initial_populations; + for (auto&& vec : initial_populations) { + flat_initial_populations.insert(flat_initial_populations.end(), vec.begin(), vec.end()); + } + mio::Vector pop(LctState::Count); + for (size_t i = 0; i < LctState::Count; i++) { + pop[i] = flat_initial_populations[i]; + } + + // Compare the result of get_derivatives() with a hand calculated result. + mio::Vector dydt(LctState::Count); + model.get_derivatives(pop, pop, 0, dydt); + // This vector is the equivalent of the result defined in the test suite testEvalRightHandSide of the LCT model. + mio::Vector compare(LctState::Count); + compare << -15.3409, -3.4091, 6.25, -17.5 * 0.91, 15 * 0.91, 0 * 0.91, -17.5 * 0.09, 15 * 0.09, 0 * 0.09, + 3.3052 * 0.2, 3.4483 * 0.2, 3.3052 * 0.8, 3.4483 * 0.8, -7.0417 * 0.25, 6.3158 * 0.25, -7.0417 * 0.75, + 6.3158 * 0.75, -2.2906 * 0.3, -2.8169 * 0.3, -2.2906 * 0.7, -2.8169 * 0.7, 12.3899, 1.6901; + + for (size_t i = 0; i < LctState::Count; i++) { + EXPECT_NEAR(compare[i], dydt[i], 1e-3) << "Condition failed at index: " << i; + } +} + +// Model setup to compare result with a previous output of an LCT model. +class ModelTestGLCTSecir : public testing::Test +{ +public: + using Model = mio::glsecir::Model<2, 6, 2, 2, 10>; + using LctState = Model::LctState; + using InfectionState = LctState::InfectionState; + +protected: + virtual void SetUp() + { + model = new Model(); + // --- Set parameters. --- + // Exposed. + model->parameters.get() = + mio::glsecir::StartingProbabilitiesExposed().get_default( + LctState::get_num_subcompartments()); + model->parameters.get() = + mio::glsecir::TransitionMatrixExposedToInfectedNoSymptoms().get_default( + LctState::get_num_subcompartments(), 3.2); + // InfectedNoSymptoms. + mio::Vector StartingProbabilitiesInfectedNoSymptoms = mio::Vector::Zero( + (Eigen::Index)LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedNoSymptoms[0] = 1 - 0.09; + StartingProbabilitiesInfectedNoSymptoms[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 0.09; + model->parameters.get() = + StartingProbabilitiesInfectedNoSymptoms; + model->parameters.get() = + mio::glsecir::TransitionMatrixInfectedNoSymptomsToInfectedSymptoms().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 2.); + model->parameters.get() = + mio::glsecir::TransitionMatrixInfectedNoSymptomsToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 2.); + // InfectedSymptoms. + mio::Vector StartingProbabilitiesInfectedSymptoms = mio::Vector::Zero( + (Eigen::Index)LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedSymptoms[0] = 0.2; + StartingProbabilitiesInfectedSymptoms[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 1 - 0.2; + model->parameters.get() = + StartingProbabilitiesInfectedSymptoms; + model->parameters.get() = + mio::glsecir::TransitionMatrixInfectedSymptomsToInfectedSevere().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 5.8); + model->parameters.get() = + mio::glsecir::TransitionMatrixInfectedSymptomsToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 5.8); + // InfectedSevere. + mio::Vector StartingProbabilitiesInfectedSevere = mio::Vector::Zero( + (Eigen::Index)LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedSevere[0] = 0.25; + StartingProbabilitiesInfectedSevere[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 1 - 0.25; + model->parameters.get() = + StartingProbabilitiesInfectedSevere; + model->parameters.get() = + mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 9.5); + model->parameters.get() = + mio::glsecir::TransitionMatrixInfectedSevereToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 9.5); + // InfectedCritical. + mio::Vector StartingProbabilitiesInfectedCritical = mio::Vector::Zero( + (Eigen::Index)LctState::get_num_subcompartments()); + StartingProbabilitiesInfectedCritical[0] = 0.3; + StartingProbabilitiesInfectedCritical[(Eigen::Index)( + LctState::get_num_subcompartments() / 2.)] = 1 - 0.3; + model->parameters.get() = + StartingProbabilitiesInfectedCritical; + model->parameters.get() = + mio::glsecir::TransitionMatrixInfectedCriticalToDead().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 7.1); + model->parameters.get() = + mio::glsecir::TransitionMatrixInfectedCriticalToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 7.1); + + model->parameters.get() = 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.7; + model->parameters.get() = 0.25; + model->parameters.get() = 0.; + model->parameters.get() = 0; + + // --- Define initial distribution of the population in the subcompartments. --- + std::vector> initial_populations = { + {750}, + {30, 20}, + {20 * StartingProbabilitiesInfectedNoSymptoms[0], 10 * StartingProbabilitiesInfectedNoSymptoms[0], + 10 * StartingProbabilitiesInfectedNoSymptoms[0], 20 * (1 - StartingProbabilitiesInfectedNoSymptoms[0]), + 10 * (1 - StartingProbabilitiesInfectedNoSymptoms[0]), + 10 * (1 - StartingProbabilitiesInfectedNoSymptoms[0])}, + {50 * StartingProbabilitiesInfectedSymptoms[0], 50 * (1 - StartingProbabilitiesInfectedSymptoms[0])}, + {50 * StartingProbabilitiesInfectedSevere[0], 50 * (1 - StartingProbabilitiesInfectedSevere[0])}, + {10 * StartingProbabilitiesInfectedCritical[0], 10 * StartingProbabilitiesInfectedCritical[0], + 5 * StartingProbabilitiesInfectedCritical[0], 3 * StartingProbabilitiesInfectedCritical[0], + 2 * StartingProbabilitiesInfectedCritical[0], 10 * (1 - StartingProbabilitiesInfectedCritical[0]), + 10 * (1 - StartingProbabilitiesInfectedCritical[0]), 5 * (1 - StartingProbabilitiesInfectedCritical[0]), + 3 * (1 - StartingProbabilitiesInfectedCritical[0]), 2 * (1 - StartingProbabilitiesInfectedCritical[0])}, + {20}, + {10}}; + // 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[mio::Index(i)] = flat_initial_populations[i]; + } + } + + virtual void TearDown() + { + delete model; + } + +public: + Model* model = nullptr; +}; + +// Test compares a simulation with a previous output of an LCT model stored in a .csv file. +// This tests that the GLCT model reduces to an LCT model with the default parameters as well as nothing has changed +// in the implementation such that we still obtain an previous result. +TEST_F(ModelTestGLCTSecir, compareWithPreviousRun) +{ + ScalarType tmax = 3; + mio::TimeSeries result = mio::simulate( + 0, tmax, 0.5, *model, + std::make_shared>()); + + // Compare InfectionState compartments. + mio::TimeSeries population = model->calculate_compartments(result); + auto compare_population = load_test_data_csv("lct-secir-compartments-compare.csv"); + + ASSERT_EQ(compare_population.size(), static_cast(population.get_num_time_points())); + for (size_t i = 0; i < compare_population.size(); i++) { + ASSERT_EQ(compare_population[i].size(), static_cast(population.get_num_elements()) + 1) + << "at row " << i; + EXPECT_NEAR(population.get_time(i), compare_population[i][0], 1e-3) << "at row " << i; + for (size_t j = 1; j < compare_population[i].size(); j++) { + EXPECT_NEAR(population.get_value(i)[j - 1], compare_population[i][j], 1e-3) << " at row " << i; + } + } +} + +// Check constraints of Parameters class. +TEST_F(ModelTestGLCTSecir, testConstraintsModel) +{ + // Deactivate temporarily log output for next tests. + mio::set_log_level(mio::LogLevel::off); + + // --- Check correct setup. --- + bool constraint_check = model->check_constraints(); + EXPECT_FALSE(constraint_check); + + // Check if the number of subcompartments does not match the dimension of the vector with StartingProbabilities. + mio::Vector wrong_size = mio::Vector::Zero(3); + wrong_size[0] = 1; + // Exposed. + model->parameters.get() = wrong_size; + constraint_check = model->check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::StartingProbabilitiesExposed().get_default( + LctState::get_num_subcompartments()); + // InfectedNoSymptoms. + model->parameters.get() = wrong_size; + constraint_check = model->check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::StartingProbabilitiesInfectedNoSymptoms().get_default( + LctState::get_num_subcompartments()); + // InfectedSymptoms. + model->parameters.get() = wrong_size; + constraint_check = model->check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::StartingProbabilitiesInfectedSymptoms().get_default( + LctState::get_num_subcompartments()); + // InfectedSevere. + model->parameters.get() = wrong_size; + constraint_check = model->check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::StartingProbabilitiesInfectedSevere().get_default( + LctState::get_num_subcompartments()); + // InfectedCritical. + model->parameters.get() = wrong_size; + constraint_check = model->check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::StartingProbabilitiesInfectedCritical().get_default( + LctState::get_num_subcompartments()); +} + +// Check constraints of Parameters class. +TEST_F(ModelTestGLCTSecir, testConstraintsParameters) +{ + // Deactivate temporarily log output for next tests. + mio::set_log_level(mio::LogLevel::off); + + // --- Check correct setup. --- + bool constraint_check = model->parameters.check_constraints(); + EXPECT_FALSE(constraint_check); + + // --- Parameters affecting the transmission of the virus. --- + // Check TransmissionProbabilityOnContact. + model->parameters.get() = 5.1; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = 0.05; + + // Check RelativeTransmissionNoSymptoms. + model->parameters.get() = -0.05; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = 0.7; + + // Check RiskOfInfectionFromSymptomatic. + model->parameters.get() = 1.1; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = 0.25; + + // Check Seasonality. + model->parameters.get() = 0.6; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = 0.; + + // --- Check with incorrect dimensions. --- + // Check non-quadratic TransitionMatrixInfectedCriticalToRecovered. + model->parameters.get() = Eigen::MatrixXd::Zero(2, 3); + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::TransitionMatrixInfectedCriticalToRecovered().get_default( + (size_t)(LctState::get_num_subcompartments() / 2.), 7.1); + + // Check non matching dimensions of TransitionMatrix and vector with StartingProbabilities. + mio::Vector wrong_size = mio::Vector::Zero(3); + wrong_size[0] = 1; + // Exposed. + model->parameters.get() = wrong_size; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::StartingProbabilitiesExposed().get_default( + LctState::get_num_subcompartments()); + // InfectedNoSymptoms. + model->parameters.get() = wrong_size; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::StartingProbabilitiesInfectedNoSymptoms().get_default( + LctState::get_num_subcompartments()); + // InfectedSymptoms. + model->parameters.get() = wrong_size; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::StartingProbabilitiesInfectedSymptoms().get_default( + LctState::get_num_subcompartments()); + // InfectedSevere. + model->parameters.get() = wrong_size; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::StartingProbabilitiesInfectedSevere().get_default( + LctState::get_num_subcompartments()); + // InfectedCritical. + model->parameters.get() = wrong_size; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get() = + mio::glsecir::StartingProbabilitiesInfectedCritical().get_default( + LctState::get_num_subcompartments()); + + // --- Check constraints of the starting probability vectors. --- + model->parameters.get()[1] = 1.5; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()[0] = 1.1; + model->parameters.get()[1] = -0.1; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()[0] = 1.; + model->parameters.get()[1] = 0.; + + // --- Check with invalid transition matrices. --- + // ExposedToInfectedNoSymptoms. + model->parameters.get()(1, 0) = 10; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()(1, 0) = 0.1; + // InfectedNoSymptomsToInfectedSymptoms. + model->parameters.get()(2, 1) = 50; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()(2, 1) = 0.1; + // InfectedNoSymptomsToRecovered. + model->parameters.get()(1, 1) = -1.45; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()(1, 1) = -1.5; + // InfectedSymptomsToInfectedSevere. + model->parameters.get()(0, 0) = 1.; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()(0, 0) = -1.; + // InfectedSymptomsToRecovered. + model->parameters.get()(0, 0) = 0.1; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()(0, 0) = -0.1; + // InfectedSevereToInfectedCritical. + model->parameters.get()(0, 0) = 0.01; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()(0, 0) = -0.01; + // InfectedSevereToRecovered. + model->parameters.get()(0, 0) = 50; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()(0, 0) = -0.1; + // InfectedCriticalToDead. + model->parameters.get()(3, 1) = 6; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()(3, 1) = 0.; + // InfectedCriticalToRecovered. + model->parameters.get()(0, 4) = 3; + constraint_check = model->parameters.check_constraints(); + EXPECT_TRUE(constraint_check); + model->parameters.get()(0, 0) = -3.; + model->parameters.get()(0, 4) = 1.2; + + // --- Check with correct parameters. --- + constraint_check = model->parameters.check_constraints(); + EXPECT_FALSE(constraint_check); + + // Reactive log output. + mio::set_log_level(mio::LogLevel::warn); +} + +// Test calculate_compartments with a TimeSeries that has an incorrect number of elements. +TEST_F(ModelTestGLCTSecir, 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); + mio::Vector vec_wrong_size = mio::Vector::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); +}