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..7231356c33 --- /dev/null +++ b/cpp/examples/lct_secir.cpp @@ -0,0 +1,92 @@ +/* +* 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 run a simulation using an LCT SECIR model. + Parameters, initial values and subcompartments are not meant to represent a realistic scenario. */ + + // 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 infection_state(num_subcompartments); + + ScalarType tmax = 20; + + // Define initial distribution of the population in the subcompartments. + Eigen::VectorXd init(infection_state.get_count()); + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::Susceptible)] = 750; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::Exposed)] = 30; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::Exposed) + 1] = 20; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms)] = 20; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms) + 1] = 10; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms) + 2] = 10; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSymptoms)] = 50; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSevere)] = 50; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical)] = 10; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 1] = 10; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 2] = 5; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 3] = 3; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 4] = 2; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::Recovered)] = 20; + init[infection_state.get_firstindex(mio::lsecir::InfectionStateBase::Dead)] = 10; + + // Initialize model. + mio::lsecir::Model model(std::move(init), infection_state); + + // Set Parameters. + model.parameters.get() = 3.2; + model.parameters.get() = 2; + model.parameters.get() = 5.8; + model.parameters.get() = 9.5; + // It is also possible to change values with the set function. + 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)); + // From SimulationTime 5, the contact pattern is reduced to 30% of the initial value. + contact_matrix[0].add_damping(0.7, mio::SimulationTime(5.)); + + model.parameters.get() = 0.7; + model.parameters.get() = 0.25; + 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 the InfectionState%s without subcompartments of the result and print it. + mio::TimeSeries population_no_subcompartments = model.calculate_populations(result); + population_no_subcompartments.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..b0ffc5a6f7 --- /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 only) +- P. J. Hurtado und A. S. Kirosingh, "Generalizations of the ‘Linear Chain Trick’: incorporating more flexible dwell time distributions into mean field ODE models“, 2019. (https://doi.org/10.1007/s00285-019-01412-w) + +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..02cdc8d6ef --- /dev/null +++ b/cpp/models/lct_secir/infection_state.h @@ -0,0 +1,232 @@ +/* +* 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 +}; + +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_subcompartment_numbers((int)InfectionStateBase::Count, 1) + , m_subcompartment_indexfirst((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] subcompartment_numbers Vector which defines the number of subcompartments for each infection state of InfectionStateBase. + */ + InfectionState(std::vector subcompartment_numbers) + : m_subcompartment_numbers(std::move(subcompartment_numbers)) + , m_subcompartment_indexfirst((int)InfectionStateBase::Count, 1) + { + check_constraints(); + set_compartment_index(); + } + + /** + * @brief Setter for the number of subcompartments. + * + * The number of subcompartments is only updated if the vector is valid. + * @param[in] subcompartment_numbers Vector which defines the number of subcompartments for each infection state of InfectionStateBase. + * @return Returns true if the function works as intended and false if the vector is not valid. + */ + bool set_subcompartment_numbers(std::vector subcompartment_numbers) + { + std::vector copy_m_subcompartment_numbers(m_subcompartment_numbers); + m_subcompartment_numbers = std::move(subcompartment_numbers); + if (check_constraints()) { + // Case where the vector is not valid. + log_warning("The vector you tried to set as the number of subcompartments is invalid. The previous vector " + "is kept."); + m_subcompartment_numbers = copy_m_subcompartment_numbers; + return false; + } + else { + set_compartment_index(); + return true; + } + } + + /** + * @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_subcompartment_numbers[(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. + * If the index does not match an infectionstate, the return value will be -1. + * @return Number of subcompartments for infectionstatebase or -1. + */ + int get_number(int infectionstatebaseindex) const + { + if ((0 <= infectionstatebaseindex) && (infectionstatebaseindex < (int)InfectionStateBase::Count)) { + return m_subcompartment_numbers[infectionstatebaseindex]; + } + else { + // Invalid index. + log_warning("The index you tried to get the number of subcompartments for was not valid."); + return -1; + } + } + + /** + * @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_subcompartment_indexfirst[(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. + * If the index does not match an infectionstate, the return value will be -1. + * @return Index of the first subcompartment for a vector with one entry per subcompartment or -1. + */ + int get_firstindex(int infectionstatebaseindex) const + { + if ((0 <= infectionstatebaseindex) && (infectionstatebaseindex < (int)InfectionStateBase::Count)) { + return m_subcompartment_indexfirst[infectionstatebaseindex]; + } + else { + return -1; + } + } + + /** + * @brief Gets the total number of subcompartments of all 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_subcompartment_numbers.size() == (int)InfectionStateBase::Count)) { + log_error("Vector for number of subcompartments has the wrong size."); + return true; + } + if (!(m_subcompartment_numbers[(int)InfectionStateBase::Susceptible] == 1)) { + log_error("Susceptible compartment can not have subcompartments."); + return true; + } + if (!(m_subcompartment_numbers[(int)InfectionStateBase::Recovered] == 1)) { + log_error("Recovered compartment can not have subcompartments."); + return true; + } + if (!(m_subcompartment_numbers[(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_subcompartment_numbers[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_subcompartment_indexfirst[i] = index; + index = index + m_subcompartment_numbers[i]; + } + m_count = index; + } + + std::vector + m_subcompartment_numbers; ///< Vector which defines the number of subcompartments for each infection state of InfectionStateBase. + std::vector + m_subcompartment_indexfirst; ///< 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 subcompartments of all 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..13c606a24d --- /dev/null +++ b/cpp/models/lct_secir/model.cpp @@ -0,0 +1,194 @@ +/* +* 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" + +namespace mio +{ +namespace lsecir +{ + +Model::Model(Eigen::VectorXd init, 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 result 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..6f0f4e1023 --- /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 differential equations of the form dydt = f(y, t). + * y is a vector containing number of individuals for each (sub-) compartment. + * This function evaluates the right-hand-side f of the ODE and can be used in an ODE solver. + * @param[in] 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 Cumulates a simulation result with subcompartments to produce a result that divides the population only into the infection states defined in InfectionStateBase. + * + * 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 the model. + 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..a2651f1480 --- /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, const Model& 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..d38e7ec3bf --- /dev/null +++ b/cpp/tests/test_lct_secir.cpp @@ -0,0 +1,567 @@ +/* +* 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. + mio::lsecir::InfectionState InfState1(std::vector((int)mio::lsecir::InfectionStateBase::Count + 1, 1)); + bool constraint_check = InfState1.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; + + std::vector SubcompartmentNumbers_copy1(SubcompartmentNumbers); + mio::lsecir::InfectionState InfState2(SubcompartmentNumbers_copy1); + constraint_check = InfState2.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; + std::vector SubcompartmentNumbers_copy2(SubcompartmentNumbers); + mio::lsecir::InfectionState InfState3(SubcompartmentNumbers_copy2); + constraint_check = InfState3.check_constraints(); + EXPECT_TRUE(constraint_check); + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Recovered] = 1; + + // For Dead. + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Dead] = 3; + std::vector SubcompartmentNumbers_copy3(SubcompartmentNumbers); + mio::lsecir::InfectionState InfState4(SubcompartmentNumbers_copy3); + constraint_check = InfState4.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; + std::vector SubcompartmentNumbers_copy4(SubcompartmentNumbers); + mio::lsecir::InfectionState InfState5(SubcompartmentNumbers_copy4); + constraint_check = InfState5.check_constraints(); + EXPECT_TRUE(constraint_check); + SubcompartmentNumbers[(int)mio::lsecir::InfectionStateBase::Exposed] = 2; + + // Check with correct parameters. + mio::lsecir::InfectionState InfState(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); +} + +// Test some setter and getter of the class InfectionState. +TEST(TestLCTSecir, testInfectionState) +{ + // Deactivate temporarily log output for next tests. + mio::set_log_level(mio::LogLevel::off); + + // Initialize InfectionState. + mio::lsecir::InfectionState InfState; + + // Try to set new number of subcompartments with an invalid vector. + EXPECT_FALSE( + InfState.set_subcompartment_numbers(std::vector((int)mio::lsecir::InfectionStateBase::Count - 1, 1))); + // Try to set new number of subcompartments with valid vector. + EXPECT_TRUE(InfState.set_subcompartment_numbers(std::vector((int)mio::lsecir::InfectionStateBase::Count, 1))); + + // Try to get the number of subcompartments for an invalid index. + EXPECT_EQ(-1, InfState.get_number(-1)); + // Valid index. + EXPECT_NE(-1, InfState.get_number(0)); + + // Try to get the index of the first subcompartment in a vector for an invalid index. + EXPECT_EQ(-1, InfState.get_firstindex(-1)); + // Valid index. + EXPECT_NE(-1, InfState.get_firstindex(0)); + + // Reactive log output. + mio::set_log_level(mio::LogLevel::warn); +} \ No newline at end of file