From 43663a2ba2597c5785068527149385f10d31567d Mon Sep 17 00:00:00 2001 From: lenaploetzke Date: Fri, 2 Feb 2024 14:31:46 +0100 Subject: [PATCH 01/10] added lct model --- cpp/CMakeLists.txt | 1 + cpp/examples/CMakeLists.txt | 4 + cpp/examples/lct_secir.cpp | 91 +++ cpp/models/lct_secir/CMakeLists.txt | 14 + cpp/models/lct_secir/README.md | 23 + cpp/models/lct_secir/infection_state.h | 229 ++++++++ cpp/models/lct_secir/model.cpp | 195 +++++++ cpp/models/lct_secir/model.h | 102 ++++ cpp/models/lct_secir/parameters.h | 392 +++++++++++++ cpp/models/lct_secir/simulation.cpp | 58 ++ cpp/models/lct_secir/simulation.h | 177 ++++++ cpp/tests/CMakeLists.txt | 3 +- .../data/lct-secir-compartments-compare.csv | 14 + .../lct-secir-subcompartments-compare.csv | 14 + cpp/tests/test_lct_secir.cpp | 534 ++++++++++++++++++ 15 files changed, 1850 insertions(+), 1 deletion(-) create mode 100644 cpp/examples/lct_secir.cpp create mode 100644 cpp/models/lct_secir/CMakeLists.txt create mode 100644 cpp/models/lct_secir/README.md create mode 100644 cpp/models/lct_secir/infection_state.h create mode 100644 cpp/models/lct_secir/model.cpp create mode 100644 cpp/models/lct_secir/model.h create mode 100644 cpp/models/lct_secir/parameters.h create mode 100644 cpp/models/lct_secir/simulation.cpp create mode 100644 cpp/models/lct_secir/simulation.h create mode 100644 cpp/tests/data/lct-secir-compartments-compare.csv create mode 100644 cpp/tests/data/lct-secir-subcompartments-compare.csv create mode 100644 cpp/tests/test_lct_secir.cpp diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index a3ed0331df..de092bbc13 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -119,6 +119,7 @@ if(MEMILIO_BUILD_MODELS) add_subdirectory(models/abm) add_subdirectory(models/ode_secir) add_subdirectory(models/ode_secirvvs) + add_subdirectory(models/lct_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 3c5e4e7369..b9bf815c66 100644 --- a/cpp/examples/CMakeLists.txt +++ b/cpp/examples/CMakeLists.txt @@ -65,6 +65,10 @@ add_executable(ide_secir_example ide_secir.cpp) target_link_libraries(ide_secir_example PRIVATE memilio ide_secir) target_compile_options(ide_secir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) +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(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/lct_secir.cpp b/cpp/examples/lct_secir.cpp new file mode 100644 index 0000000000..b056d36fc3 --- /dev/null +++ b/cpp/examples/lct_secir.cpp @@ -0,0 +1,91 @@ +/* +* 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 "lct_secir/model.h" +#include "lct_secir/infection_state.h" +#include "lct_secir/simulation.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/math/eigen.h" +#include + +int main() +{ + /** Simple example to demonstrate how to simulate using an LCT SECIR model. + Parameters, initial values and subcompartments are not realistic. */ + + // Set vector that specifies the number of subcompartments. + std::vector num_subcompartments((int)mio::lsecir::InfectionStateBase::Count, 1); + num_subcompartments[(int)mio::lsecir::InfectionStateBase::Exposed] = 2; + num_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedNoSymptoms] = 3; + num_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedCritical] = 5; + mio::lsecir::InfectionState infectionState(num_subcompartments); + + ScalarType tmax = 20; + + // Define initial distribution of the population in the subcompartments. + Eigen::VectorXd init(infectionState.get_count()); + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::Susceptible)] = 750; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::Exposed)] = 30; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::Exposed) + 1] = 20; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms)] = 20; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms) + 1] = 10; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms) + 2] = 10; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSymptoms)] = 50; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSevere)] = 50; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical)] = 10; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 1] = 10; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 2] = 5; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 3] = 3; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 4] = 2; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::Recovered)] = 20; + init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::Dead)] = 10; + + // Initialize model. + mio::lsecir::Model model(std::move(init), infectionState); + + // Set Parameters. + model.parameters.get() = 2 * 4.2 - 5.2; + model.parameters.get() = 2 * (5.2 - 4.2); + model.parameters.get() = 5.8; + model.parameters.get() = 9.5; + // Also possible to change values with setter. + model.parameters.set(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(5.)); + + model.parameters.get() = 0.7; + model.parameters.get() = 0.25; + model.parameters.get() = 0.09; + model.parameters.get() = 0.2; + model.parameters.get() = 0.25; + model.parameters.set(0.3); + + // Perform a simulation. + mio::TimeSeries result = mio::lsecir::simulate(0, tmax, 0.5, model); + // Calculate the distribution in infectionState without subcompartments of the result and print it. + mio::TimeSeries populations = model.calculate_populations(result); + populations.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8); +} \ No newline at end of file diff --git a/cpp/models/lct_secir/CMakeLists.txt b/cpp/models/lct_secir/CMakeLists.txt new file mode 100644 index 0000000000..e93fbe8bed --- /dev/null +++ b/cpp/models/lct_secir/CMakeLists.txt @@ -0,0 +1,14 @@ +add_library(lct_secir + infection_state.h + model.h + model.cpp + simulation.h + simulation.cpp + parameters.h +) +target_link_libraries(lct_secir PUBLIC memilio) +target_include_directories(lct_secir PUBLIC + $ + $ +) +target_compile_options(lct_secir PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS}) \ No newline at end of file diff --git a/cpp/models/lct_secir/README.md b/cpp/models/lct_secir/README.md new file mode 100644 index 0000000000..462fad5626 --- /dev/null +++ b/cpp/models/lct_secir/README.md @@ -0,0 +1,23 @@ +# LCT SECIR model + +This model is based on the Linear Chain Trick. +For the concept see +- Lena Plötzke, "Der Linear Chain Trick in der epidemiologischen Modellierung als Kompromiss zwischen gewöhnlichen und Integro-Differentialgleichungen", 2023. (https://elib.dlr.de/200381/, german) +- P. J. Hurtado und A. S. Kirosingh, "Generalizations of the ‘Linear Chain Trick’: incorporating more flexible dwell time distributions into mean field ODE models“, 2019. (https://doi.org/10.1007/s00285-019-01412-w) + +The eight compartments +- Susceptible, may become exposed at any time +- Exposed, becomes infected after some time +- InfectedNoSymptoms, becomes InfectedSymptoms or Recovered after some time +- InfectedSymptoms, becomes InfectedSevere or Recovered after some time +- InfectedSevere, becomes InfectedCritical or Recovered after some time +- InfectedCritical, becomes Recovered or Dead after some time +- Recovered +- Dead + +are used to simulate the spread of the disease. +It is possible to include subcompartments for the five compartments Exposed, InfectedNoSymptoms, InfectedSymptoms, InfectedSevere and InfectedCritical. + +A simple example can be found at: examples/lct_secir.cpp. + + diff --git a/cpp/models/lct_secir/infection_state.h b/cpp/models/lct_secir/infection_state.h new file mode 100644 index 0000000000..ffca8d9871 --- /dev/null +++ b/cpp/models/lct_secir/infection_state.h @@ -0,0 +1,229 @@ +/* +* 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 LCTSECIR_INFECTIONSTATE_H +#define LCTSECIR_INFECTIONSTATE_H + +#include "memilio/utils/logging.h" +#include + +namespace mio +{ +namespace lsecir +{ + +/** + * @brief The InfectionStateBase enum describes the possible basic + * categories for the infection state of persons + */ +enum class InfectionStateBase +{ + Susceptible = 0, + Exposed = 1, + InfectedNoSymptoms = 2, + InfectedSymptoms = 3, + InfectedSevere = 4, + InfectedCritical = 5, + Recovered = 6, + Dead = 7, + Count = 8 +}; + +/** + * @brief The InfectionTransition enum describes the possible + * transitions of the infectious state of persons. + */ +enum class InfectionTransition +{ + SusceptibleToExposed = 0, + ExposedToInfectedNoSymptoms = 1, + InfectedNoSymptomsToInfectedSymptoms = 2, + InfectedNoSymptomsToRecovered = 3, + InfectedSymptomsToInfectedSevere = 4, + InfectedSymptomsToRecovered = 5, + InfectedSevereToInfectedCritical = 6, + InfectedSevereToRecovered = 7, + InfectedCriticalToDead = 8, + InfectedCriticalToRecovered = 9, + Count = 10 +}; + +class InfectionState +{ +public: + /** + * @brief Constructor for the InfectionState class. + * + * InfectionState class defines the possible InfectionState%s with the number of Subcompartments for the LCT model. + * With the default constructor, the class is defined without subcompartments, i.e. only the subdivision in InfectionStateBase + * is used. + */ + InfectionState() + : m_SubcompartmentNumbers(std::vector((int)InfectionStateBase::Count, 1)) + , m_SubcompartmentNumbersindexfirst(std::vector((int)InfectionStateBase::Count, 1)) + { + set_compartment_index(); + } + + /** + * @brief Constructor for the InfectionState class. + * + * InfectionState class defines the possible InfectionState%s with the number of Subcompartments for the LCT model. + * @param[in] SubcompartmentNumbers Vector which defines the number of Subcompartments for each infection state of InfectionStateBase. + */ + InfectionState(std::vector SubcompartmentNumbers) + : m_SubcompartmentNumbers(std::move(SubcompartmentNumbers)) + , m_SubcompartmentNumbersindexfirst(std::vector((int)InfectionStateBase::Count, 1)) + { + bool constraint_check = check_constraints(); + if (!constraint_check) { + set_compartment_index(); + } + } + + /** + * @brief Setter for the number of Subcompartments. + * + * @param[in] SubcompartmentNumbers Vector which defines the number of Subcompartments for each infection state of InfectionStateBase. + */ + void set_SubcompartmentNumbers(std::vector SubcompartmentNumbers) + { + m_SubcompartmentNumbers = SubcompartmentNumbers; + bool constraint_check = check_constraints(); + if (!constraint_check) { + set_compartment_index(); + } + } + + /** + * @brief Gets the number of subcompartments in an infection state. + * + * @param[in] infectionstatebase Infection state for which the number of subcompartments should be returned. + * @return Number of Subcompartments for infectionstatebase. + */ + int get_number(InfectionStateBase infectionstatebase) const + { + return m_SubcompartmentNumbers[(int)infectionstatebase]; + } + + /** + * @brief Gets the number of subcompartments in an infection state. + * + * @param[in] infectionstatebase Index of an infection state for which the number of subcompartments should be returned. + * @return Number of Subcompartments for infectionstatebase. + */ + int get_number(int infectionstatebaseindex) const + { + return m_SubcompartmentNumbers[infectionstatebaseindex]; + } + + /** + * @brief Gets the index of the first subcompartment in an vector with all subcompartments for an infection state. + * + * In a simulation, the number of individuals in the subcompartments are stored in vectors. + * Accordingly, the index in such a vector of the first subcompartment of an infection state is given. + * @param[in] infectionstatebase Infection state for which the index should be returned. + * @return Index of the first Subcompartment for a vector with one entry per subcompartment. + */ + int get_firstindex(InfectionStateBase infectionstatebase) const + { + return m_SubcompartmentNumbersindexfirst[(int)infectionstatebase]; + } + + /** + * @brief Gets the index of the first subcompartment in an vector with all subcompartments for an infection state. + * + * In a simulation, the number of individuals in the subcompartments are stored in vectors. + * Accordingly, the index in such a vector of the first subcompartment of an infection state is given. + * @param[in] infectionstatebase Index of an infection state for which the index of a vector should be returned. + * @return Index of the first Subcompartment for a vector with one entry per subcompartment. + */ + int get_firstindex(int infectionstatebaseindex) const + { + return m_SubcompartmentNumbersindexfirst[infectionstatebaseindex]; + } + + /** + * @brief Gets the total number of (sub-)compartments of infection states. + */ + int get_count() const + { + return m_Count; + } + + /** + * @brief Checks constraints on infection states. + * + * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false. + */ + bool check_constraints() const + { + if (!(m_SubcompartmentNumbers.size() == (int)InfectionStateBase::Count)) { + log_error("Vector for number of subcompartments has the wrong size."); + return true; + } + if (!(m_SubcompartmentNumbers[(int)InfectionStateBase::Susceptible] == 1)) { + log_error("Susceptible compartment can not have Subcompartments."); + return true; + } + if (!(m_SubcompartmentNumbers[(int)InfectionStateBase::Recovered] == 1)) { + log_error("Recovered compartment can not have Subcompartments."); + return true; + } + if (!(m_SubcompartmentNumbers[(int)InfectionStateBase::Dead] == 1)) { + log_error("Dead compartment can not have Subcompartments."); + return true; + } + for (int i = 0; i < (int)InfectionStateBase::Count; ++i) { + if (m_SubcompartmentNumbers[i] < 1) { + log_error("All compartments should have at least one Subcompartment."); + return true; + } + } + return false; + } + +private: + /** + * @brief Calculates Index of the first Subcompartment for a vector with one entry per subcompartment. + * + * Therefore the vector with number of subcompartments per infection state is used. + */ + void set_compartment_index() + { + int index = 0; + for (int i = 0; i < (int)(InfectionStateBase::Count); i++) { + m_SubcompartmentNumbersindexfirst[i] = index; + index = index + m_SubcompartmentNumbers[i]; + } + m_Count = index; + } + + std::vector + m_SubcompartmentNumbers; ///< Vector which defines the number of Subcompartments for each infection state of InfectionStateBase. + std::vector + m_SubcompartmentNumbersindexfirst; ///< Vector with Indexes for all infection states of the first Subcompartment for a vector with one entry per subcompartment. + int m_Count; ///< Total number of (sub-)compartments of infection states. +}; + +} // namespace lsecir +} // namespace mio + +#endif \ No newline at end of file diff --git a/cpp/models/lct_secir/model.cpp b/cpp/models/lct_secir/model.cpp new file mode 100644 index 0000000000..ca963ee800 --- /dev/null +++ b/cpp/models/lct_secir/model.cpp @@ -0,0 +1,195 @@ +/* +* 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 "lct_secir/model.h" +#include "infection_state.h" +#include "memilio/config.h" +#include "memilio/math/eigen.h" +#include "memilio/utils/logging.h" +#include "memilio/math/eigen.h" + +namespace mio +{ +namespace lsecir +{ + +Model::Model(Eigen::VectorXd init, const InfectionState infectionState_init, Parameters&& parameters_init) + : parameters{parameters_init} + , infectionState{infectionState_init} + , m_initial_values{std::move(init)} +{ + m_N0 = m_initial_values.sum(); +} + +bool Model::check_constraints() const +{ + if (!(infectionState.get_count() == m_initial_values.size())) { + log_error("Size of the initial values does not match Subcompartments."); + return true; + } + for (int i = 0; i < infectionState.get_count(); i++) { + if (m_initial_values[i] < 0) { + log_warning( + "Initial values for one subcompartment are less than zero. Simulation results are not realistic."); + return true; + } + } + return parameters.check_constraints(); +} + +void Model::eval_right_hand_side(Eigen::Ref y, ScalarType t, + Eigen::Ref dydt) const +{ + dydt.setZero(); + + ScalarType C = 0; + ScalarType I = 0; + ScalarType dummy = 0; + + // Calculate sum of all subcompartments for InfectedNoSymptoms. + C = y.segment(infectionState.get_firstindex(InfectionStateBase::InfectedNoSymptoms), + infectionState.get_number(InfectionStateBase::InfectedNoSymptoms)) + .sum(); + // Calculate sum of all subcompartments for InfectedSymptoms. + I = y.segment(infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms), + infectionState.get_number(InfectionStateBase::InfectedSymptoms)) + .sum(); + + // S' + ScalarType season_val = + 1 + parameters.get() * + sin(3.141592653589793 * (std::fmod((parameters.get() + t), 365.0) / 182.5 + 0.5)); + dydt[0] = + -y[0] / (m_N0 - y[infectionState.get_firstindex(InfectionStateBase::Dead)]) * season_val * + parameters.get() * + parameters.get().get_cont_freq_mat().get_matrix_at(t)(0, 0) * + (parameters.get() * C + parameters.get() * I); + + // E' + dydt[1] = -dydt[0]; + for (int i = 0; i < infectionState.get_number(InfectionStateBase::Exposed); i++) { + // Dummy stores the value of the flow from dydt[1 + i] to dydt[2 + i]. + // 1+i is always the index of a (sub-)compartment of E and 2+i can also be the index of the first (sub-)compartment of C. + dummy = infectionState.get_number(InfectionStateBase::Exposed) * (1 / parameters.get()) * y[1 + i]; + // Subtract flow from dydt[1 + i] and add to dydt[2 + i]. + dydt[1 + i] = dydt[1 + i] - dummy; + dydt[2 + i] = dummy; + } + + // C' + for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedNoSymptoms); i++) { + dummy = infectionState.get_number(InfectionStateBase::InfectedNoSymptoms) * + (1 / parameters.get()) * + y[infectionState.get_firstindex(InfectionStateBase::InfectedNoSymptoms) + i]; + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedNoSymptoms) + i] = + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedNoSymptoms) + i] - dummy; + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedNoSymptoms) + i + 1] = dummy; + } + + // I' + // Flow from last (sub-) compartment of C must be split between I_1 and R. + dydt[infectionState.get_firstindex(InfectionStateBase::Recovered)] = + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms)] * + parameters.get(); + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms)] = + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms)] * + (1 - parameters.get()); + + for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedSymptoms); i++) { + dummy = infectionState.get_number(InfectionStateBase::InfectedSymptoms) * + (1 / parameters.get()) * + y[infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms) + i]; + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms) + i] = + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms) + i] - dummy; + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms) + i + 1] = dummy; + } + + // H' + dydt[infectionState.get_firstindex(InfectionStateBase::Recovered)] = + dydt[infectionState.get_firstindex(InfectionStateBase::Recovered)] + + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSevere)] * + (1 - parameters.get()); + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSevere)] = + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSevere)] * + parameters.get(); + for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedSevere); i++) { + dummy = infectionState.get_number(InfectionStateBase::InfectedSevere) * + (1 / parameters.get()) * + y[infectionState.get_firstindex(InfectionStateBase::InfectedSevere) + i]; + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSevere) + i] = + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSevere) + i] - dummy; + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedSevere) + i + 1] = dummy; + } + + // U' + dydt[infectionState.get_firstindex(InfectionStateBase::Recovered)] = + dydt[infectionState.get_firstindex(InfectionStateBase::Recovered)] + + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedCritical)] * + (1 - parameters.get()); + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedCritical)] = + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedCritical)] * parameters.get(); + for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedCritical) - 1; i++) { + dummy = infectionState.get_number(InfectionStateBase::InfectedCritical) * + (1 / parameters.get()) * + y[infectionState.get_firstindex(InfectionStateBase::InfectedCritical) + i]; + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedCritical) + i] = + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedCritical) + i] - dummy; + dydt[infectionState.get_firstindex(InfectionStateBase::InfectedCritical) + i + 1] = dummy; + } + // Last flow from U has to be divided between R and D. + // Must be calculated separately in order not to overwrite the already calculated values ​​for R. + dummy = infectionState.get_number(InfectionStateBase::InfectedCritical) * + (1 / parameters.get()) * + y[infectionState.get_firstindex(InfectionStateBase::Recovered) - 1]; + dydt[infectionState.get_firstindex(InfectionStateBase::Recovered) - 1] = + dydt[infectionState.get_firstindex(InfectionStateBase::Recovered) - 1] - dummy; + dydt[infectionState.get_firstindex(InfectionStateBase::Recovered)] = + dydt[infectionState.get_firstindex(InfectionStateBase::Recovered)] + + (1 - parameters.get()) * dummy; + dydt[infectionState.get_firstindex(InfectionStateBase::Dead)] = parameters.get() * dummy; +} + +TimeSeries Model::calculate_populations(const TimeSeries& result) const +{ + if (!(infectionState.get_count() == result.get_num_elements())) { + log_error("Result does not match infectionState of the Model."); + TimeSeries populations((int)InfectionStateBase::Count); + Eigen::VectorXd wrong_size = Eigen::VectorXd::Constant((int)InfectionStateBase::Count, -1); + populations.add_time_point(-1, wrong_size); + return populations; + } + TimeSeries populations((int)InfectionStateBase::Count); + Eigen::VectorXd dummy((int)InfectionStateBase::Count); + for (Eigen::Index i = 0; i < result.get_num_time_points(); ++i) { + for (int j = 0; j < dummy.size(); ++j) { + // Use segment of vector of the rsult with subcompartments of InfectionStateBase with index j and sum up values of subcompartments. + dummy[j] = + result[i] + .segment(Eigen::Index(infectionState.get_firstindex(j)), Eigen::Index(infectionState.get_number(j))) + .sum(); + } + populations.add_time_point(result.get_time(i), dummy); + } + + return populations; +} + +} // namespace lsecir +} // namespace mio \ No newline at end of file diff --git a/cpp/models/lct_secir/model.h b/cpp/models/lct_secir/model.h new file mode 100644 index 0000000000..6778b571dc --- /dev/null +++ b/cpp/models/lct_secir/model.h @@ -0,0 +1,102 @@ +/* +* 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 LCT_SECIR_MODEL_H +#define LCT_SECIR_MODEL_H + +#include "lct_secir/parameters.h" +#include "lct_secir/infection_state.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "memilio/math/eigen.h" +#include + +namespace mio +{ +namespace lsecir +{ +class Model +{ + +public: + /** + * @brief Constructor to create an LCT SECIR Model. + * + * @param[in] init Vector with initial values for all infection states inclusive subcompartments. + * @param[in, out] infectionState_init infectionState for the Model, specifies number of Subcompartments for each infection state. + * @param[in, out] parameters_init Specifies Parameters necessary for the Model. + */ + Model(Eigen::VectorXd init, InfectionState infectionState_init = InfectionState(), + Parameters&& parameters_init = Parameters()); + + /** + * @brief Checks constraints of the model inclusive check for model parameters. + */ + bool check_constraints() const; + + /** + * @brief Evaulates the right-hand-side f of the LCT dydt = f(y, t). + * + * The LCT-SECIR model is defined through ordinary differetial 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] y the current state of the model + * @param[in] t the current time + * @param[out] dydt a reference to the calculated output + */ + void eval_right_hand_side(Eigen::Ref y, ScalarType t, + Eigen::Ref dydt) const; + + /** + * @brief Calculates the population divided in states defined in InfectionStateBase out of a result with subcompartments. + * + * If the model is used for simulation, we will get a result in form of a TimeSeries with infection states divided in Subcompartments. + * Function transforms this TimeSeries in another TimeSeries with just the Basic infectionState. + * This is done by summing up the numbers in the Subcompartments. + * @param[in] result result of a simulation with the model. + * @return result of the simulation divided in the Base infection states. + * Returns TimeSeries with values -1 if calculation is not possible. + */ + TimeSeries calculate_populations(const TimeSeries& result) const; + + /** + * @brief Returns the initial values for the model. + * + * This can be used as initial conditions in an ODE solver. + * @return Vector with initial values for all (sub-)compartments. + */ + Eigen::VectorXd get_initial_values() + { + return m_initial_values; + } + + Parameters parameters{}; ///< Parameters of Model Parameters. + InfectionState infectionState; ///< InfectionState specifies number of subcompartments. + +private: + Eigen::VectorXd m_initial_values; ///< Initial values of the model. + ScalarType m_N0{ + 0}; ///< Total population size at time t_0 for the considered region (inclusive initial value for Dead). +}; + +} // namespace lsecir +} // namespace mio + +#endif // LCTSECIR_MODEL_H \ No newline at end of file diff --git a/cpp/models/lct_secir/parameters.h b/cpp/models/lct_secir/parameters.h new file mode 100644 index 0000000000..24661fc251 --- /dev/null +++ b/cpp/models/lct_secir/parameters.h @@ -0,0 +1,392 @@ +/* +* 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 LCT_SECIR_PARAMS_H +#define LCT_SECIR_PARAMS_H + +#include "memilio/config.h" +#include "memilio/utils/parameter_set.h" +#include "memilio/math/eigen.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/utils/logging.h" + +namespace mio +{ +namespace lsecir +{ + +/********************************************** +* Define Parameters of the LCT-SECIHURD model * +**********************************************/ + +/** + * @brief Average Time spent in the Exposed compartment. + */ +struct TimeExposed { + using Type = ScalarType; + static Type get_default() + { + return 2.0; + } + static std::string name() + { + return "TimeExposed"; + } +}; + +/** + * @brief Average time spent in the TimeInfectedNoSymptoms before developing + * Symptoms or recover in the SECIR model in day unit. + */ +struct TimeInfectedNoSymptoms { + using Type = ScalarType; + static Type get_default() + { + return 1.0; + } + static std::string name() + { + return "TimeInfectedNoSymptoms"; + } +}; + +/** + * @brief Average time spent in the TimeInfectedSymptoms before going to Hospital + * or recover in the SECIR model in day unit. + */ +struct TimeInfectedSymptoms { + using Type = ScalarType; + static Type get_default() + { + return 1.5; + } + static std::string name() + { + return "TimeInfectedNoSymptoms"; + } +}; + +/** + * @brief Average time being in the Hospital before treated by ICU or recover in the + * SECIR model in day unit. + */ +struct TimeInfectedSevere { + using Type = ScalarType; + static Type get_default() + { + return 1.0; + } + static std::string name() + { + return "TimeInfectedSevere"; + } +}; + +/** + * @brief Average time treated by ICU before dead or recover in the SECIR model in day unit. + */ +struct TimeInfectedCritical { + using Type = ScalarType; + static Type get_default() + { + return 1.0; + } + static std::string name() + { + return "TimeInfectedCritical"; + } +}; + +/** + * @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 SECIR model. + */ +struct RiskOfInfectionFromSymptomatic { + using Type = ScalarType; + static Type get_default() + { + return 0.5; + } + static std::string name() + { + return "RiskOfInfectionFromSymptomatic"; + } +}; + +/** + * @brief The percentage of asymptomatic cases in the SECIR model. + */ +struct RecoveredPerInfectedNoSymptoms { + using Type = ScalarType; + static Type get_default() + { + return 0.5; + } + static std::string name() + { + return "RecoveredPerInfectedNoSymptoms"; + } +}; + +/** + * @brief The percentage of hospitalized patients per infected patients in the SECIR model. + */ +struct SeverePerInfectedSymptoms { + using Type = ScalarType; + static Type get_default() + { + return 0.5; + } + static std::string name() + { + return "SeverePerInfectedSymptoms"; + } +}; + +/** + * @brief The percentage of ICU patients per hospitalized patients in the SECIR model. + */ +struct CriticalPerSevere { + using Type = ScalarType; + static Type get_default() + { + return 0.5; + } + static std::string name() + { + return "CriticalPerSevere"; + } +}; + +/** + * @brief The percentage of dead patients per ICU patients in the SECIR model. + */ +struct DeathsPerCritical { + using Type = ScalarType; + static Type get_default() + { + return 0.1; + } + static std::string name() + { + return "DeathsPerCritical"; + } +}; + +/** + * @brief The start day in the LCT 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 LCT-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 LCT-SECIR model. + */ +class Parameters : public ParametersBase +{ +public: + /** + * @brief Default constructor. + */ + Parameters() + : ParametersBase() + { + } + + /** + * @brief checks whether all Parameters satisfy their corresponding constraints and throws errors, if they do not. + * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false. + */ + bool check_constraints() const + { + if (this->get() < 1.0) { + log_error("Constraint check: Parameter TimeExposed is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get() < 1.0) { + log_error("Constraint check: Parameter TimeInfectedNoSymptoms is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get() < 1.0) { + log_error("Constraint check: Parameter TimeInfectedSymptoms is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get() < 1.0) { + log_error("Constraint check: Parameter TimeInfectedSevere is smaller than {:.4f}", 1.0); + return true; + } + + if (this->get() < 1.0) { + log_error("Constraint check: Parameter TimeInfectedCritical is smaller than {:.4f}", 1.0); + return true; + } + + 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() > 1.0) { + log_error("Constraint check: Parameter RecoveredPerInfectedNoSymptoms smaller {:d} or larger {:d}", 0, 1); + return true; + } + + if (this->get() < 0.0 || this->get() > 1.0) { + log_error("Constraint check: Parameter SeverePerInfectedSymptoms smaller {:d} or larger {:d}", 0, 1); + return true; + } + + if (this->get() < 0.0 || this->get() > 1.0) { + log_error("Constraint check: Parameter CriticalPerSevere smaller {:d} or larger {:d}", 0, 1); + return true; + } + + if (this->get() < 0.0 || this->get() > 1.0) { + log_error("Constraint check: Parameter DeathsPerCritical 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; + } + + 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(base, ParametersBase::deserialize(io)); + return success(Parameters(std::move(base))); + } +}; + +} // namespace lsecir +} // namespace mio + +#endif // LCT_SECIR_PARAMS_H \ No newline at end of file diff --git a/cpp/models/lct_secir/simulation.cpp b/cpp/models/lct_secir/simulation.cpp new file mode 100644 index 0000000000..c17dd429c1 --- /dev/null +++ b/cpp/models/lct_secir/simulation.cpp @@ -0,0 +1,58 @@ +/* +* 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 "lct_secir/simulation.h" +#include "lct_secir/model.h" +#include "lct_secir/parameters.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "memilio/math/stepper_wrapper.h" +#include "memilio/math/eigen.h" +#include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp" + +namespace mio +{ +namespace lsecir +{ + +Simulation::Simulation(Model const& model, ScalarType t0, ScalarType dt) + : m_integratorCore( + std::make_shared>()) + , m_model(std::make_unique(model)) + , m_integrator(m_integratorCore) + , m_result(t0, m_model->get_initial_values()) + , m_dt(dt) +{ +} + +TimeSeries simulate(ScalarType t0, ScalarType tmax, ScalarType dt, Model const& model, + std::shared_ptr integrator) +{ + model.check_constraints(); + Simulation sim(model, t0, dt); + if (integrator) { + sim.set_integrator(integrator); + } + sim.advance(tmax); + return sim.get_result(); +} + +} // namespace lsecir +} // namespace mio \ No newline at end of file diff --git a/cpp/models/lct_secir/simulation.h b/cpp/models/lct_secir/simulation.h new file mode 100644 index 0000000000..b66b1e4cb7 --- /dev/null +++ b/cpp/models/lct_secir/simulation.h @@ -0,0 +1,177 @@ +/* +* 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 LCT_SECIR_SIMULATION_H +#define LCT_SECIR_SIMULATION_H + +#include "lct_secir/model.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "memilio/utils/metaprogramming.h" +#include "memilio/math/stepper_wrapper.h" +#include "memilio/math/eigen.h" + +namespace mio +{ +namespace lsecir +{ +/** + * @brief A class for the simulation of a LCT model. + */ +class Simulation +{ +public: + /** + * @brief Set up the simulation with an ODE solver. + * + * Per default Runge Kutta Cash Karp 54 solver is used. + * Use function set_integrator() to set a different integrator. + * @param[in] model An instance of a LCT model. + * @param[in] t0 The Start time, usually 0. + * @param[in] dt Initial step size of integration. + */ + Simulation(Model const& model, ScalarType t0 = 0., ScalarType dt = 0.1); + + /** + * @brief Set the core integrator used in the simulation. + * + * @param[in] integrator Core integrator that should be used for simulation. + */ + void set_integrator(std::shared_ptr integrator) + { + m_integratorCore = std::move(integrator); + m_integrator.set_integrator(m_integratorCore); + } + + /** + * @brief Get a reference to the used integrator. + * + * @return Reference to the core integrator used in the simulation + */ + IntegratorCore& get_integrator() + { + return *m_integratorCore; + } + + /** + * @brief Get a reference to the used integrator. + * + * @return Reference to the core integrator used in the simulation + */ + IntegratorCore const& get_integrator() const + { + return *m_integratorCore; + } + + /** + * @brief Advance simulation to tmax. + * + * tmax must be greater than get_result().get_last_time_point(). + * @param tmax Stopping point of the simulation. + */ + Eigen::Ref advance(ScalarType tmax) + { + return m_integrator.advance( + [this](auto&& y, auto&& t, auto&& dydt) { + get_model().eval_right_hand_side(y, t, dydt); + }, + tmax, m_dt, m_result); + } + + /** + * @brief Get the result of the simulation. + * + * Return the number of persons in all InfectionState%s (inclusive subcompartments). + * @return The result of the simulation in form of a TimeSeries. + */ + TimeSeries& get_result() + { + return m_result; + } + + /** + * @brief Get the result of the simulation. + * + * Return the number of persons in all InfectionState%s (inclusive subcompartments). + * @return The result of the simulation in form of a TimeSeries. + */ + const TimeSeries& get_result() const + { + return m_result; + } + + /** + * @brief Returns a reference to the simulation model used in simulation. + */ + const Model& get_model() const + { + return *m_model; + } + + /** + * @brief Returns a reference to the simulation model used in simulation. + */ + Model& get_model() + { + return *m_model; + } + + /** + * @brief Returns the next time step chosen by integrator. + */ + ScalarType& get_dt() + { + return m_dt; + } + + /** + * @brief Returns the next time step chosen by integrator. + */ + const ScalarType& get_dt() const + { + return m_dt; + } + +private: + std::shared_ptr m_integratorCore; ///< InteratorCore used for Simulation. + std::unique_ptr m_model; ///< LCT-model the simulation should be performed with. + OdeIntegrator m_integrator; ///< OdeIntegrator used to perform simulation. + TimeSeries m_result; ///< The simulation results. + ScalarType m_dt; ///< The time step used (and possibly set) by m_integratorCore::step. +}; + +/** + * @brief Performs a simulation with specified parameters for given model. + * + * @param[in] t0 Start time of the simulation. + * @param[in] tmax End time of the simulation. + * @param[in] dt Initial step size of integration. + * @param[in] model An instance of a LCT model for which the simulation should be performed. + * @param[in] integrator An integrator that should be used to discretize the model equations. + * If default value is used the simulation will be performed with the runge_kutta_cash_karp54 method. + * @return A TimeSeries with the result of the simulation. + */ +TimeSeries simulate(ScalarType t0, ScalarType tmax, ScalarType dt, Model const& model, + std::shared_ptr integrator = nullptr); + +} // namespace lsecir +} // namespace mio + +#endif // LCT_SECIR_SIMULATION_H \ No newline at end of file diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index c19b9fcec2..e30b6f52e2 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -65,6 +65,7 @@ set(TESTSOURCES matchers.h temp_file_register.h sanitizers.cpp + test_lct_secir.cpp ) if(MEMILIO_HAS_JSONCPP) @@ -83,7 +84,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 ide_seir ide_secir abm gtest_main) +target_link_libraries(memilio-test PRIVATE memilio ode_secir ode_seir ode_secirvvs ide_seir ide_secir lct_secir abm gtest_main) 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/data/lct-secir-compartments-compare.csv b/cpp/tests/data/lct-secir-compartments-compare.csv new file mode 100644 index 0000000000..f6c90b6ce8 --- /dev/null +++ b/cpp/tests/data/lct-secir-compartments-compare.csv @@ -0,0 +1,14 @@ +# time | S | E | C | I | H | U | R | D +0.00000000 750.00000000 50.00000000 40.00000000 50.00000000 50.00000000 30.00000000 20.00000000 10.00000000 +0.10000000 748.46824343 50.26284184 39.76391713 50.50314806 49.64881726 29.98773082 21.32228601 10.04301546 +0.29072295 745.55860018 50.65687853 39.32541531 51.51432605 48.99408015 29.94667771 23.87442948 10.12959258 +0.48144590 742.66416112 50.92893708 38.82680068 52.64226763 48.35930851 29.88013044 26.47554743 10.22284710 +0.69951506 739.37529863 51.11312649 38.15900489 54.06952788 47.65846992 29.76917609 29.51661451 10.33878159 +0.94011476 735.77502876 51.18413896 37.33702852 55.75155953 46.91644175 29.59910565 32.95707813 10.47961870 +1.20982607 731.85282708 51.05289643 36.38618932 57.66826350 46.12324676 29.34392108 36.91743422 10.65522161 +1.49825088 728.46517521 50.03000624 35.41581910 59.64744582 45.31873551 28.99124062 41.26645170 10.86512581 +1.85517405 725.99585366 47.01366564 34.32697685 61.88928798 44.38205556 28.43822802 46.79653735 11.15739494 +2.12963864 724.78054190 44.06842389 33.54397453 63.42228150 43.70283544 27.92703310 51.14873009 11.40617955 +2.45390654 723.38317300 40.69846259 32.61593743 65.01365771 42.94248900 27.23273570 56.38817869 11.72536588 +2.89287892 721.51261986 36.46886578 31.26065299 66.80715050 41.97894676 26.15458682 63.62123679 12.19594050 +3.00000000 721.06078616 35.50223737 30.90410239 67.18501599 41.75434540 25.87051417 65.40643018 12.31656834 \ No newline at end of file diff --git a/cpp/tests/data/lct-secir-subcompartments-compare.csv b/cpp/tests/data/lct-secir-subcompartments-compare.csv new file mode 100644 index 0000000000..b1a7a9fa24 --- /dev/null +++ b/cpp/tests/data/lct-secir-subcompartments-compare.csv @@ -0,0 +1,14 @@ +# time | S | E1 | E2 | C1 | C2 | C3 | I | H | U1 | U2 | U3 | U4 | U5 | R | D +0 750.00000000 30.00000000 20.00000000 20.00000000 10.00000000 10.00000000 50.00000000 50.00000000 10.00000000 10.00000000 5.00000000 3.00000000 2.00000000 20.00000000 10.00000000 +0.10000000 748.46824343 29.66723995 20.59560189 18.39292262 11.27497032 10.09602418 50.50314806 49.64881726 9.44659952 9.98074799 5.33955113 3.14782359 2.07300858 21.32228601 10.04301546 +0.29072295 745.55860018 29.07626753 21.58061100 16.00693881 12.71846120 10.60001530 51.51432605 48.99408015 8.49092451 9.85090613 5.91568674 3.46129489 2.22786544 23.87442948 10.12959258 +0.48144590 742.66416112 28.53736239 22.39157469 14.30701953 13.30306581 11.21671534 52.64226763 48.35930851 7.65233409 9.62495724 6.39661907 3.80166291 2.40455713 26.47554743 10.22284710 +0.69951506 739.37529863 27.97565895 23.13746754 12.96585914 13.37340723 11.81973852 54.06952788 47.65846992 6.81828949 9.28178048 6.83204179 4.20428323 2.63278111 29.51661451 10.33878159 +0.94011476 735.77502876 27.41213738 23.77200158 12.00281102 13.08420856 12.25000894 55.75155953 46.91644175 6.03092978 8.83306189 7.17864315 4.64336713 2.91310371 32.95707813 10.47961870 +1.20982607 731.85282708 26.76494962 24.28794681 11.34519845 12.59793545 12.44305542 57.66826350 46.12324676 5.28826093 8.28010536 7.41601376 5.10528556 3.25425548 36.91743422 10.65522161 +1.49825088 728.46517521 25.43323249 24.59677375 10.94270800 12.07556224 12.39754886 59.64744582 45.31873551 4.63012700 7.66567073 7.51659762 5.54181441 3.63703086 41.26645170 10.86512581 +1.85517405 725.99585366 22.53301437 24.48065127 10.65304719 11.53623453 12.13769513 61.88928798 44.38205556 3.97335382 6.91208660 7.46202330 5.97764331 4.11312098 46.79653735 11.15739494 +2.12963864 724.78054190 20.09633142 23.97209247 10.46592242 11.20475910 11.87329300 63.42228150 43.70283544 3.56418839 6.35739777 7.31464373 6.22631004 4.46449317 51.14873009 11.40617955 +2.45390654 723.38317300 17.67435181 23.02411078 10.20372528 10.86770945 11.54450270 65.01365771 42.94248900 3.16689926 5.74360799 7.05333091 6.42218033 4.84671722 56.38817869 11.72536588 +2.89287892 721.51261986 15.06890402 21.39996177 9.73272421 10.42689561 11.10103317 66.80715050 41.97894676 2.74638652 4.99700433 6.59844111 6.52840395 5.28435092 63.62123679 12.19594050 +3.00000000 721.06078616 14.53009114 20.97214623 9.59804298 10.31367815 10.99238125 67.18501599 41.75434540 2.66050262 4.83027461 6.47590662 6.52900540 5.37482492 65.40643018 12.31656834 \ No newline at end of file diff --git a/cpp/tests/test_lct_secir.cpp b/cpp/tests/test_lct_secir.cpp new file mode 100644 index 0000000000..0bd6d397bb --- /dev/null +++ b/cpp/tests/test_lct_secir.cpp @@ -0,0 +1,534 @@ +/* +* 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 "lct_secir/model.h" +#include "lct_secir/infection_state.h" +#include "lct_secir/simulation.h" +#include "lct_secir/parameters.h" +#include "ode_secir/model.h" +#include "memilio/config.h" +#include "memilio/utils/time_series.h" +#include "memilio/epidemiology/uncertain_matrix.h" +#include "memilio/math/eigen.h" +#include "load_test_data.h" + +#include +#include +#include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp" + +// Test confirms that default construction of an LCT model works. +TEST(TestLCTSecir, simulateDefault) +{ + ScalarType t0 = 0; + ScalarType tmax = 1; + ScalarType dt = 0.1; + + Eigen::VectorXd init = Eigen::VectorXd::Constant((int)mio::lsecir::InfectionStateBase::Count, 15); + init[0] = 200; + init[3] = 50; + init[5] = 30; + + mio::lsecir::Model model(init); + mio::TimeSeries result = mio::lsecir::simulate(t0, tmax, dt, model); + + EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); + ScalarType sum_pop = init.sum(); + for (Eigen::Index i = 0; i < result.get_num_time_points(); i++) { + ASSERT_NEAR(sum_pop, result[i].sum(), 1e-5); + } +} + +/* Test compares the result for an LCT SECIR model with one single Subcompartment for each infection state + with the result of the equivalent ODE SECIR model. */ +TEST(TestLCTSecir, compareWithOdeSecir) +{ + ScalarType t0 = 0; + ScalarType tmax = 5; + ScalarType dt = 0.1; + + // Initialization vector for both models. + Eigen::VectorXd init = Eigen::VectorXd::Constant((int)mio::lsecir::InfectionStateBase::Count, 15); + init[0] = 200; + init[3] = 50; + init[5] = 30; + + // Define LCT model. + mio::lsecir::Model model_lct(init); + // Set Parameters. + model_lct.parameters.get() = 2 * 4.2 - 5.2; + model_lct.parameters.get() = 2 * (5.2 - 4.2); + model_lct.parameters.get() = 5.8; + model_lct.parameters.get() = 9.5; + model_lct.parameters.get() = 7.1; + + model_lct.parameters.get() = 0.05; + + mio::ContactMatrixGroup& contact_matrix_lct = model_lct.parameters.get(); + contact_matrix_lct[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10)); + contact_matrix_lct[0].add_damping(0.7, mio::SimulationTime(2.)); + + model_lct.parameters.get() = 0.7; + model_lct.parameters.get() = 0.25; + model_lct.parameters.get() = 50; + model_lct.parameters.get() = 0.1; + model_lct.parameters.get() = 0.09; + model_lct.parameters.get() = 0.2; + model_lct.parameters.get() = 0.25; + model_lct.parameters.get() = 0.3; + + // Simulate. + mio::TimeSeries result_lct = mio::lsecir::simulate( + t0, tmax, dt, model_lct, + std::make_shared>()); + + // Initialize ODE model with one age group. + mio::osecir::Model model_ode(1); + // Set initial distribution of the population. + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = + init[Eigen::Index(mio::lsecir::InfectionStateBase::Exposed)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] = + init[Eigen::Index(mio::lsecir::InfectionStateBase::InfectedNoSymptoms)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptomsConfirmed}] = 0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] = + init[Eigen::Index(mio::lsecir::InfectionStateBase::InfectedSymptoms)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptomsConfirmed}] = 0; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSevere}] = + init[Eigen::Index(mio::lsecir::InfectionStateBase::InfectedSevere)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}] = + init[Eigen::Index(mio::lsecir::InfectionStateBase::InfectedCritical)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Recovered}] = + init[Eigen::Index(mio::lsecir::InfectionStateBase::Recovered)]; + model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Dead}] = + init[Eigen::Index(mio::lsecir::InfectionStateBase::Dead)]; + model_ode.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}, + init.sum()); + + // Set parameters according to the parameters of the LCT model. + // No restrictions by additional parameters. + model_ode.parameters.get() = std::numeric_limits::max(); + model_ode.parameters.get() = std::numeric_limits::max(); + + model_ode.parameters.set(50); + model_ode.parameters.set(0.1); + model_ode.parameters.get()[(mio::AgeGroup)0] = + 5.2; // TimeExposed = 2 * SerialInterval - IncubationTime. + model_ode.parameters.get()[(mio::AgeGroup)0] = + 4.2; // TimeInfectedNoSymptoms = 2* (IncubationTime - SerialInterval). + model_ode.parameters.get()[(mio::AgeGroup)0] = 5.8; + model_ode.parameters.get()[(mio::AgeGroup)0] = 9.5; + model_ode.parameters.get()[(mio::AgeGroup)0] = 7.1; + + mio::ContactMatrixGroup& contact_matrix_ode = model_ode.parameters.get(); + contact_matrix_ode[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 10)); + contact_matrix_ode[0].add_damping(0.7, mio::SimulationTime(2.)); + + model_ode.parameters.get()[(mio::AgeGroup)0] = 0.05; + model_ode.parameters.get()[(mio::AgeGroup)0] = 0.7; + model_ode.parameters.get()[(mio::AgeGroup)0] = 0.09; + model_ode.parameters.get()[(mio::AgeGroup)0] = 0.25; + model_ode.parameters.get()[(mio::AgeGroup)0] = 0.2; + model_ode.parameters.get()[(mio::AgeGroup)0] = 0.25; + model_ode.parameters.get()[(mio::AgeGroup)0] = 0.3; + + // Simulate. + mio::TimeSeries result_ode = + simulate(t0, tmax, dt, model_ode, + std::make_shared>()); + + // Simulation results should be equal. + ASSERT_EQ(result_lct.get_num_time_points(), result_ode.get_num_time_points()); + for (int i = 0; i < 4; ++i) { + ASSERT_NEAR(result_lct.get_time(i), result_ode.get_time(i), 1e-5); + + ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::Susceptible], + result_ode[i][(int)mio::osecir::InfectionState::Susceptible], 1e-5); + ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::Exposed], + result_ode[i][(int)mio::osecir::InfectionState::Exposed], 1e-5); + ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::InfectedNoSymptoms], + result_ode[i][(int)mio::osecir::InfectionState::InfectedNoSymptoms], 1e-5); + ASSERT_NEAR(0, result_ode[i][(int)mio::osecir::InfectionState::InfectedNoSymptomsConfirmed], 1e-5); + ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::InfectedSymptoms], + result_ode[i][(int)mio::osecir::InfectionState::InfectedSymptoms], 1e-5); + ASSERT_NEAR(0, result_ode[i][(int)mio::osecir::InfectionState::InfectedSymptomsConfirmed], 1e-5); + ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::InfectedCritical], + result_ode[i][(int)mio::osecir::InfectionState::InfectedCritical], 1e-5); + ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::InfectedSevere], + result_ode[i][(int)mio::osecir::InfectionState::InfectedSevere], 1e-5); + ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::Recovered], + result_ode[i][(int)mio::osecir::InfectionState::Recovered], 1e-5); + ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::Dead], + result_ode[i][(int)mio::osecir::InfectionState::Dead], 1e-5); + } +} + +// Test if the function eval_right_hand_side() is working using a hand calculated result. +TEST(TestLCTSecir, testEvalRightHandSide) +{ + // Define model. + /* Number of subcompartments, chose more than one subcompartment for all compartments except S, R, D + so that the function is correct for all selections. */ + std::vector SubcompartmentNumbers((int)mio::lsecir::InfectionStateBase::Count, 1); + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Exposed] = 2; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::InfectedNoSymptoms] = 3; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::InfectedSymptoms] = 2; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::InfectedSevere] = 2; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::InfectedCritical] = 2; + mio::lsecir::InfectionState InfState(SubcompartmentNumbers); + + // Define initial population distribution in infection states, one entry per subcompartment. + Eigen::VectorXd init(InfState.get_count()); + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::Susceptible)] = 750; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::Exposed)] = 30; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::Exposed) + 1] = 20; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms)] = 20; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms) + 1] = 10; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms) + 2] = 10; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSymptoms)] = 30; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSymptoms) + 1] = 20; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSevere)] = 40; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSevere) + 1] = 10; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical)] = 10; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 1] = 20; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::Recovered)] = 20; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::Dead)] = 10; + + mio::lsecir::Model model(std::move(init), InfState); + + // Set parameters. + model.parameters.set(2 * 4.2 - 5.2); + model.parameters.get() = 2 * (5.2 - 4.2); + model.parameters.get() = 5.8; + model.parameters.get() = 9.5; + model.parameters.get() = 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.09; + model.parameters.get() = 0.; + model.parameters.get() = 0; + model.parameters.get() = 0.2; + model.parameters.get() = 0.25; + model.parameters.get() = 0.3; + + // Compare the result of eval_right_hand_side() with a hand calculated result. + size_t num_subcompartments = model.infectionState.get_count(); + Eigen::VectorXd dydt(num_subcompartments); + model.eval_right_hand_side(model.get_initial_values(), 0, dydt); + + Eigen::VectorXd compare(num_subcompartments); + compare << -15.3409, -3.4091, 6.25, -17.5, 15, 0, 3.3052, 3.4483, -7.0417, 6.3158, -2.2906, -2.8169, 12.3899, + 1.6901; + + for (size_t i = 0; i < num_subcompartments; i++) { + ASSERT_NEAR(compare[i], dydt[i], 1e-3); + } +} + +// Model setup to compare result with a previous output. +class ModelTestLCTSecir : public testing::Test +{ +protected: + virtual void SetUp() + { + // Define number of subcompartments. + std::vector SubcompartmentNumbers((int)mio::lsecir::InfectionStateBase::Count, 1); + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Exposed] = 2; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::InfectedNoSymptoms] = 3; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::InfectedCritical] = 5; + mio::lsecir::InfectionState InfState(SubcompartmentNumbers); + + // Define initial population distribution in infection states, one entry per Subcompartment. + Eigen::VectorXd init(InfState.get_count()); + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::Susceptible)] = 750; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::Exposed)] = 30; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::Exposed) + 1] = 20; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms)] = 20; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms) + 1] = 10; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms) + 2] = 10; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSymptoms)] = 50; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSevere)] = 50; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical)] = 10; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 1] = 10; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 2] = 5; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 3] = 3; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 4] = 2; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::Recovered)] = 20; + init[InfState.get_firstindex(mio::lsecir::InfectionStateBase::Dead)] = 10; + + // Initialize model and set parameters. + model = new mio::lsecir::Model(std::move(init), InfState); + model->parameters.get() = 2 * 4.2 - 5.2; + model->parameters.get() = 2 * (5.2 - 4.2); + model->parameters.get() = 5.8; + model->parameters.get() = 9.5; + model->parameters.get() = 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.09; + model->parameters.get() = 0.2; + model->parameters.get() = 0.25; + model->parameters.get() = 0.3; + } + + virtual void TearDown() + { + delete model; + } + +public: + mio::lsecir::Model* model = nullptr; +}; + +// Test compares a simulation with the result of a previous run stored in a .csv file. +TEST_F(ModelTestLCTSecir, compareWithPreviousRun) +{ + ScalarType tmax = 3; + mio::TimeSeries result = mio::lsecir::simulate( + 0, tmax, 0.5, *model, + std::make_shared>()); + + // Compare subcompartments. + auto compare = load_test_data_csv("lct-secir-subcompartments-compare.csv"); + + ASSERT_EQ(compare.size(), static_cast(result.get_num_time_points())); + for (size_t i = 0; i < compare.size(); i++) { + ASSERT_EQ(compare[i].size(), static_cast(result.get_num_elements()) + 1) << "at row " << i; + ASSERT_NEAR(result.get_time(i), compare[i][0], 1e-7) << "at row " << i; + for (size_t j = 1; j < compare[i].size(); j++) { + ASSERT_NEAR(result.get_value(i)[j - 1], compare[i][j], 1e-7) << " at row " << i; + } + } + + // Compare base compartments. + mio::TimeSeries population = model->calculate_populations(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; + ASSERT_NEAR(population.get_time(i), compare_population[i][0], 1e-7) << "at row " << i; + for (size_t j = 1; j < compare_population[i].size(); j++) { + ASSERT_NEAR(population.get_value(i)[j - 1], compare_population[i][j], 1e-7) << " at row " << i; + } + } +} + +// Test calculate_populations with a vector of a wrong size. +TEST_F(ModelTestLCTSecir, testCalculatePopWrongSize) +{ + // Deactivate temporarily log output because an error is expected. + mio::set_log_level(mio::LogLevel::off); + mio::TimeSeries init_wrong_size((int)mio::lsecir::InfectionStateBase::Count); + Eigen::VectorXd vec_wrong_size = Eigen::VectorXd::Ones((int)mio::lsecir::InfectionStateBase::Count); + init_wrong_size.add_time_point(-10, vec_wrong_size); + init_wrong_size.add_time_point(-9, vec_wrong_size); + mio::TimeSeries population = model->calculate_populations(init_wrong_size); + EXPECT_EQ(1, population.get_num_time_points()); + for (int i = 0; i < population.get_num_elements(); i++) { + EXPECT_EQ(-1, population.get_last_value()[i]); + } + // Reactive log output. + mio::set_log_level(mio::LogLevel::warn); +} + +// Check constraints of InfectionState and Parameters. +TEST(TestLCTSecir, testConstraints) +{ + // Deactivate temporarily log output for next tests. + mio::set_log_level(mio::LogLevel::off); + + // Check InfectionState with a wrong size of the initialization vector. + std::vector SubcompartmentNumbers1((int)mio::lsecir::InfectionStateBase::Count - 1, 1); + mio::lsecir::InfectionState InfState(SubcompartmentNumbers1); + bool constraint_check = InfState.check_constraints(); + EXPECT_TRUE(constraint_check); + + // Check with right size but wrong number of Subcompartments for Susceptibles. + std::vector SubcompartmentNumbers((int)mio::lsecir::InfectionStateBase::Count, 1); + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Susceptible] = 2; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Exposed] = 2; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::InfectedNoSymptoms] = 3; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::InfectedSymptoms] = 4; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::InfectedSevere] = 5; + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::InfectedCritical] = 6; + InfState.set_SubcompartmentNumbers(SubcompartmentNumbers); + constraint_check = InfState.check_constraints(); + EXPECT_TRUE(constraint_check); + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Susceptible] = 1; + + // Wrong number of subcompartments for Recovered. + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Recovered] = 5; + InfState.set_SubcompartmentNumbers(SubcompartmentNumbers); + constraint_check = InfState.check_constraints(); + EXPECT_TRUE(constraint_check); + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Recovered] = 1; + + // For Dead. + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Dead] = 3; + InfState.set_SubcompartmentNumbers(SubcompartmentNumbers); + constraint_check = InfState.check_constraints(); + EXPECT_TRUE(constraint_check); + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Dead] = 1; + + // Check if number of Subcompartments is zero. + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Exposed] = 0; + InfState.set_SubcompartmentNumbers(SubcompartmentNumbers); + constraint_check = InfState.check_constraints(); + EXPECT_TRUE(constraint_check); + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Exposed] = 2; + + // Check with correct parameters. + InfState.set_SubcompartmentNumbers(SubcompartmentNumbers); + constraint_check = InfState.check_constraints(); + EXPECT_FALSE(constraint_check); + + // Check for exceptions of parameters. + mio::lsecir::Parameters parameters_lct; + parameters_lct.get() = 0; + parameters_lct.get() = 3.1; + parameters_lct.get() = 6.1; + parameters_lct.get() = 11.1; + parameters_lct.get() = 17.1; + parameters_lct.get() = 0.01; + mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, 1); + parameters_lct.get() = mio::UncertainContactMatrix(contact_matrix); + + parameters_lct.get() = 1; + parameters_lct.get() = 1; + parameters_lct.get() = 0.1; + parameters_lct.get() = 0.1; + parameters_lct.get() = 0.1; + parameters_lct.get() = 0.1; + + // Check improper TimeExposed. + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 3.1; + + // Check TimeInfectedNoSymptoms. + parameters_lct.get() = 0.1; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 3.1; + + // Check TimeInfectedSymptoms. + parameters_lct.get() = -0.1; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 6.1; + + // Check TimeInfectedSevere. + parameters_lct.get() = 0.5; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 11.1; + + // Check TimeInfectedCritical. + parameters_lct.get() = 0.; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 17.1; + + // Check TransmissionProbabilityOnContact. + parameters_lct.get() = -1; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 0.01; + + // Check RelativeTransmissionNoSymptoms. + parameters_lct.get() = 1.5; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 1; + + // Check RiskOfInfectionFromSymptomatic. + parameters_lct.get() = 1.5; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 1; + + // Check RecoveredPerInfectedNoSymptoms. + parameters_lct.get() = 1.5; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 0.1; + + // Check SeverePerInfectedSymptoms. + parameters_lct.get() = -1; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 0.1; + + // Check CriticalPerSevere. + parameters_lct.get() = -1; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 0.1; + + // Check DeathsPerCritical. + parameters_lct.get() = -1; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 0.1; + + // Check Seasonality. + parameters_lct.get() = 1; + constraint_check = parameters_lct.check_constraints(); + EXPECT_TRUE(constraint_check); + parameters_lct.get() = 0.1; + + // Check with correct parameters. + constraint_check = parameters_lct.check_constraints(); + EXPECT_FALSE(constraint_check); + + // Check for model. + // Check wrong size of initial value vector. + mio::lsecir::Model model1(std::move(Eigen::VectorXd::Ones(InfState.get_count() - 1)), InfState, + std::move(parameters_lct)); + constraint_check = model1.check_constraints(); + EXPECT_TRUE(constraint_check); + + // Check with values smaller than zero. + mio::lsecir::Model model2(std::move(Eigen::VectorXd::Constant(InfState.get_count(), -1)), InfState, + std::move(parameters_lct)); + constraint_check = model2.check_constraints(); + EXPECT_TRUE(constraint_check); + + // Check with correct conditions. + mio::lsecir::Model model3(std::move(Eigen::VectorXd::Constant(InfState.get_count(), 100)), InfState, + std::move(parameters_lct)); + constraint_check = model3.check_constraints(); + EXPECT_FALSE(constraint_check); + + // Reactive log output. + mio::set_log_level(mio::LogLevel::warn); +} \ No newline at end of file From dee2f39c7b58d7ab8e41c735cfb90ff2c17b9852 Mon Sep 17 00:00:00 2001 From: lenaploetzke Date: Mon, 5 Feb 2024 10:19:30 +0100 Subject: [PATCH 02/10] Added ConfirmedCasesNoAgeEntry to epi_data --- cpp/memilio/io/epi_data.h | 52 +++++++++++++++ cpp/tests/data/cases_all_germany.json | 92 +++++++++++++++++++++++++++ cpp/tests/test_epi_data_io.cpp | 58 +++++++++++++++++ 3 files changed, 202 insertions(+) create mode 100644 cpp/tests/data/cases_all_germany.json diff --git a/cpp/memilio/io/epi_data.h b/cpp/memilio/io/epi_data.h index 587d3e4761..fd1c4d024d 100644 --- a/cpp/memilio/io/epi_data.h +++ b/cpp/memilio/io/epi_data.h @@ -65,6 +65,58 @@ class StringDate : public Date } }; +/** + * Represents the entries of a confirmed cases data file, e.g., from RKI. + * Number of confirmed, recovered and deceased in a Germany on a specific date. + * ConfirmedCasesNoAgeEntry is a simplified version of ConfirmedCasesDataEntry without agegroups. + */ +class ConfirmedCasesNoAgeEntry +{ +public: + double num_confirmed; + double num_recovered; + double num_deaths; + Date date; + + template + static IOResult deserialize(IOContext& io) + { + auto obj = io.expect_object("ConfirmedCasesNoAgeEntry"); + auto num_confirmed = obj.expect_element("Confirmed", Tag{}); + auto num_recovered = obj.expect_element("Recovered", Tag{}); + auto num_deaths = obj.expect_element("Deaths", Tag{}); + auto date = obj.expect_element("Date", Tag{}); + return apply( + io, + [](auto&& nc, auto&& nr, auto&& nd, auto&& d) { + return ConfirmedCasesNoAgeEntry{nc, nr, nd, d}; + }, + num_confirmed, num_recovered, num_deaths, date); + } +}; + +/** + * Deserialize a list of ConfirmedCasesNoAgeEntry from json. + * @param jsvalue json value, must be an array of objects, objects must match ConfirmedCasesNoAgeEntry. + * Exactly one entry per date should be provided, for example no entries per age group. + * @return list of ConfirmedCasesNoAgeEntry. + */ +inline IOResult> deserialize_confirmed_cases_noage(const Json::Value& jsvalue) +{ + return deserialize_json(jsvalue, Tag>{}); +} + +/** + * Deserialize a list of ConfirmedCasesNoAgeEntry from json. + * @param jsvalue json value, must be an array of objects, objects must match ConfirmedCasesNoAgeEntry. + * Exactly one entry per date should be provided, for example no entries per age group + * @return list of ConfirmedCasesNoAgeEntry. + */ +inline IOResult> read_confirmed_cases_noage(const std::string& filename) +{ + return read_json(filename, Tag>{}); +} + /** * Represents the entries of a confirmed cases data file, e.g., from RKI. * Number of confirmed, recovered and deceased in a region on a specific date. diff --git a/cpp/tests/data/cases_all_germany.json b/cpp/tests/data/cases_all_germany.json new file mode 100644 index 0000000000..f4859ae673 --- /dev/null +++ b/cpp/tests/data/cases_all_germany.json @@ -0,0 +1,92 @@ +[ + { + "Date": "2020-05-24", + "Confirmed": 1.2, + "Deaths": 0.5, + "Recovered": 0 + }, + { + "Date": "2020-05-25", + "Confirmed": 5, + "Deaths": 2.0, + "Recovered": 2.5 + }, + { + "Date": "2020-05-26", + "Confirmed": 20.2, + "Deaths": 5.01, + "Recovered": 3.0 + }, + { + "Date": "2020-05-27", + "Confirmed": 25.6, + "Deaths": 6, + "Recovered": 3.0 + }, + { + "Date": "2020-05-28", + "Confirmed": 27.6, + "Deaths": 6, + "Recovered": 3.0 + }, + { + "Date": "2020-05-29", + "Confirmed": 30.6, + "Deaths": 6, + "Recovered": 3.0 + }, + { + "Date": "2020-05-30", + "Confirmed": 30.6, + "Deaths": 6, + "Recovered": 3.0 + }, + { + "Date": "2020-05-31", + "Confirmed": 35, + "Deaths": 6, + "Recovered": 3.0 + }, + { + "Date": "2020-06-01", + "Confirmed": 44, + "Deaths": 8, + "Recovered": 3.0 + }, + { + "Date": "2020-06-02", + "Confirmed": 44.5, + "Deaths": 8.999, + "Recovered": 6 + }, + { + "Date": "2020-06-03", + "Confirmed": 70, + "Deaths": 9, + "Recovered": 6 + }, + { + "Date": "2020-06-04", + "Confirmed": 100, + "Deaths": 20, + "Recovered": 6 + }, + { + "Date": "2020-06-05", + "Confirmed": 100.3, + "Deaths": 23.56, + "Recovered": 6 + }, + { + "Date": "2020-06-06", + "Confirmed": 115, + "Deaths": 24, + "Recovered": 6 + }, + { + "Date": "2020-06-07", + "Confirmed": 120.6, + "Deaths": 30, + "Recovered": 6 + } +] \ No newline at end of file diff --git a/cpp/tests/test_epi_data_io.cpp b/cpp/tests/test_epi_data_io.cpp index 7f197fa4a3..69334dafce 100644 --- a/cpp/tests/test_epi_data_io.cpp +++ b/cpp/tests/test_epi_data_io.cpp @@ -94,6 +94,36 @@ TEST(TestEpiDataIo, read_rki_error_age) ASSERT_THAT(print_wrap(result), IsFailure(mio::StatusCode::InvalidValue)); } +TEST(TestEpiDataIo, read_confirmed_cases_noage) +{ + Json::Value js(Json::arrayValue); + js[0]["Date"] = "2021-12-01"; + js[0]["Confirmed"] = 1; + js[0]["Deaths"] = 2; + js[0]["Recovered"] = 3; + + js[1]["Date"] = "2021-12-02"; + js[1]["Confirmed"] = 4; + js[1]["Deaths"] = 5; + js[1]["Recovered"] = 6; + + auto result = mio::deserialize_confirmed_cases_noage(js); + ASSERT_THAT(print_wrap(result), IsSuccess()); + + auto rki_data_noage = result.value(); + ASSERT_EQ(rki_data_noage.size(), 2); + + ASSERT_EQ(rki_data_noage[0].date, mio::Date(2021, 12, 1)); + ASSERT_EQ(rki_data_noage[0].num_confirmed, 1); + ASSERT_EQ(rki_data_noage[0].num_deaths, 2); + ASSERT_EQ(rki_data_noage[0].num_recovered, 3); + + ASSERT_EQ(rki_data_noage[1].date, mio::Date(2021, 12, 2)); + ASSERT_EQ(rki_data_noage[1].num_confirmed, 4); + ASSERT_EQ(rki_data_noage[1].num_deaths, 5); + ASSERT_EQ(rki_data_noage[1].num_recovered, 6); +} + TEST(TestEpiDataIo, read_divi) { Json::Value js(Json::arrayValue); @@ -301,6 +331,34 @@ TEST(TestEpiDataIo, read_confirmed_cases_data) ASSERT_EQ(case_data[2].state_id, boost::none); } +TEST(TestEpiDataIo, read_confirmed_cases_noage_data) +{ + auto rki_data_noage = + mio::read_confirmed_cases_noage(mio::path_join(TEST_DATA_DIR, "cases_all_germany.json")).value(); + + ASSERT_EQ(rki_data_noage.size(), 15); + + ASSERT_EQ(rki_data_noage[0].date, mio::Date(2020, 05, 24)); + ASSERT_EQ(rki_data_noage[0].num_confirmed, 1.2); + ASSERT_EQ(rki_data_noage[0].num_deaths, 0.5); + ASSERT_EQ(rki_data_noage[0].num_recovered, 0); + + ASSERT_EQ(rki_data_noage[8].date, mio::Date(2020, 06, 01)); + ASSERT_EQ(rki_data_noage[8].num_confirmed, 44); + ASSERT_EQ(rki_data_noage[8].num_deaths, 8); + ASSERT_EQ(rki_data_noage[8].num_recovered, 3); + + ASSERT_EQ(rki_data_noage[9].date, mio::Date(2020, 06, 02)); + ASSERT_EQ(rki_data_noage[9].num_confirmed, 44.5); + ASSERT_EQ(rki_data_noage[9].num_deaths, 8.999); + ASSERT_EQ(rki_data_noage[9].num_recovered, 6); + + ASSERT_EQ(rki_data_noage[14].date, mio::Date(2020, 06, 07)); + ASSERT_EQ(rki_data_noage[14].num_confirmed, 120.6); + ASSERT_EQ(rki_data_noage[14].num_deaths, 30); + ASSERT_EQ(rki_data_noage[14].num_recovered, 6); +} + TEST(TestEpiDataIO, read_vaccination_data) { auto vacc_data = mio::read_vaccination_data(mio::path_join(TEST_DATA_DIR, "test_all_ageinf_vacc.json")).value(); From 21f3e4c9de8fd0727a443901dc82f31236664e44 Mon Sep 17 00:00:00 2001 From: lenaploetzke Date: Mon, 5 Feb 2024 14:07:57 +0100 Subject: [PATCH 03/10] added functionality to initialize LCT with real data --- cpp/models/lct_secir/CMakeLists.txt | 2 + cpp/models/lct_secir/README.md | 22 +- cpp/models/lct_secir/parameters_io.cpp | 275 +++++++++++++++++++++++++ cpp/models/lct_secir/parameters_io.h | 67 ++++++ cpp/tests/CMakeLists.txt | 1 + cpp/tests/data/test_empty_file.json | 1 + cpp/tests/test_lct_parameters_io.cpp | 145 +++++++++++++ 7 files changed, 504 insertions(+), 9 deletions(-) create mode 100644 cpp/models/lct_secir/parameters_io.cpp create mode 100644 cpp/models/lct_secir/parameters_io.h create mode 100644 cpp/tests/data/test_empty_file.json create mode 100644 cpp/tests/test_lct_parameters_io.cpp diff --git a/cpp/models/lct_secir/CMakeLists.txt b/cpp/models/lct_secir/CMakeLists.txt index e93fbe8bed..f69a80acfa 100644 --- a/cpp/models/lct_secir/CMakeLists.txt +++ b/cpp/models/lct_secir/CMakeLists.txt @@ -5,6 +5,8 @@ add_library(lct_secir simulation.h simulation.cpp parameters.h + parameters_io.h + parameters_io.cpp ) target_link_libraries(lct_secir PUBLIC memilio) target_include_directories(lct_secir PUBLIC diff --git a/cpp/models/lct_secir/README.md b/cpp/models/lct_secir/README.md index 462fad5626..1b14152a70 100644 --- a/cpp/models/lct_secir/README.md +++ b/cpp/models/lct_secir/README.md @@ -1,23 +1,27 @@ # LCT SECIR model This model is based on the Linear Chain Trick. +The Linear Chain Trick provides the Option to use Erlang-distributed sojourn times in the compartments through the use of subcompartments. +The normal ODE models have (possibly unrealistic) exponentially distributed sojourn times. +The LCT model can still be described by an ordinary differential equation system. For the concept see - Lena Plötzke, "Der Linear Chain Trick in der epidemiologischen Modellierung als Kompromiss zwischen gewöhnlichen und Integro-Differentialgleichungen", 2023. (https://elib.dlr.de/200381/, german) - P. J. Hurtado und A. S. Kirosingh, "Generalizations of the ‘Linear Chain Trick’: incorporating more flexible dwell time distributions into mean field ODE models“, 2019. (https://doi.org/10.1007/s00285-019-01412-w) The eight compartments -- Susceptible, may become exposed at any time -- Exposed, becomes infected after some time -- InfectedNoSymptoms, becomes InfectedSymptoms or Recovered after some time -- InfectedSymptoms, becomes InfectedSevere or Recovered after some time -- InfectedSevere, becomes InfectedCritical or Recovered after some time -- InfectedCritical, becomes Recovered or Dead after some time -- Recovered -- Dead +- `Susceptible`, may become exposed at any time +- `Exposed`, becomes infected after some time +- `InfectedNoSymptoms`, becomes InfectedSymptoms or Recovered after some time +- `InfectedSymptoms`, becomes InfectedSevere or Recovered after some time +- `InfectedSevere`, becomes InfectedCritical or Recovered after some time +- `InfectedCritical`, becomes Recovered or Dead after some time +- `Recovered` +- `Dead` are used to simulate the spread of the disease. It is possible to include subcompartments for the five compartments Exposed, InfectedNoSymptoms, InfectedSymptoms, InfectedSevere and InfectedCritical. -A simple example can be found at: examples/lct_secir.cpp. +A simple example can be found at [LCT minimal example](../../examples/lct_secir.cpp). +The file [parameters_io](parameters_io.h) provides functionality to compute an initial value vector for the LCT-SECIR model based on real data. diff --git a/cpp/models/lct_secir/parameters_io.cpp b/cpp/models/lct_secir/parameters_io.cpp new file mode 100644 index 0000000000..2b3259d0e0 --- /dev/null +++ b/cpp/models/lct_secir/parameters_io.cpp @@ -0,0 +1,275 @@ +/* +* 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 "lct_secir/parameters_io.h" +#include "memilio/config.h" + +#ifdef MEMILIO_HAS_JSONCPP + +#include "lct_secir/infection_state.h" +#include "lct_secir/parameters.h" +#include "memilio/io/epi_data.h" +#include "memilio/io/io.h" +#include "memilio/utils/date.h" +#include "memilio/utils/logging.h" +#include "memilio/math/eigen.h" + +#include + +namespace mio +{ +namespace lsecir +{ + +IOResult get_initial_data_from_file(std::string const& path, Date date, InfectionState infectionState, + Parameters&& parameters, ScalarType total_population, + ScalarType scale_confirmed_cases) +{ + // Try to get rki data from path. + BOOST_OUTCOME_TRY(rki_data, mio::read_confirmed_cases_noage(path)); + auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { + return a.date < b.date; + }); + if (max_date_entry == rki_data.end()) { + log_error("RKI data file is empty."); + return failure(StatusCode::InvalidFileFormat, path + ", file is empty."); + } + auto max_date = max_date_entry->date; + if (max_date < date) { + log_error("Specified date does not exist in RKI data."); + return failure(StatusCode::OutOfRange, path + ", specified date does not exist in RKI data."); + } + // Compute initial values for all subcompartments. + Eigen::VectorXd init = Eigen::VectorXd::Zero(infectionState.get_count()); + // Define variables for parameters that are often needed. + ScalarType timeExposed = parameters.get(); + ScalarType timeInfectedNoSymptoms = parameters.get(); + ScalarType timeInfectedSymptoms = parameters.get(); + ScalarType timeInfectedSevere = parameters.get(); + ScalarType timeInfectedCritical = parameters.get(); + ScalarType scale_mu_CR = 1 / (1 - parameters.get()); + + ScalarType min_offset_needed = std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical); + ScalarType max_offset_needed = std::ceil(timeExposed + timeInfectedNoSymptoms); + + bool min_offset_needed_avail = false; + bool max_offset_needed_avail = false; + + // Go through the entries of rki_data and check if date is needed for calculation. Confirmed cases are scaled. + for (auto&& entry : rki_data) { + auto offset = get_offset_in_days(entry.date, date); + if (!(offset < min_offset_needed) || !(offset > max_offset_needed)) { + // Add confirmed cases at date to compartment Recovered. + if (offset == 0) { + init[infectionState.get_firstindex(InfectionStateBase::Recovered)] += + scale_confirmed_cases * entry.num_confirmed; + } + + // Compute initial values for compartment InfectedNoSymptoms. + if (offset >= 0 && offset <= std::ceil(timeInfectedNoSymptoms)) { + ScalarType TCi = + timeInfectedNoSymptoms / infectionState.get_number(InfectionStateBase::InfectedNoSymptoms); + int idx_Cn = infectionState.get_firstindex(InfectionStateBase::InfectedNoSymptoms) + + infectionState.get_number(InfectionStateBase::InfectedNoSymptoms) - 1; + for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedNoSymptoms); i++) { + if (offset == std::floor(i * TCi)) { + init[idx_Cn - i] -= (1 - (i * TCi - std::floor(i * TCi))) * scale_mu_CR * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(i * TCi)) { + init[idx_Cn - i] -= + (i * TCi - std::floor(i * TCi)) * scale_mu_CR * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::floor((i + 1) * TCi)) { + init[idx_Cn - i] += (1 - ((i + 1) * TCi - std::floor((i + 1) * TCi))) * scale_mu_CR * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil((i + 1) * TCi)) { + init[idx_Cn - i] += ((i + 1) * TCi - std::floor((i + 1) * TCi)) * scale_mu_CR * + scale_confirmed_cases * entry.num_confirmed; + } + } + } + + // Compute initial values for compartment Exposed. + if (offset >= std::floor(timeInfectedNoSymptoms) && offset <= max_offset_needed) { + ScalarType TEi = timeExposed / infectionState.get_number(InfectionStateBase::Exposed); + int idx_En = infectionState.get_firstindex(InfectionStateBase::Exposed) + + infectionState.get_number(InfectionStateBase::Exposed) - 1; + for (int i = 0; i < infectionState.get_number(InfectionStateBase::Exposed); i++) { + if (offset == std::floor(timeInfectedNoSymptoms + i * TEi)) { + init[idx_En - i] -= + (1 - (timeInfectedNoSymptoms + i * TEi - std::floor(timeInfectedNoSymptoms + i * TEi))) * + scale_mu_CR * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(timeInfectedNoSymptoms + i * TEi)) { + init[idx_En - i] -= + (timeInfectedNoSymptoms + i * TEi - std::floor(timeInfectedNoSymptoms + i * TEi)) * + scale_mu_CR * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::floor(timeInfectedNoSymptoms + (i + 1) * TEi)) { + init[idx_En - i] += (1 - (timeInfectedNoSymptoms + (i + 1) * TEi - + std::floor(timeInfectedNoSymptoms + (i + 1) * TEi))) * + scale_mu_CR * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(timeInfectedNoSymptoms + (i + 1) * TEi)) { + init[idx_En - i] += (timeInfectedNoSymptoms + (i + 1) * TEi - + std::floor(timeInfectedNoSymptoms + (i + 1) * TEi)) * + scale_mu_CR * scale_confirmed_cases * entry.num_confirmed; + } + } + } + + // Compute initial values for compartment InfectedSymptoms. + if (offset >= std::floor(-timeInfectedSymptoms) && offset <= 0) { + ScalarType TIi = timeInfectedSymptoms / infectionState.get_number(InfectionStateBase::InfectedSymptoms); + int idx_I1 = infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms); + for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedSymptoms); i++) { + if (offset == std::floor(-TIi * (i + 1))) { + init[idx_I1 + i] -= (1 - (-(i + 1) * TIi - std::floor(-(i + 1) * TIi))) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-TIi * (i + 1))) { + init[idx_I1 + i] -= + (-(i + 1) * TIi - std::floor(-(i + 1) * TIi)) * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::floor(-TIi * i)) { + init[idx_I1 + i] += + (1 - (-i * TIi - std::floor(-i * TIi))) * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-TIi * i)) { + init[idx_I1 + i] += + (-i * TIi - std::floor(-i * TIi)) * scale_confirmed_cases * entry.num_confirmed; + } + } + } + + // Compute initial values for compartment InfectedSevere. + if (offset >= std::floor(-timeInfectedSymptoms - timeInfectedSevere) && + offset <= std::ceil(-timeInfectedSymptoms)) { + ScalarType THi = timeInfectedSevere / infectionState.get_number(InfectionStateBase::InfectedSevere); + ScalarType mu_IH = parameters.get(); + int idx_H1 = infectionState.get_firstindex(InfectionStateBase::InfectedSevere); + for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedSevere); i++) { + if (offset == std::floor(-timeInfectedSymptoms - THi * (i + 1))) { + init[idx_H1 + i] -= mu_IH * + (1 - (-timeInfectedSymptoms - (i + 1) * THi - + std::floor(-timeInfectedSymptoms - (i + 1) * THi))) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-timeInfectedSymptoms - THi * (i + 1))) { + init[idx_H1 + i] -= mu_IH * + (-timeInfectedSymptoms - (i + 1) * THi - + std::floor(-timeInfectedSymptoms - (i + 1) * THi)) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::floor(-timeInfectedSymptoms - THi * i)) { + init[idx_H1 + i] += + mu_IH * + (1 - (-timeInfectedSymptoms - i * THi - std::floor(-timeInfectedSymptoms - i * THi))) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-timeInfectedSymptoms - THi * i)) { + init[idx_H1 + i] += + mu_IH * (-timeInfectedSymptoms - i * THi - std::floor(-timeInfectedSymptoms - i * THi)) * + scale_confirmed_cases * entry.num_confirmed; + } + } + } + + // Compute initial values for compartment InfectedCritical. + if (offset >= min_offset_needed && offset <= std::ceil(-timeInfectedSymptoms - timeInfectedSevere)) { + ScalarType TUi = timeInfectedCritical / infectionState.get_number(InfectionStateBase::InfectedCritical); + ScalarType mu_IH = parameters.get(); + ScalarType mu_HU = parameters.get(); + int idx_U1 = infectionState.get_firstindex(InfectionStateBase::InfectedCritical); + for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedCritical); i++) { + if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere - TUi * (i + 1))) { + init[idx_U1 + i] -= + mu_IH * mu_HU * + (1 - (-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi))) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - TUi * (i + 1))) { + init[idx_U1 + i] -= mu_IH * mu_HU * + (-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi)) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere - TUi * i)) { + init[idx_U1 + i] += mu_IH * mu_HU * + (1 - (-timeInfectedSymptoms - timeInfectedSevere - i * TUi - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - i * TUi))) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - TUi * i)) { + init[idx_U1 + i] += mu_IH * mu_HU * + (-timeInfectedSymptoms - timeInfectedSevere - i * TUi - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - i * TUi)) * + scale_confirmed_cases * entry.num_confirmed; + } + } + } + + // Compute Dead. + if (offset == min_offset_needed) { + min_offset_needed_avail = true; + init[infectionState.get_firstindex(InfectionStateBase::Dead)] += + (1 - (-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical))) * + entry.num_deaths; + } + if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical)) { + init[infectionState.get_firstindex(InfectionStateBase::Dead)] += + (-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical)) * + entry.num_deaths; + } + if (offset == max_offset_needed) { + max_offset_needed_avail = true; + } + } + } + + // Compute Recovered. + init[infectionState.get_firstindex(InfectionStateBase::Recovered)] -= + init.segment(infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms), + infectionState.get_number(InfectionStateBase::InfectedSymptoms) + + infectionState.get_number(InfectionStateBase::InfectedSevere) + + infectionState.get_number(InfectionStateBase::InfectedCritical)) + .sum(); + init[infectionState.get_firstindex(InfectionStateBase::Recovered)] -= + init[infectionState.get_firstindex(InfectionStateBase::Dead)]; + + // Compute Susceptibles + init[infectionState.get_firstindex(InfectionStateBase::Susceptible)] = + total_population - init.segment(1, infectionState.get_count() - 1).sum(); + if (!max_offset_needed_avail || !min_offset_needed_avail) { + log_error("Necessary range of dates needed to compute initial values does not exist in RKI data."); + return failure(StatusCode::OutOfRange, path + ", necessary range of dates does not exist in RKI data."); + } + + return init; +} + +} // namespace lsecir +} // namespace mio +#endif // MEMILIO_HAS_JSONCPP \ No newline at end of file diff --git a/cpp/models/lct_secir/parameters_io.h b/cpp/models/lct_secir/parameters_io.h new file mode 100644 index 0000000000..b52b6cc20e --- /dev/null +++ b/cpp/models/lct_secir/parameters_io.h @@ -0,0 +1,67 @@ +/* +* 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 LCTSECIR_PARAMETERS_IO_H +#define LCTSECIR_PARAMETERS_IO_H + +#include "memilio/config.h" + +#ifdef MEMILIO_HAS_JSONCPP + +#include "lct_secir/parameters.h" +#include "lct_secir/infection_state.h" +#include "memilio/math/eigen.h" + +#include + +namespace mio +{ +namespace lsecir +{ + +/** +* @brief Computes an initialization vector for a LCT model with data from RKI. +* +* Computes an initial value vector for an LCT model with the defined infectionState and the parameters. +* For the computation expected sojourntime in the subcompartments are used. To calculate the initial values, +* we assume for simplicity that individuals stay in the subcompartment for exactly the expected time. +* The RKI data are linearly interpolated within one day. +* The RKI data should contain data for each needed day without division of age groups, the completeness of the dates is not verified. +* Data could be downloaded eg with the file pycode/memilio-epidata/memilio/epidata/getCaseData.py, +* which creates a file named cases_all_germany.json or a similar name. One should set impute_dates=True so that missing dates are imputed. +* +* @param path Path to the RKI file. +* @param date Date for which the initial values should be computed. date is the start date of the simulation. +* @param infectionState InfectionState used to calculate the initial values. Defines eg the size of the calculated vector. +* @param parameters Parameters used to calculate the initial values. Defines eg the expected sojourntimes. +* @param total_population Total size of the population of the considered region. +* The sum of all values in the returned vector will be equal to that value. +* @param scale_confirmed_cases Factor by which to scale the confirmed cases of rki data to consider unreported cases. +* @returns Initialization Vector or any io errors that happen during reading of the files. +*/ +IOResult get_initial_data_from_file(std::string const& path, Date date, InfectionState infectionState, + Parameters&& parameters, ScalarType total_population, + ScalarType scale_confirmed_cases = 1.); + +} // namespace lsecir +} // namespace mio +#endif // MEMILIO_HAS_JSONCPP + +#endif // LCTSECIR_PARAMETERS_IO_H \ No newline at end of file diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index e30b6f52e2..e23087a860 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -66,6 +66,7 @@ set(TESTSOURCES temp_file_register.h sanitizers.cpp test_lct_secir.cpp + test_lct_parameters_io.cpp ) if(MEMILIO_HAS_JSONCPP) diff --git a/cpp/tests/data/test_empty_file.json b/cpp/tests/data/test_empty_file.json new file mode 100644 index 0000000000..0637a088a0 --- /dev/null +++ b/cpp/tests/data/test_empty_file.json @@ -0,0 +1 @@ +[] \ No newline at end of file diff --git a/cpp/tests/test_lct_parameters_io.cpp b/cpp/tests/test_lct_parameters_io.cpp new file mode 100644 index 0000000000..615996e388 --- /dev/null +++ b/cpp/tests/test_lct_parameters_io.cpp @@ -0,0 +1,145 @@ +/* +* 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 "memilio/config.h" + +#include "lct_secir/parameters_io.h" +#include "lct_secir/parameters.h" +#include "lct_secir/infection_state.h" +#include "memilio/math/eigen.h" +#include "memilio/utils/date.h" +#include "test_data_dir.h" +#include "memilio/io/io.h" +#include + +#include + +// Check that Initialization based on synthetic RKI data match previous result. +TEST(TestLCTParametersIo, ReadPopulationDataRKI) +{ + ScalarType total_population = 1000.0; + auto start_date = mio::Date(2020, 6, 1); + + // Define parameters used for simulation and initialization. + mio::lsecir::Parameters parameters; + parameters.get() = 2.3; + parameters.get() = 3.3; + parameters.get() = 2.4; + parameters.get() = 1.8; + parameters.get() = 3.0; + + parameters.get() = 0.2; + parameters.get() = 0.1; + parameters.get() = 0.3; + parameters.get() = 0.2; + + // Define number of subcompartments. + std::vector vec_subcompartments((int)mio::lsecir::InfectionStateBase::Count, 1); + // Use subcompartments with a soujourn time of approximately one day in each subcompartment. + vec_subcompartments[(int)mio::lsecir::InfectionStateBase::Exposed] = + (int)round(parameters.get()); + vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedNoSymptoms] = + (int)round(parameters.get()); + vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedSymptoms] = + (int)round(parameters.get()); + vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedSevere] = + (int)round(parameters.get()); + // Both realistic distributions for times corresponding to InfectedCritical of the IDE model are exponential distributions. + vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedCritical] = 1; + mio::lsecir::InfectionState infectionState(vec_subcompartments); + + // Calculate initial value vector for subcompartments with RKI data. + auto read_result = + mio::lsecir::get_initial_data_from_file(mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date, + infectionState, std::move(parameters), total_population, 1.); + + ASSERT_THAT(print_wrap(read_result), IsSuccess()); + + auto init_subcompartments = read_result.value(); + // Previous result. + Eigen::VectorXd compare(infectionState.get_count()); + compare << 863.05, 14.30625, 8.53125, 30.1125, 36.1875, 3.8125, 9.88, 3.52, 0.09, 0.25, 0.6888, 27.8712, 1.7; + + for (int i = 0; i < infectionState.get_count(); i++) { + EXPECT_NEAR(init_subcompartments[i], compare[i], 1e-4) << "at subcompartment number " << i; + } +} + +// Check some cases where computation of initial values for an LCT model based on RKI data should fail. +TEST(TestLCTParametersIo, ReadPopulationDataRKIFailure) +{ + ScalarType total_population = 1000.0; + + // Define parameters used for simulation and initialization. + mio::lsecir::Parameters parameters; + parameters.get() = 2.3; + parameters.get() = 3.3; + parameters.get() = 2.4; + parameters.get() = 1.8; + parameters.get() = 3.0; + + parameters.get() = 0.2; + parameters.get() = 0.08; + parameters.get() = 0.15; + parameters.get() = 0.22; + + // Define number of subcompartments. + std::vector vec_subcompartments((int)mio::lsecir::InfectionStateBase::Count, 1); + // Use subcompartments with a soujourn time of approximately one day in each subcompartment. + vec_subcompartments[(int)mio::lsecir::InfectionStateBase::Exposed] = + (int)round(parameters.get()); + vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedNoSymptoms] = + (int)round(parameters.get()); + vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedSymptoms] = + (int)round(parameters.get()); + vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedSevere] = + (int)round(parameters.get()); + // Both realistic distributions for times corresponding to InfectedCritical of the IDE model are exponential distributions. + vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedCritical] = 1; + mio::lsecir::InfectionState infectionState(vec_subcompartments); + + // Deactivate temporarily log output for next tests. + mio::set_log_level(mio::LogLevel::off); + + // Case where start_date is later than maximal provided date in file. + auto start_date = mio::Date(2020, 6, 9); + auto read_result1 = + mio::lsecir::get_initial_data_from_file(mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date, + infectionState, std::move(parameters), total_population, 1.); + + ASSERT_THAT(print_wrap(read_result1), IsFailure(mio::StatusCode::OutOfRange)); + // Case where not all needed dates are provided. + start_date = mio::Date(2020, 6, 6); + auto read_result2 = + mio::lsecir::get_initial_data_from_file(mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date, + infectionState, std::move(parameters), total_population, 1.); + + ASSERT_THAT(print_wrap(read_result2), IsFailure(mio::StatusCode::OutOfRange)); + + // Case with empty RKI data file. + auto read_result3 = + mio::lsecir::get_initial_data_from_file(mio::path_join(TEST_DATA_DIR, "test_empty_file.json"), start_date, + infectionState, std::move(parameters), total_population, 1.); + + ASSERT_THAT(print_wrap(read_result3), IsFailure(mio::StatusCode::InvalidFileFormat)); + + // Reactive log output. + mio::set_log_level(mio::LogLevel::warn); +} \ No newline at end of file From a6a208a7ce2e3bbf8dd5020c19646dbb214aa948 Mon Sep 17 00:00:00 2001 From: lenaploetzke Date: Mon, 5 Feb 2024 14:28:49 +0100 Subject: [PATCH 04/10] Improved incorrect test file placement in cmake --- cpp/tests/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index e23087a860..47deb4195e 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -66,13 +66,13 @@ set(TESTSOURCES temp_file_register.h sanitizers.cpp test_lct_secir.cpp - test_lct_parameters_io.cpp ) if(MEMILIO_HAS_JSONCPP) set(TESTSOURCES ${TESTSOURCES} test_json_serializer.cpp test_epi_data_io.cpp +test_lct_parameters_io.cpp ) endif() From 88b7bb47ba304faaa7171d46c51633675db0fc27 Mon Sep 17 00:00:00 2001 From: lenaploetzke Date: Thu, 29 Feb 2024 10:01:27 +0100 Subject: [PATCH 05/10] corrected small error in initialization --- cpp/models/lct_secir/parameters_io.cpp | 2 +- cpp/models/lct_secir/parameters_io.h | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/cpp/models/lct_secir/parameters_io.cpp b/cpp/models/lct_secir/parameters_io.cpp index 2b3259d0e0..b82adcc6e8 100644 --- a/cpp/models/lct_secir/parameters_io.cpp +++ b/cpp/models/lct_secir/parameters_io.cpp @@ -75,7 +75,7 @@ IOResult get_initial_data_from_file(std::string const& path, Da // Go through the entries of rki_data and check if date is needed for calculation. Confirmed cases are scaled. for (auto&& entry : rki_data) { auto offset = get_offset_in_days(entry.date, date); - if (!(offset < min_offset_needed) || !(offset > max_offset_needed)) { + if (!(offset < min_offset_needed) && !(offset > max_offset_needed)) { // Add confirmed cases at date to compartment Recovered. if (offset == 0) { init[infectionState.get_firstindex(InfectionStateBase::Recovered)] += diff --git a/cpp/models/lct_secir/parameters_io.h b/cpp/models/lct_secir/parameters_io.h index b52b6cc20e..8c5477073b 100644 --- a/cpp/models/lct_secir/parameters_io.h +++ b/cpp/models/lct_secir/parameters_io.h @@ -47,13 +47,13 @@ namespace lsecir * Data could be downloaded eg with the file pycode/memilio-epidata/memilio/epidata/getCaseData.py, * which creates a file named cases_all_germany.json or a similar name. One should set impute_dates=True so that missing dates are imputed. * -* @param path Path to the RKI file. -* @param date Date for which the initial values should be computed. date is the start date of the simulation. -* @param infectionState InfectionState used to calculate the initial values. Defines eg the size of the calculated vector. -* @param parameters Parameters used to calculate the initial values. Defines eg the expected sojourntimes. -* @param total_population Total size of the population of the considered region. +* @param[in] path Path to the RKI file. +* @param[in] date Date for which the initial values should be computed. date is the start date of the simulation. +* @param[in] infectionState InfectionState used to calculate the initial values. Defines eg the size of the calculated vector. +* @param[in] parameters Parameters used to calculate the initial values. Defines eg the expected sojourntimes. +* @param[in] total_population Total size of the population of the considered region. * The sum of all values in the returned vector will be equal to that value. -* @param scale_confirmed_cases Factor by which to scale the confirmed cases of rki data to consider unreported cases. +* @param[in] scale_confirmed_cases Factor by which to scale the confirmed cases of rki data to consider unreported cases. * @returns Initialization Vector or any io errors that happen during reading of the files. */ IOResult get_initial_data_from_file(std::string const& path, Date date, InfectionState infectionState, From 017a334e54d10ea93fe999c8211c7d1a0defa104 Mon Sep 17 00:00:00 2001 From: lenaploetzke Date: Wed, 6 Mar 2024 16:23:27 +0100 Subject: [PATCH 06/10] adaption with template --- cpp/models/lct_secir/CMakeLists.txt | 1 - cpp/models/lct_secir/README.md | 4 +- cpp/models/lct_secir/model.h | 10 + cpp/models/lct_secir/parameters_io.cpp | 275 ------------------------- cpp/models/lct_secir/parameters_io.h | 274 ++++++++++++++++++++++-- cpp/tests/test_lct_parameters_io.cpp | 70 +++---- 6 files changed, 296 insertions(+), 338 deletions(-) delete mode 100644 cpp/models/lct_secir/parameters_io.cpp diff --git a/cpp/models/lct_secir/CMakeLists.txt b/cpp/models/lct_secir/CMakeLists.txt index 35a78f248d..362d09b1ee 100644 --- a/cpp/models/lct_secir/CMakeLists.txt +++ b/cpp/models/lct_secir/CMakeLists.txt @@ -5,7 +5,6 @@ add_library(lct_secir simulation.h parameters.h parameters_io.h - parameters_io.cpp ) target_link_libraries(lct_secir PUBLIC memilio) target_include_directories(lct_secir PUBLIC diff --git a/cpp/models/lct_secir/README.md b/cpp/models/lct_secir/README.md index b73a1f015d..8c3c963d7a 100644 --- a/cpp/models/lct_secir/README.md +++ b/cpp/models/lct_secir/README.md @@ -2,8 +2,8 @@ This model is based on the Linear Chain Trick. -The Linear Chain Trick provides the option to use Erlang-distributed sojourn times in the compartments through the use of subcompartments. -The normal ODE models have (possibly unrealistic) exponentially distributed sojourn times. +The Linear Chain Trick provides the option to use Erlang-distributed stay times in the compartments through the use of subcompartments. +The normal ODE models have (possibly unrealistic) exponentially distributed stay times. The LCT model can still be described by an ordinary differential equation system. For the concept see diff --git a/cpp/models/lct_secir/model.h b/cpp/models/lct_secir/model.h index 0c2f67a769..0d98c4dec1 100644 --- a/cpp/models/lct_secir/model.h +++ b/cpp/models/lct_secir/model.h @@ -283,6 +283,16 @@ class Model return m_initial_values; } + /** + * @brief Sets the initial values for the model. + * + * @param[in] init Vector with initial values for all infection states inclusive subcompartments. + */ + void set_initial_values(Eigen::VectorXd init) + { + m_initial_values = init; + } + Parameters parameters{}; ///< Parameters of the model. private: diff --git a/cpp/models/lct_secir/parameters_io.cpp b/cpp/models/lct_secir/parameters_io.cpp deleted file mode 100644 index b82adcc6e8..0000000000 --- a/cpp/models/lct_secir/parameters_io.cpp +++ /dev/null @@ -1,275 +0,0 @@ -/* -* 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 "lct_secir/parameters_io.h" -#include "memilio/config.h" - -#ifdef MEMILIO_HAS_JSONCPP - -#include "lct_secir/infection_state.h" -#include "lct_secir/parameters.h" -#include "memilio/io/epi_data.h" -#include "memilio/io/io.h" -#include "memilio/utils/date.h" -#include "memilio/utils/logging.h" -#include "memilio/math/eigen.h" - -#include - -namespace mio -{ -namespace lsecir -{ - -IOResult get_initial_data_from_file(std::string const& path, Date date, InfectionState infectionState, - Parameters&& parameters, ScalarType total_population, - ScalarType scale_confirmed_cases) -{ - // Try to get rki data from path. - BOOST_OUTCOME_TRY(rki_data, mio::read_confirmed_cases_noage(path)); - auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { - return a.date < b.date; - }); - if (max_date_entry == rki_data.end()) { - log_error("RKI data file is empty."); - return failure(StatusCode::InvalidFileFormat, path + ", file is empty."); - } - auto max_date = max_date_entry->date; - if (max_date < date) { - log_error("Specified date does not exist in RKI data."); - return failure(StatusCode::OutOfRange, path + ", specified date does not exist in RKI data."); - } - // Compute initial values for all subcompartments. - Eigen::VectorXd init = Eigen::VectorXd::Zero(infectionState.get_count()); - // Define variables for parameters that are often needed. - ScalarType timeExposed = parameters.get(); - ScalarType timeInfectedNoSymptoms = parameters.get(); - ScalarType timeInfectedSymptoms = parameters.get(); - ScalarType timeInfectedSevere = parameters.get(); - ScalarType timeInfectedCritical = parameters.get(); - ScalarType scale_mu_CR = 1 / (1 - parameters.get()); - - ScalarType min_offset_needed = std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical); - ScalarType max_offset_needed = std::ceil(timeExposed + timeInfectedNoSymptoms); - - bool min_offset_needed_avail = false; - bool max_offset_needed_avail = false; - - // Go through the entries of rki_data and check if date is needed for calculation. Confirmed cases are scaled. - for (auto&& entry : rki_data) { - auto offset = get_offset_in_days(entry.date, date); - if (!(offset < min_offset_needed) && !(offset > max_offset_needed)) { - // Add confirmed cases at date to compartment Recovered. - if (offset == 0) { - init[infectionState.get_firstindex(InfectionStateBase::Recovered)] += - scale_confirmed_cases * entry.num_confirmed; - } - - // Compute initial values for compartment InfectedNoSymptoms. - if (offset >= 0 && offset <= std::ceil(timeInfectedNoSymptoms)) { - ScalarType TCi = - timeInfectedNoSymptoms / infectionState.get_number(InfectionStateBase::InfectedNoSymptoms); - int idx_Cn = infectionState.get_firstindex(InfectionStateBase::InfectedNoSymptoms) + - infectionState.get_number(InfectionStateBase::InfectedNoSymptoms) - 1; - for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedNoSymptoms); i++) { - if (offset == std::floor(i * TCi)) { - init[idx_Cn - i] -= (1 - (i * TCi - std::floor(i * TCi))) * scale_mu_CR * - scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::ceil(i * TCi)) { - init[idx_Cn - i] -= - (i * TCi - std::floor(i * TCi)) * scale_mu_CR * scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::floor((i + 1) * TCi)) { - init[idx_Cn - i] += (1 - ((i + 1) * TCi - std::floor((i + 1) * TCi))) * scale_mu_CR * - scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::ceil((i + 1) * TCi)) { - init[idx_Cn - i] += ((i + 1) * TCi - std::floor((i + 1) * TCi)) * scale_mu_CR * - scale_confirmed_cases * entry.num_confirmed; - } - } - } - - // Compute initial values for compartment Exposed. - if (offset >= std::floor(timeInfectedNoSymptoms) && offset <= max_offset_needed) { - ScalarType TEi = timeExposed / infectionState.get_number(InfectionStateBase::Exposed); - int idx_En = infectionState.get_firstindex(InfectionStateBase::Exposed) + - infectionState.get_number(InfectionStateBase::Exposed) - 1; - for (int i = 0; i < infectionState.get_number(InfectionStateBase::Exposed); i++) { - if (offset == std::floor(timeInfectedNoSymptoms + i * TEi)) { - init[idx_En - i] -= - (1 - (timeInfectedNoSymptoms + i * TEi - std::floor(timeInfectedNoSymptoms + i * TEi))) * - scale_mu_CR * scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::ceil(timeInfectedNoSymptoms + i * TEi)) { - init[idx_En - i] -= - (timeInfectedNoSymptoms + i * TEi - std::floor(timeInfectedNoSymptoms + i * TEi)) * - scale_mu_CR * scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::floor(timeInfectedNoSymptoms + (i + 1) * TEi)) { - init[idx_En - i] += (1 - (timeInfectedNoSymptoms + (i + 1) * TEi - - std::floor(timeInfectedNoSymptoms + (i + 1) * TEi))) * - scale_mu_CR * scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::ceil(timeInfectedNoSymptoms + (i + 1) * TEi)) { - init[idx_En - i] += (timeInfectedNoSymptoms + (i + 1) * TEi - - std::floor(timeInfectedNoSymptoms + (i + 1) * TEi)) * - scale_mu_CR * scale_confirmed_cases * entry.num_confirmed; - } - } - } - - // Compute initial values for compartment InfectedSymptoms. - if (offset >= std::floor(-timeInfectedSymptoms) && offset <= 0) { - ScalarType TIi = timeInfectedSymptoms / infectionState.get_number(InfectionStateBase::InfectedSymptoms); - int idx_I1 = infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms); - for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedSymptoms); i++) { - if (offset == std::floor(-TIi * (i + 1))) { - init[idx_I1 + i] -= (1 - (-(i + 1) * TIi - std::floor(-(i + 1) * TIi))) * - scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::ceil(-TIi * (i + 1))) { - init[idx_I1 + i] -= - (-(i + 1) * TIi - std::floor(-(i + 1) * TIi)) * scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::floor(-TIi * i)) { - init[idx_I1 + i] += - (1 - (-i * TIi - std::floor(-i * TIi))) * scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::ceil(-TIi * i)) { - init[idx_I1 + i] += - (-i * TIi - std::floor(-i * TIi)) * scale_confirmed_cases * entry.num_confirmed; - } - } - } - - // Compute initial values for compartment InfectedSevere. - if (offset >= std::floor(-timeInfectedSymptoms - timeInfectedSevere) && - offset <= std::ceil(-timeInfectedSymptoms)) { - ScalarType THi = timeInfectedSevere / infectionState.get_number(InfectionStateBase::InfectedSevere); - ScalarType mu_IH = parameters.get(); - int idx_H1 = infectionState.get_firstindex(InfectionStateBase::InfectedSevere); - for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedSevere); i++) { - if (offset == std::floor(-timeInfectedSymptoms - THi * (i + 1))) { - init[idx_H1 + i] -= mu_IH * - (1 - (-timeInfectedSymptoms - (i + 1) * THi - - std::floor(-timeInfectedSymptoms - (i + 1) * THi))) * - scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::ceil(-timeInfectedSymptoms - THi * (i + 1))) { - init[idx_H1 + i] -= mu_IH * - (-timeInfectedSymptoms - (i + 1) * THi - - std::floor(-timeInfectedSymptoms - (i + 1) * THi)) * - scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::floor(-timeInfectedSymptoms - THi * i)) { - init[idx_H1 + i] += - mu_IH * - (1 - (-timeInfectedSymptoms - i * THi - std::floor(-timeInfectedSymptoms - i * THi))) * - scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::ceil(-timeInfectedSymptoms - THi * i)) { - init[idx_H1 + i] += - mu_IH * (-timeInfectedSymptoms - i * THi - std::floor(-timeInfectedSymptoms - i * THi)) * - scale_confirmed_cases * entry.num_confirmed; - } - } - } - - // Compute initial values for compartment InfectedCritical. - if (offset >= min_offset_needed && offset <= std::ceil(-timeInfectedSymptoms - timeInfectedSevere)) { - ScalarType TUi = timeInfectedCritical / infectionState.get_number(InfectionStateBase::InfectedCritical); - ScalarType mu_IH = parameters.get(); - ScalarType mu_HU = parameters.get(); - int idx_U1 = infectionState.get_firstindex(InfectionStateBase::InfectedCritical); - for (int i = 0; i < infectionState.get_number(InfectionStateBase::InfectedCritical); i++) { - if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere - TUi * (i + 1))) { - init[idx_U1 + i] -= - mu_IH * mu_HU * - (1 - (-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi - - std::floor(-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi))) * - scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - TUi * (i + 1))) { - init[idx_U1 + i] -= mu_IH * mu_HU * - (-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi - - std::floor(-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi)) * - scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere - TUi * i)) { - init[idx_U1 + i] += mu_IH * mu_HU * - (1 - (-timeInfectedSymptoms - timeInfectedSevere - i * TUi - - std::floor(-timeInfectedSymptoms - timeInfectedSevere - i * TUi))) * - scale_confirmed_cases * entry.num_confirmed; - } - if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - TUi * i)) { - init[idx_U1 + i] += mu_IH * mu_HU * - (-timeInfectedSymptoms - timeInfectedSevere - i * TUi - - std::floor(-timeInfectedSymptoms - timeInfectedSevere - i * TUi)) * - scale_confirmed_cases * entry.num_confirmed; - } - } - } - - // Compute Dead. - if (offset == min_offset_needed) { - min_offset_needed_avail = true; - init[infectionState.get_firstindex(InfectionStateBase::Dead)] += - (1 - (-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical - - std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical))) * - entry.num_deaths; - } - if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical)) { - init[infectionState.get_firstindex(InfectionStateBase::Dead)] += - (-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical - - std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical)) * - entry.num_deaths; - } - if (offset == max_offset_needed) { - max_offset_needed_avail = true; - } - } - } - - // Compute Recovered. - init[infectionState.get_firstindex(InfectionStateBase::Recovered)] -= - init.segment(infectionState.get_firstindex(InfectionStateBase::InfectedSymptoms), - infectionState.get_number(InfectionStateBase::InfectedSymptoms) + - infectionState.get_number(InfectionStateBase::InfectedSevere) + - infectionState.get_number(InfectionStateBase::InfectedCritical)) - .sum(); - init[infectionState.get_firstindex(InfectionStateBase::Recovered)] -= - init[infectionState.get_firstindex(InfectionStateBase::Dead)]; - - // Compute Susceptibles - init[infectionState.get_firstindex(InfectionStateBase::Susceptible)] = - total_population - init.segment(1, infectionState.get_count() - 1).sum(); - if (!max_offset_needed_avail || !min_offset_needed_avail) { - log_error("Necessary range of dates needed to compute initial values does not exist in RKI data."); - return failure(StatusCode::OutOfRange, path + ", necessary range of dates does not exist in RKI data."); - } - - return init; -} - -} // namespace lsecir -} // namespace mio -#endif // MEMILIO_HAS_JSONCPP \ No newline at end of file diff --git a/cpp/models/lct_secir/parameters_io.h b/cpp/models/lct_secir/parameters_io.h index 8c5477073b..1b049c1da3 100644 --- a/cpp/models/lct_secir/parameters_io.h +++ b/cpp/models/lct_secir/parameters_io.h @@ -25,43 +25,289 @@ #ifdef MEMILIO_HAS_JSONCPP -#include "lct_secir/parameters.h" +#include "lct_secir/model.h" #include "lct_secir/infection_state.h" +#include "lct_secir/parameters.h" +#include "memilio/io/epi_data.h" +#include "memilio/io/io.h" +#include "memilio/utils/date.h" +#include "memilio/utils/logging.h" #include "memilio/math/eigen.h" #include - +#include namespace mio { namespace lsecir { /** -* @brief Computes an initialization vector for a LCT model with data from RKI. +* @brief Computes an initialization vector for a LCT model with case data from RKI. * * Computes an initial value vector for an LCT model with the defined infectionState and the parameters. -* For the computation expected sojourntime in the subcompartments are used. To calculate the initial values, +* For the computation expected stay times in the subcompartments are used. To calculate the initial values, * we assume for simplicity that individuals stay in the subcompartment for exactly the expected time. * The RKI data are linearly interpolated within one day. * The RKI data should contain data for each needed day without division of age groups, the completeness of the dates is not verified. -* Data could be downloaded eg with the file pycode/memilio-epidata/memilio/epidata/getCaseData.py, +* Data can be downloaded e.g. with the file pycode/memilio-epidata/memilio/epidata/getCaseData.py, * which creates a file named cases_all_germany.json or a similar name. One should set impute_dates=True so that missing dates are imputed. * +* @param[in, out] model The model for which the inital data should be computed and set. * @param[in] path Path to the RKI file. * @param[in] date Date for which the initial values should be computed. date is the start date of the simulation. -* @param[in] infectionState InfectionState used to calculate the initial values. Defines eg the size of the calculated vector. -* @param[in] parameters Parameters used to calculate the initial values. Defines eg the expected sojourntimes. -* @param[in] total_population Total size of the population of the considered region. -* The sum of all values in the returned vector will be equal to that value. +* @param[in] total_population Total size of the population of Germany. +* The sum of all values in the vector of subcompartments will be equal to that value. * @param[in] scale_confirmed_cases Factor by which to scale the confirmed cases of rki data to consider unreported cases. -* @returns Initialization Vector or any io errors that happen during reading of the files. +* @returns Any io errors that happen during reading of the files. */ -IOResult get_initial_data_from_file(std::string const& path, Date date, InfectionState infectionState, - Parameters&& parameters, ScalarType total_population, - ScalarType scale_confirmed_cases = 1.); +template +IOResult set_initial_data_from_confirmed_cases(Model& model, const std::string& path, Date date, + ScalarType total_population, ScalarType scale_confirmed_cases = 1.) +{ + using LctState = typename Model::LctState; + + // Try to get rki data from path. + BOOST_OUTCOME_TRY(rki_data, mio::read_confirmed_cases_noage(path)); + auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) { + return a.date < b.date; + }); + if (max_date_entry == rki_data.end()) { + log_error("RKI data file is empty."); + return failure(StatusCode::InvalidFileFormat, path + ", file is empty."); + } + auto max_date = max_date_entry->date; + if (max_date < date) { + log_error("Specified date does not exist in RKI data."); + return failure(StatusCode::OutOfRange, path + ", specified date does not exist in RKI data."); + } + // Compute initial values for all subcompartments. + Eigen::VectorXd init = Eigen::VectorXd::Zero(LctState::Count); + // Define variables for parameters that are often needed. + Parameters parameters = model.parameters; + + ScalarType timeExposed = parameters.get(); + ScalarType timeInfectedNoSymptoms = parameters.get(); + ScalarType timeInfectedSymptoms = parameters.get(); + ScalarType timeInfectedSevere = parameters.get(); + ScalarType timeInfectedCritical = parameters.get(); + ScalarType scale_asymptomatic = 1 / (1 - parameters.get()); + + ScalarType min_offset_needed = std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical); + ScalarType max_offset_needed = std::ceil(timeExposed + timeInfectedNoSymptoms); + + bool min_offset_needed_avail = false; + bool max_offset_needed_avail = false; + + // Go through the entries of rki_data and check if date is needed for calculation. Confirmed cases are scaled by scale_confirmed_cases. + for (auto&& entry : rki_data) { + int offset = get_offset_in_days(entry.date, date); + if ((offset >= min_offset_needed) && (offset <= max_offset_needed)) { + // Add confirmed cases at date to compartment Recovered. + if (offset == 0) { + init[LctState::template get_first_index()] += + scale_confirmed_cases * entry.num_confirmed; + } + + // Compute initial values for compartment InfectedNoSymptoms. + if (offset >= 0 && offset <= std::ceil(timeInfectedNoSymptoms)) { + ScalarType TCi = timeInfectedNoSymptoms / + LctState::template get_num_subcompartments(); + int idx_Cn = LctState::template get_first_index() + + LctState::template get_num_subcompartments() - 1; + for (int i = 0; + i < (int)LctState::template get_num_subcompartments(); i++) { + if (offset == std::floor(i * TCi)) { + init[idx_Cn - i] -= (1 - (i * TCi - std::floor(i * TCi))) * scale_asymptomatic * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(i * TCi)) { + init[idx_Cn - i] -= (i * TCi - std::floor(i * TCi)) * scale_asymptomatic * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::floor((i + 1) * TCi)) { + init[idx_Cn - i] += (1 - ((i + 1) * TCi - std::floor((i + 1) * TCi))) * scale_asymptomatic * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil((i + 1) * TCi)) { + init[idx_Cn - i] += ((i + 1) * TCi - std::floor((i + 1) * TCi)) * scale_asymptomatic * + scale_confirmed_cases * entry.num_confirmed; + } + } + } + + // Compute initial values for compartment Exposed. + if (offset >= std::floor(timeInfectedNoSymptoms) && offset <= max_offset_needed) { + ScalarType TEi = timeExposed / LctState::template get_num_subcompartments(); + int idx_En = LctState::template get_first_index() + + LctState::template get_num_subcompartments() - 1; + for (int i = 0; i < (int)LctState::template get_num_subcompartments(); i++) { + if (offset == std::floor(timeInfectedNoSymptoms + i * TEi)) { + init[idx_En - i] -= + (1 - (timeInfectedNoSymptoms + i * TEi - std::floor(timeInfectedNoSymptoms + i * TEi))) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(timeInfectedNoSymptoms + i * TEi)) { + init[idx_En - i] -= + (timeInfectedNoSymptoms + i * TEi - std::floor(timeInfectedNoSymptoms + i * TEi)) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::floor(timeInfectedNoSymptoms + (i + 1) * TEi)) { + init[idx_En - i] += (1 - (timeInfectedNoSymptoms + (i + 1) * TEi - + std::floor(timeInfectedNoSymptoms + (i + 1) * TEi))) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(timeInfectedNoSymptoms + (i + 1) * TEi)) { + init[idx_En - i] += (timeInfectedNoSymptoms + (i + 1) * TEi - + std::floor(timeInfectedNoSymptoms + (i + 1) * TEi)) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; + } + } + } + + // Compute initial values for compartment InfectedSymptoms. + if (offset >= std::floor(-timeInfectedSymptoms) && offset <= 0) { + ScalarType TIi = + timeInfectedSymptoms / + (ScalarType)LctState::template get_num_subcompartments(); + int idx_I1 = LctState::template get_first_index(); + for (int i = 0; i < (int)LctState::template get_num_subcompartments(); + i++) { + if (offset == std::floor(-TIi * (i + 1))) { + init[idx_I1 + i] -= (1 - (-(i + 1.) * TIi - std::floor(-(i + 1) * TIi))) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-TIi * (i + 1))) { + init[idx_I1 + i] -= + (-(i + 1) * TIi - std::floor(-(i + 1) * TIi)) * scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::floor(-TIi * i)) { + init[idx_I1 + i] += (1 - (-1 * i * TIi - std::floor(-1 * i * TIi))) * scale_confirmed_cases * + entry.num_confirmed; + } + if (offset == std::ceil(-TIi * i)) { + init[idx_I1 + i] += + (-i * TIi - std::floor(-i * TIi)) * scale_confirmed_cases * entry.num_confirmed; + } + } + } + + // Compute initial values for compartment InfectedSevere. + if (offset >= std::floor(-timeInfectedSymptoms - timeInfectedSevere) && + offset <= std::ceil(-timeInfectedSymptoms)) { + ScalarType THi = + timeInfectedSevere / + (ScalarType)LctState::template get_num_subcompartments(); + ScalarType mu_IH = parameters.get(); + int idx_H1 = LctState::template get_first_index(); + for (int i = 0; i < (int)LctState::template get_num_subcompartments(); + i++) { + if (offset == std::floor(-timeInfectedSymptoms - THi * (i + 1))) { + init[idx_H1 + i] -= mu_IH * + (1 - (-timeInfectedSymptoms - (i + 1) * THi - + std::floor(-timeInfectedSymptoms - (i + 1) * THi))) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-timeInfectedSymptoms - THi * (i + 1))) { + init[idx_H1 + i] -= mu_IH * + (-timeInfectedSymptoms - (i + 1) * THi - + std::floor(-timeInfectedSymptoms - (i + 1) * THi)) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::floor(-timeInfectedSymptoms - THi * i)) { + init[idx_H1 + i] += + mu_IH * + (1 - (-timeInfectedSymptoms - i * THi - std::floor(-timeInfectedSymptoms - i * THi))) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-timeInfectedSymptoms - THi * i)) { + init[idx_H1 + i] += + mu_IH * (-timeInfectedSymptoms - i * THi - std::floor(-timeInfectedSymptoms - i * THi)) * + scale_confirmed_cases * entry.num_confirmed; + } + } + } + + // Compute initial values for compartment InfectedCritical. + if (offset >= min_offset_needed && offset <= std::ceil(-timeInfectedSymptoms - timeInfectedSevere)) { + ScalarType TUi = timeInfectedCritical / + LctState::template get_num_subcompartments(); + ScalarType mu_IH = parameters.get(); + ScalarType mu_HU = parameters.get(); + int idx_U1 = LctState::template get_first_index(); + for (int i = 0; i < (int)LctState::template get_num_subcompartments(); + i++) { + if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere - TUi * (i + 1))) { + init[idx_U1 + i] -= + mu_IH * mu_HU * + (1 - (-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi))) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - TUi * (i + 1))) { + init[idx_U1 + i] -= mu_IH * mu_HU * + (-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi)) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere - TUi * i)) { + init[idx_U1 + i] += mu_IH * mu_HU * + (1 - (-timeInfectedSymptoms - timeInfectedSevere - i * TUi - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - i * TUi))) * + scale_confirmed_cases * entry.num_confirmed; + } + if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - TUi * i)) { + init[idx_U1 + i] += mu_IH * mu_HU * + (-timeInfectedSymptoms - timeInfectedSevere - i * TUi - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - i * TUi)) * + scale_confirmed_cases * entry.num_confirmed; + } + } + } + + // Compute Dead. + if (offset == min_offset_needed) { + min_offset_needed_avail = true; + init[LctState::template get_first_index()] += + (1 - (-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical))) * + entry.num_deaths; + } + if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical)) { + init[LctState::template get_first_index()] += + (-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical)) * + entry.num_deaths; + } + if (offset == max_offset_needed) { + max_offset_needed_avail = true; + } + } + } + + // Compute Recovered. + init[LctState::template get_first_index()] -= + init.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments() + + LctState::template get_num_subcompartments() + + LctState::template get_num_subcompartments()) + .sum(); + init[LctState::template get_first_index()] -= + init[LctState::template get_first_index()]; + + // Compute Susceptibles + init[LctState::template get_first_index()] = + total_population - init.segment(1, LctState::Count - 1).sum(); + if (!max_offset_needed_avail || !min_offset_needed_avail) { + log_error("Necessary range of dates needed to compute initial values does not exist in RKI data."); + return failure(StatusCode::OutOfRange, path + ", necessary range of dates does not exist in RKI data."); + } + + model.set_initial_values(std::move(init)); + return mio::success(); +} } // namespace lsecir } // namespace mio #endif // MEMILIO_HAS_JSONCPP -#endif // LCTSECIR_PARAMETERS_IO_H \ No newline at end of file +#endif // LCTSECIR_PARAMETERS_IO_H diff --git a/cpp/tests/test_lct_parameters_io.cpp b/cpp/tests/test_lct_parameters_io.cpp index 615996e388..318ebe057e 100644 --- a/cpp/tests/test_lct_parameters_io.cpp +++ b/cpp/tests/test_lct_parameters_io.cpp @@ -50,35 +50,24 @@ TEST(TestLCTParametersIo, ReadPopulationDataRKI) parameters.get() = 0.3; parameters.get() = 0.2; - // Define number of subcompartments. - std::vector vec_subcompartments((int)mio::lsecir::InfectionStateBase::Count, 1); - // Use subcompartments with a soujourn time of approximately one day in each subcompartment. - vec_subcompartments[(int)mio::lsecir::InfectionStateBase::Exposed] = - (int)round(parameters.get()); - vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedNoSymptoms] = - (int)round(parameters.get()); - vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedSymptoms] = - (int)round(parameters.get()); - vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedSevere] = - (int)round(parameters.get()); - // Both realistic distributions for times corresponding to InfectedCritical of the IDE model are exponential distributions. - vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedCritical] = 1; - mio::lsecir::InfectionState infectionState(vec_subcompartments); + // Use rounded stay times for the number of subcompartments except for InfectedCritical. + // Template parameters must be compile-time constants, which is why the values are not calculated with parameters. + using Model = mio::lsecir::Model<2, 3, 2, 2, 1>; + using LctState = Model::LctState; + Model model(Eigen::VectorXd::Zero(LctState::Count), std::move(parameters)); // Calculate initial value vector for subcompartments with RKI data. - auto read_result = - mio::lsecir::get_initial_data_from_file(mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date, - infectionState, std::move(parameters), total_population, 1.); + auto read_result = mio::lsecir::set_initial_data_from_confirmed_cases( + model, mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date, total_population, 1.); ASSERT_THAT(print_wrap(read_result), IsSuccess()); - auto init_subcompartments = read_result.value(); // Previous result. - Eigen::VectorXd compare(infectionState.get_count()); + Eigen::VectorXd compare(LctState::Count); compare << 863.05, 14.30625, 8.53125, 30.1125, 36.1875, 3.8125, 9.88, 3.52, 0.09, 0.25, 0.6888, 27.8712, 1.7; - for (int i = 0; i < infectionState.get_count(); i++) { - EXPECT_NEAR(init_subcompartments[i], compare[i], 1e-4) << "at subcompartment number " << i; + for (unsigned int i = 0; i < LctState::Count; i++) { + EXPECT_NEAR(model.get_initial_values()[i], compare[i], 1e-4) << "at subcompartment number " << i; } } @@ -100,43 +89,32 @@ TEST(TestLCTParametersIo, ReadPopulationDataRKIFailure) parameters.get() = 0.15; parameters.get() = 0.22; - // Define number of subcompartments. - std::vector vec_subcompartments((int)mio::lsecir::InfectionStateBase::Count, 1); - // Use subcompartments with a soujourn time of approximately one day in each subcompartment. - vec_subcompartments[(int)mio::lsecir::InfectionStateBase::Exposed] = - (int)round(parameters.get()); - vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedNoSymptoms] = - (int)round(parameters.get()); - vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedSymptoms] = - (int)round(parameters.get()); - vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedSevere] = - (int)round(parameters.get()); - // Both realistic distributions for times corresponding to InfectedCritical of the IDE model are exponential distributions. - vec_subcompartments[(int)mio::lsecir::InfectionStateBase::InfectedCritical] = 1; - mio::lsecir::InfectionState infectionState(vec_subcompartments); + /// Use rounded stay times for the number of subcompartments except for InfectedCritical. + // Template parameters must be compile-time constants, which is why the values are not calculated with parameters. + using Model = mio::lsecir::Model<2, 3, 2, 2, 1>; + using LctState = Model::LctState; + Model model(Eigen::VectorXd::Zero(LctState::Count), std::move(parameters)); // Deactivate temporarily log output for next tests. mio::set_log_level(mio::LogLevel::off); // Case where start_date is later than maximal provided date in file. - auto start_date = mio::Date(2020, 6, 9); - auto read_result1 = - mio::lsecir::get_initial_data_from_file(mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date, - infectionState, std::move(parameters), total_population, 1.); + auto start_date = mio::Date(2020, 6, 9); + auto read_result1 = mio::lsecir::set_initial_data_from_confirmed_cases( + model, mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date, total_population, 1.); ASSERT_THAT(print_wrap(read_result1), IsFailure(mio::StatusCode::OutOfRange)); + // Case where not all needed dates are provided. - start_date = mio::Date(2020, 6, 6); - auto read_result2 = - mio::lsecir::get_initial_data_from_file(mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date, - infectionState, std::move(parameters), total_population, 1.); + start_date = mio::Date(2020, 6, 6); + auto read_result2 = mio::lsecir::set_initial_data_from_confirmed_cases( + model, mio::path_join(TEST_DATA_DIR, "cases_all_germany.json"), start_date, total_population, 1.); ASSERT_THAT(print_wrap(read_result2), IsFailure(mio::StatusCode::OutOfRange)); // Case with empty RKI data file. - auto read_result3 = - mio::lsecir::get_initial_data_from_file(mio::path_join(TEST_DATA_DIR, "test_empty_file.json"), start_date, - infectionState, std::move(parameters), total_population, 1.); + auto read_result3 = mio::lsecir::set_initial_data_from_confirmed_cases( + model, mio::path_join(TEST_DATA_DIR, "test_empty_file.json"), start_date, total_population, 1.); ASSERT_THAT(print_wrap(read_result3), IsFailure(mio::StatusCode::InvalidFileFormat)); From 4f7c2f78900949e48d4778eb54638bb2d7320335 Mon Sep 17 00:00:00 2001 From: lenaploetzke Date: Thu, 7 Mar 2024 11:41:48 +0100 Subject: [PATCH 07/10] included review comments Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com> --- cpp/memilio/io/epi_data.h | 16 +++++++------- cpp/models/lct_secir/parameters_io.h | 32 ++++++++++++++++------------ 2 files changed, 26 insertions(+), 22 deletions(-) diff --git a/cpp/memilio/io/epi_data.h b/cpp/memilio/io/epi_data.h index fd1c4d024d..412a27ab4d 100644 --- a/cpp/memilio/io/epi_data.h +++ b/cpp/memilio/io/epi_data.h @@ -67,7 +67,7 @@ class StringDate : public Date /** * Represents the entries of a confirmed cases data file, e.g., from RKI. - * Number of confirmed, recovered and deceased in a Germany on a specific date. + * Number of confirmed, recovered and deceased in Germany on a specific date. * ConfirmedCasesNoAgeEntry is a simplified version of ConfirmedCasesDataEntry without agegroups. */ class ConfirmedCasesNoAgeEntry @@ -97,9 +97,9 @@ class ConfirmedCasesNoAgeEntry /** * Deserialize a list of ConfirmedCasesNoAgeEntry from json. - * @param jsvalue json value, must be an array of objects, objects must match ConfirmedCasesNoAgeEntry. - * Exactly one entry per date should be provided, for example no entries per age group. - * @return list of ConfirmedCasesNoAgeEntry. + * @param jsvalue json value, must be an array of objects, objects must match ConfirmedCasesNoAgeEntry. + * Accordingly, there should be one entry per date with the values for Confirmed, Recovered and Deaths. In addition, no age groups should be specified. + * @return List of ConfirmedCasesNoAgeEntry. */ inline IOResult> deserialize_confirmed_cases_noage(const Json::Value& jsvalue) { @@ -107,10 +107,10 @@ inline IOResult> deserialize_confirmed_cas } /** - * Deserialize a list of ConfirmedCasesNoAgeEntry from json. - * @param jsvalue json value, must be an array of objects, objects must match ConfirmedCasesNoAgeEntry. - * Exactly one entry per date should be provided, for example no entries per age group - * @return list of ConfirmedCasesNoAgeEntry. + * Read list of ConfirmedCasesNoAgeEntry from a json file. + * @param filename name of the json file. File content must be an array of objects, objects must match ConfirmedCasesNoAgeEntry. + * Accordingly, there should be one entry per date with the values for Confirmed, Recovered and Deaths. In addition, no age groups should be specified. + * @return List of ConfirmedCasesNoAgeEntry. */ inline IOResult> read_confirmed_cases_noage(const std::string& filename) { diff --git a/cpp/models/lct_secir/parameters_io.h b/cpp/models/lct_secir/parameters_io.h index 1b049c1da3..bebb2f2b87 100644 --- a/cpp/models/lct_secir/parameters_io.h +++ b/cpp/models/lct_secir/parameters_io.h @@ -44,10 +44,10 @@ namespace lsecir /** * @brief Computes an initialization vector for a LCT model with case data from RKI. * -* Computes an initial value vector for an LCT model with the defined infectionState and the parameters. +* Calculates an initial value vector for an LCT model and updates the initial value vector in the model. * For the computation expected stay times in the subcompartments are used. To calculate the initial values, * we assume for simplicity that individuals stay in the subcompartment for exactly the expected time. -* The RKI data are linearly interpolated within one day. +* The RKI data are linearly interpolated within one day to match the expected stay time in a subcompartment. * The RKI data should contain data for each needed day without division of age groups, the completeness of the dates is not verified. * Data can be downloaded e.g. with the file pycode/memilio-epidata/memilio/epidata/getCaseData.py, * which creates a file named cases_all_germany.json or a similar name. One should set impute_dates=True so that missing dates are imputed. @@ -58,6 +58,7 @@ namespace lsecir * @param[in] total_population Total size of the population of Germany. * The sum of all values in the vector of subcompartments will be equal to that value. * @param[in] scale_confirmed_cases Factor by which to scale the confirmed cases of rki data to consider unreported cases. +* @tparam Model is expected to be an LCT-SECIR model defined in models/lct_secir/model.h. * @returns Any io errors that happen during reading of the files. */ template @@ -83,14 +84,13 @@ IOResult set_initial_data_from_confirmed_cases(Model& model, const std::st // Compute initial values for all subcompartments. Eigen::VectorXd init = Eigen::VectorXd::Zero(LctState::Count); // Define variables for parameters that are often needed. - Parameters parameters = model.parameters; - ScalarType timeExposed = parameters.get(); - ScalarType timeInfectedNoSymptoms = parameters.get(); - ScalarType timeInfectedSymptoms = parameters.get(); - ScalarType timeInfectedSevere = parameters.get(); - ScalarType timeInfectedCritical = parameters.get(); - ScalarType scale_asymptomatic = 1 / (1 - parameters.get()); + ScalarType timeExposed = model.parameters.template get(); + ScalarType timeInfectedNoSymptoms = model.parameters.template get(); + ScalarType timeInfectedSymptoms = model.parameters.template get(); + ScalarType timeInfectedSevere = model.parameters.template get(); + ScalarType timeInfectedCritical = model.parameters.template get(); + ScalarType scale_asymptomatic = 1 / (1 - model.parameters.template get()); ScalarType min_offset_needed = std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical); ScalarType max_offset_needed = std::ceil(timeExposed + timeInfectedNoSymptoms); @@ -181,8 +181,8 @@ IOResult set_initial_data_from_confirmed_cases(Model& model, const std::st (-(i + 1) * TIi - std::floor(-(i + 1) * TIi)) * scale_confirmed_cases * entry.num_confirmed; } if (offset == std::floor(-TIi * i)) { - init[idx_I1 + i] += (1 - (-1 * i * TIi - std::floor(-1 * i * TIi))) * scale_confirmed_cases * - entry.num_confirmed; + init[idx_I1 + i] += + (1 - (-i * TIi - std::floor(-i * TIi))) * scale_confirmed_cases * entry.num_confirmed; } if (offset == std::ceil(-TIi * i)) { init[idx_I1 + i] += @@ -197,7 +197,7 @@ IOResult set_initial_data_from_confirmed_cases(Model& model, const std::st ScalarType THi = timeInfectedSevere / (ScalarType)LctState::template get_num_subcompartments(); - ScalarType mu_IH = parameters.get(); + ScalarType mu_IH = model.parameters.template get(); int idx_H1 = LctState::template get_first_index(); for (int i = 0; i < (int)LctState::template get_num_subcompartments(); i++) { @@ -231,8 +231,8 @@ IOResult set_initial_data_from_confirmed_cases(Model& model, const std::st if (offset >= min_offset_needed && offset <= std::ceil(-timeInfectedSymptoms - timeInfectedSevere)) { ScalarType TUi = timeInfectedCritical / LctState::template get_num_subcompartments(); - ScalarType mu_IH = parameters.get(); - ScalarType mu_HU = parameters.get(); + ScalarType mu_IH = model.parameters.template get(); + ScalarType mu_HU = model.parameters.template get(); int idx_U1 = LctState::template get_first_index(); for (int i = 0; i < (int)LctState::template get_num_subcompartments(); i++) { @@ -297,12 +297,16 @@ IOResult set_initial_data_from_confirmed_cases(Model& model, const std::st // Compute Susceptibles init[LctState::template get_first_index()] = total_population - init.segment(1, LctState::Count - 1).sum(); + + // Check whether all necessary dates are available. if (!max_offset_needed_avail || !min_offset_needed_avail) { log_error("Necessary range of dates needed to compute initial values does not exist in RKI data."); return failure(StatusCode::OutOfRange, path + ", necessary range of dates does not exist in RKI data."); } + // Set computed initial value vector model.set_initial_values(std::move(init)); + return mio::success(); } From 4cae29c1e1c62218a12d713998bfd3090ceeac4f Mon Sep 17 00:00:00 2001 From: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> Date: Thu, 7 Mar 2024 12:03:23 +0100 Subject: [PATCH 08/10] Try if this co-authorship works. Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com> --- cpp/tests/test_lct_parameters_io.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/tests/test_lct_parameters_io.cpp b/cpp/tests/test_lct_parameters_io.cpp index 318ebe057e..87ea792d47 100644 --- a/cpp/tests/test_lct_parameters_io.cpp +++ b/cpp/tests/test_lct_parameters_io.cpp @@ -34,6 +34,7 @@ // Check that Initialization based on synthetic RKI data match previous result. TEST(TestLCTParametersIo, ReadPopulationDataRKI) { + // Define start date and the total population used for the initialization. ScalarType total_population = 1000.0; auto start_date = mio::Date(2020, 6, 1); From a3a6aed7065a481e4abbeee4326f52b2aefe82f3 Mon Sep 17 00:00:00 2001 From: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> Date: Fri, 8 Mar 2024 10:42:08 +0100 Subject: [PATCH 09/10] Update cpp/models/lct_secir/model.h --- cpp/models/lct_secir/model.h | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/models/lct_secir/model.h b/cpp/models/lct_secir/model.h index 0d98c4dec1..34c13da8af 100644 --- a/cpp/models/lct_secir/model.h +++ b/cpp/models/lct_secir/model.h @@ -291,6 +291,7 @@ class Model void set_initial_values(Eigen::VectorXd init) { m_initial_values = init; + m_N0 = m_initial_values.sum(); } Parameters parameters{}; ///< Parameters of the model. From f96f76da34c2a86407f0a533b062f4919a00b87c Mon Sep 17 00:00:00 2001 From: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> Date: Mon, 18 Mar 2024 10:47:16 +0100 Subject: [PATCH 10/10] renamed compartments --- cpp/models/lct_secir/parameters_io.h | 231 ++++++++++++++++----------- cpp/tests/test_lct_parameters_io.cpp | 4 +- 2 files changed, 137 insertions(+), 98 deletions(-) diff --git a/cpp/models/lct_secir/parameters_io.h b/cpp/models/lct_secir/parameters_io.h index bebb2f2b87..776023dcbe 100644 --- a/cpp/models/lct_secir/parameters_io.h +++ b/cpp/models/lct_secir/parameters_io.h @@ -110,83 +110,102 @@ IOResult set_initial_data_from_confirmed_cases(Model& model, const std::st // Compute initial values for compartment InfectedNoSymptoms. if (offset >= 0 && offset <= std::ceil(timeInfectedNoSymptoms)) { - ScalarType TCi = timeInfectedNoSymptoms / - LctState::template get_num_subcompartments(); - int idx_Cn = LctState::template get_first_index() + - LctState::template get_num_subcompartments() - 1; + // Mean stay time in each subcompartment. + ScalarType timeInfectedNoSymptoms_i = + timeInfectedNoSymptoms / + LctState::template get_num_subcompartments(); + // Index of the last subcompartment of InfectedNoSymptoms. + int idxInfectedNoSymptoms_last = + LctState::template get_first_index() + + LctState::template get_num_subcompartments() - 1; for (int i = 0; i < (int)LctState::template get_num_subcompartments(); i++) { - if (offset == std::floor(i * TCi)) { - init[idx_Cn - i] -= (1 - (i * TCi - std::floor(i * TCi))) * scale_asymptomatic * - scale_confirmed_cases * entry.num_confirmed; + if (offset == std::floor(i * timeInfectedNoSymptoms_i)) { + init[idxInfectedNoSymptoms_last - i] -= + (1 - (i * timeInfectedNoSymptoms_i - std::floor(i * timeInfectedNoSymptoms_i))) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::ceil(i * TCi)) { - init[idx_Cn - i] -= (i * TCi - std::floor(i * TCi)) * scale_asymptomatic * - scale_confirmed_cases * entry.num_confirmed; + if (offset == std::ceil(i * timeInfectedNoSymptoms_i)) { + init[idxInfectedNoSymptoms_last - i] -= + (i * timeInfectedNoSymptoms_i - std::floor(i * timeInfectedNoSymptoms_i)) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::floor((i + 1) * TCi)) { - init[idx_Cn - i] += (1 - ((i + 1) * TCi - std::floor((i + 1) * TCi))) * scale_asymptomatic * - scale_confirmed_cases * entry.num_confirmed; + if (offset == std::floor((i + 1) * timeInfectedNoSymptoms_i)) { + init[idxInfectedNoSymptoms_last - i] += + (1 - + ((i + 1) * timeInfectedNoSymptoms_i - std::floor((i + 1) * timeInfectedNoSymptoms_i))) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::ceil((i + 1) * TCi)) { - init[idx_Cn - i] += ((i + 1) * TCi - std::floor((i + 1) * TCi)) * scale_asymptomatic * - scale_confirmed_cases * entry.num_confirmed; + if (offset == std::ceil((i + 1) * timeInfectedNoSymptoms_i)) { + init[idxInfectedNoSymptoms_last - i] += + ((i + 1) * timeInfectedNoSymptoms_i - std::floor((i + 1) * timeInfectedNoSymptoms_i)) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; } } } // Compute initial values for compartment Exposed. if (offset >= std::floor(timeInfectedNoSymptoms) && offset <= max_offset_needed) { - ScalarType TEi = timeExposed / LctState::template get_num_subcompartments(); - int idx_En = LctState::template get_first_index() + - LctState::template get_num_subcompartments() - 1; + // Mean stay time in each subcompartment. + ScalarType timeExposed_i = + timeExposed / LctState::template get_num_subcompartments(); + // Index of the last subcompartment of Exposed. + int idxExposed_last = LctState::template get_first_index() + + LctState::template get_num_subcompartments() - 1; for (int i = 0; i < (int)LctState::template get_num_subcompartments(); i++) { - if (offset == std::floor(timeInfectedNoSymptoms + i * TEi)) { - init[idx_En - i] -= - (1 - (timeInfectedNoSymptoms + i * TEi - std::floor(timeInfectedNoSymptoms + i * TEi))) * - scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; + if (offset == std::floor(timeInfectedNoSymptoms + i * timeExposed_i)) { + init[idxExposed_last - i] -= (1 - (timeInfectedNoSymptoms + i * timeExposed_i - + std::floor(timeInfectedNoSymptoms + i * timeExposed_i))) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::ceil(timeInfectedNoSymptoms + i * TEi)) { - init[idx_En - i] -= - (timeInfectedNoSymptoms + i * TEi - std::floor(timeInfectedNoSymptoms + i * TEi)) * - scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; + if (offset == std::ceil(timeInfectedNoSymptoms + i * timeExposed_i)) { + init[idxExposed_last - i] -= (timeInfectedNoSymptoms + i * timeExposed_i - + std::floor(timeInfectedNoSymptoms + i * timeExposed_i)) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::floor(timeInfectedNoSymptoms + (i + 1) * TEi)) { - init[idx_En - i] += (1 - (timeInfectedNoSymptoms + (i + 1) * TEi - - std::floor(timeInfectedNoSymptoms + (i + 1) * TEi))) * - scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; + if (offset == std::floor(timeInfectedNoSymptoms + (i + 1) * timeExposed_i)) { + init[idxExposed_last - i] += + (1 - (timeInfectedNoSymptoms + (i + 1) * timeExposed_i - + std::floor(timeInfectedNoSymptoms + (i + 1) * timeExposed_i))) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::ceil(timeInfectedNoSymptoms + (i + 1) * TEi)) { - init[idx_En - i] += (timeInfectedNoSymptoms + (i + 1) * TEi - - std::floor(timeInfectedNoSymptoms + (i + 1) * TEi)) * - scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; + if (offset == std::ceil(timeInfectedNoSymptoms + (i + 1) * timeExposed_i)) { + init[idxExposed_last - i] += (timeInfectedNoSymptoms + (i + 1) * timeExposed_i - + std::floor(timeInfectedNoSymptoms + (i + 1) * timeExposed_i)) * + scale_asymptomatic * scale_confirmed_cases * entry.num_confirmed; } } } // Compute initial values for compartment InfectedSymptoms. if (offset >= std::floor(-timeInfectedSymptoms) && offset <= 0) { - ScalarType TIi = + // Mean stay time in each subcompartment. + ScalarType timeInfectedSymptoms_i = timeInfectedSymptoms / (ScalarType)LctState::template get_num_subcompartments(); - int idx_I1 = LctState::template get_first_index(); + // Index of the first subcompartment of InfectedSymptoms. + int idxInfectedSymptoms_first = LctState::template get_first_index(); for (int i = 0; i < (int)LctState::template get_num_subcompartments(); i++) { - if (offset == std::floor(-TIi * (i + 1))) { - init[idx_I1 + i] -= (1 - (-(i + 1.) * TIi - std::floor(-(i + 1) * TIi))) * - scale_confirmed_cases * entry.num_confirmed; + if (offset == std::floor(-timeInfectedSymptoms_i * (i + 1))) { + init[idxInfectedSymptoms_first + i] -= + (1 - (-(i + 1.) * timeInfectedSymptoms_i - std::floor(-(i + 1) * timeInfectedSymptoms_i))) * + scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::ceil(-TIi * (i + 1))) { - init[idx_I1 + i] -= - (-(i + 1) * TIi - std::floor(-(i + 1) * TIi)) * scale_confirmed_cases * entry.num_confirmed; + if (offset == std::ceil(-timeInfectedSymptoms_i * (i + 1))) { + init[idxInfectedSymptoms_first + i] -= + (-(i + 1) * timeInfectedSymptoms_i - std::floor(-(i + 1) * timeInfectedSymptoms_i)) * + scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::floor(-TIi * i)) { - init[idx_I1 + i] += - (1 - (-i * TIi - std::floor(-i * TIi))) * scale_confirmed_cases * entry.num_confirmed; + if (offset == std::floor(-timeInfectedSymptoms_i * i)) { + init[idxInfectedSymptoms_first + i] += + (1 - (-i * timeInfectedSymptoms_i - std::floor(-i * timeInfectedSymptoms_i))) * + scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::ceil(-TIi * i)) { - init[idx_I1 + i] += - (-i * TIi - std::floor(-i * TIi)) * scale_confirmed_cases * entry.num_confirmed; + if (offset == std::ceil(-timeInfectedSymptoms_i * i)) { + init[idxInfectedSymptoms_first + i] += + (-i * timeInfectedSymptoms_i - std::floor(-i * timeInfectedSymptoms_i)) * + scale_confirmed_cases * entry.num_confirmed; } } } @@ -194,34 +213,42 @@ IOResult set_initial_data_from_confirmed_cases(Model& model, const std::st // Compute initial values for compartment InfectedSevere. if (offset >= std::floor(-timeInfectedSymptoms - timeInfectedSevere) && offset <= std::ceil(-timeInfectedSymptoms)) { - ScalarType THi = + // Mean stay time in each subcompartment. + ScalarType timeInfectedSevere_i = timeInfectedSevere / (ScalarType)LctState::template get_num_subcompartments(); - ScalarType mu_IH = model.parameters.template get(); - int idx_H1 = LctState::template get_first_index(); + // Transmission probability. + ScalarType prob_SeverePerInfectedSymptoms = model.parameters.template get(); + // Index of the first subcompartment of InfectedSevere. + int idxInfectedSevere_first = LctState::template get_first_index(); for (int i = 0; i < (int)LctState::template get_num_subcompartments(); i++) { - if (offset == std::floor(-timeInfectedSymptoms - THi * (i + 1))) { - init[idx_H1 + i] -= mu_IH * - (1 - (-timeInfectedSymptoms - (i + 1) * THi - - std::floor(-timeInfectedSymptoms - (i + 1) * THi))) * - scale_confirmed_cases * entry.num_confirmed; + if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere_i * (i + 1))) { + init[idxInfectedSevere_first + i] -= + prob_SeverePerInfectedSymptoms * + (1 - (-timeInfectedSymptoms - (i + 1) * timeInfectedSevere_i - + std::floor(-timeInfectedSymptoms - (i + 1) * timeInfectedSevere_i))) * + scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::ceil(-timeInfectedSymptoms - THi * (i + 1))) { - init[idx_H1 + i] -= mu_IH * - (-timeInfectedSymptoms - (i + 1) * THi - - std::floor(-timeInfectedSymptoms - (i + 1) * THi)) * - scale_confirmed_cases * entry.num_confirmed; + if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere_i * (i + 1))) { + init[idxInfectedSevere_first + i] -= + prob_SeverePerInfectedSymptoms * + (-timeInfectedSymptoms - (i + 1) * timeInfectedSevere_i - + std::floor(-timeInfectedSymptoms - (i + 1) * timeInfectedSevere_i)) * + scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::floor(-timeInfectedSymptoms - THi * i)) { - init[idx_H1 + i] += - mu_IH * - (1 - (-timeInfectedSymptoms - i * THi - std::floor(-timeInfectedSymptoms - i * THi))) * + if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere_i * i)) { + init[idxInfectedSevere_first + i] += + prob_SeverePerInfectedSymptoms * + (1 - (-timeInfectedSymptoms - i * timeInfectedSevere_i - + std::floor(-timeInfectedSymptoms - i * timeInfectedSevere_i))) * scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::ceil(-timeInfectedSymptoms - THi * i)) { - init[idx_H1 + i] += - mu_IH * (-timeInfectedSymptoms - i * THi - std::floor(-timeInfectedSymptoms - i * THi)) * + if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere_i * i)) { + init[idxInfectedSevere_first + i] += + prob_SeverePerInfectedSymptoms * + (-timeInfectedSymptoms - i * timeInfectedSevere_i - + std::floor(-timeInfectedSymptoms - i * timeInfectedSevere_i)) * scale_confirmed_cases * entry.num_confirmed; } } @@ -229,37 +256,49 @@ IOResult set_initial_data_from_confirmed_cases(Model& model, const std::st // Compute initial values for compartment InfectedCritical. if (offset >= min_offset_needed && offset <= std::ceil(-timeInfectedSymptoms - timeInfectedSevere)) { - ScalarType TUi = timeInfectedCritical / - LctState::template get_num_subcompartments(); - ScalarType mu_IH = model.parameters.template get(); - ScalarType mu_HU = model.parameters.template get(); - int idx_U1 = LctState::template get_first_index(); + // Mean stay time in each subcompartment. + ScalarType timeInfectedCritical_i = + timeInfectedCritical / + LctState::template get_num_subcompartments(); + // Transmission probabilities. + ScalarType prob_SeverePerInfectedSymptoms = model.parameters.template get(); + ScalarType prob_CriticalPerSevere = model.parameters.template get(); + // Index of the first subcompartment of InfectedCritical. + int idxInfectedCritical_first = LctState::template get_first_index(); for (int i = 0; i < (int)LctState::template get_num_subcompartments(); i++) { - if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere - TUi * (i + 1))) { - init[idx_U1 + i] -= - mu_IH * mu_HU * - (1 - (-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi - - std::floor(-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi))) * + if (offset == + std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical_i * (i + 1))) { + init[idxInfectedCritical_first + i] -= + prob_SeverePerInfectedSymptoms * prob_CriticalPerSevere * + (1 - (-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * timeInfectedCritical_i - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - + (i + 1) * timeInfectedCritical_i))) * scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - TUi * (i + 1))) { - init[idx_U1 + i] -= mu_IH * mu_HU * - (-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi - - std::floor(-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * TUi)) * - scale_confirmed_cases * entry.num_confirmed; + if (offset == + std::ceil(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical_i * (i + 1))) { + init[idxInfectedCritical_first + i] -= + prob_SeverePerInfectedSymptoms * prob_CriticalPerSevere * + (-timeInfectedSymptoms - timeInfectedSevere - (i + 1) * timeInfectedCritical_i - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - + (i + 1) * timeInfectedCritical_i)) * + scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere - TUi * i)) { - init[idx_U1 + i] += mu_IH * mu_HU * - (1 - (-timeInfectedSymptoms - timeInfectedSevere - i * TUi - - std::floor(-timeInfectedSymptoms - timeInfectedSevere - i * TUi))) * - scale_confirmed_cases * entry.num_confirmed; + if (offset == std::floor(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical_i * i)) { + init[idxInfectedCritical_first + i] += + prob_SeverePerInfectedSymptoms * prob_CriticalPerSevere * + (1 - + (-timeInfectedSymptoms - timeInfectedSevere - i * timeInfectedCritical_i - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - i * timeInfectedCritical_i))) * + scale_confirmed_cases * entry.num_confirmed; } - if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - TUi * i)) { - init[idx_U1 + i] += mu_IH * mu_HU * - (-timeInfectedSymptoms - timeInfectedSevere - i * TUi - - std::floor(-timeInfectedSymptoms - timeInfectedSevere - i * TUi)) * - scale_confirmed_cases * entry.num_confirmed; + if (offset == std::ceil(-timeInfectedSymptoms - timeInfectedSevere - timeInfectedCritical_i * i)) { + init[idxInfectedCritical_first + i] += + prob_SeverePerInfectedSymptoms * prob_CriticalPerSevere * + (-timeInfectedSymptoms - timeInfectedSevere - i * timeInfectedCritical_i - + std::floor(-timeInfectedSymptoms - timeInfectedSevere - i * timeInfectedCritical_i)) * + scale_confirmed_cases * entry.num_confirmed; } } } @@ -294,7 +333,7 @@ IOResult set_initial_data_from_confirmed_cases(Model& model, const std::st init[LctState::template get_first_index()] -= init[LctState::template get_first_index()]; - // Compute Susceptibles + // Compute Susceptibles. init[LctState::template get_first_index()] = total_population - init.segment(1, LctState::Count - 1).sum(); @@ -304,7 +343,7 @@ IOResult set_initial_data_from_confirmed_cases(Model& model, const std::st return failure(StatusCode::OutOfRange, path + ", necessary range of dates does not exist in RKI data."); } - // Set computed initial value vector + // Set computed initial value vector. model.set_initial_values(std::move(init)); return mio::success(); diff --git a/cpp/tests/test_lct_parameters_io.cpp b/cpp/tests/test_lct_parameters_io.cpp index 87ea792d47..54ec501988 100644 --- a/cpp/tests/test_lct_parameters_io.cpp +++ b/cpp/tests/test_lct_parameters_io.cpp @@ -31,7 +31,7 @@ #include -// Check that Initialization based on synthetic RKI data match previous result. +// Check that initialization based on synthetic RKI data match previous result. TEST(TestLCTParametersIo, ReadPopulationDataRKI) { // Define start date and the total population used for the initialization. @@ -63,7 +63,7 @@ TEST(TestLCTParametersIo, ReadPopulationDataRKI) ASSERT_THAT(print_wrap(read_result), IsSuccess()); - // Previous result. + // Result of a previous simulation. Eigen::VectorXd compare(LctState::Count); compare << 863.05, 14.30625, 8.53125, 30.1125, 36.1875, 3.8125, 9.88, 3.52, 0.09, 0.25, 0.6888, 27.8712, 1.7;