From afd49822b2bcba997b32d770ca8eead2be0faaa3 Mon Sep 17 00:00:00 2001 From: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> Date: Fri, 1 Mar 2024 15:22:02 +0100 Subject: [PATCH] 936 implement the class infectionstate of the lct model more efficiently (#941) - LctInfectionState is now a class template and has only static members (variables and functions). More computation is now done at compile time. - LctInfectionState can be used by various LCT models (not just SECIR models). Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com> --- cpp/examples/lct_secir.cpp | 95 +++++-- cpp/memilio/CMakeLists.txt | 1 + .../epidemiology/lct_infection_state.h | 86 ++++++ cpp/models/lct_secir/CMakeLists.txt | 3 +- cpp/models/lct_secir/README.md | 10 +- cpp/models/lct_secir/infection_state.h | 189 +------------ cpp/models/lct_secir/model.cpp | 165 ------------ cpp/models/lct_secir/model.h | 217 ++++++++++++++- cpp/models/lct_secir/simulation.cpp | 58 ---- cpp/models/lct_secir/simulation.h | 25 +- cpp/tests/test_lct_secir.cpp | 252 ++++++------------ 11 files changed, 476 insertions(+), 625 deletions(-) create mode 100644 cpp/memilio/epidemiology/lct_infection_state.h delete mode 100644 cpp/models/lct_secir/simulation.cpp diff --git a/cpp/examples/lct_secir.cpp b/cpp/examples/lct_secir.cpp index 7231356c33..3e24719850 100644 --- a/cpp/examples/lct_secir.cpp +++ b/cpp/examples/lct_secir.cpp @@ -25,42 +25,85 @@ #include "memilio/utils/time_series.h" #include "memilio/epidemiology/uncertain_matrix.h" #include "memilio/math/eigen.h" +#include "memilio/utils/logging.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. */ + // Simple example to demonstrate how to run a simulation using an LCT SECIR model. + // Parameters, initial values and the number of 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); + using Model = mio::lsecir::Model<2, 3, 1, 1, 5>; + using LctState = Model::LctState; 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; + // Define the initial value vector init with the distribution of the population into subcompartments. + // This method of defining the vector using a vector of vectors is a bit of overhead, but should remind you how + // the entries of the initial value vector relate to the defined template parameters of the model or the number of subcompartments. + // It is also possible to define the initial value vector directly. + std::vector> initial_populations = {{750}, {30, 20}, {20, 10, 10}, {50}, + {50}, {10, 10, 5, 3, 2}, {20}, {10}}; + + // Assert that initial_populations has the right shape. + if (initial_populations.size() != (int)LctState::InfectionState::Count) { + mio::log_error("The number of vectors in initial_populations does not match the number of InfectionStates."); + return 1; + } + if ((initial_populations[(int)LctState::InfectionState::Susceptible].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(int)LctState::InfectionState::Exposed].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(int)LctState::InfectionState::InfectedNoSymptoms].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(int)LctState::InfectionState::InfectedSymptoms].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(int)LctState::InfectionState::InfectedSevere].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(int)LctState::InfectionState::InfectedCritical].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(int)LctState::InfectionState::Recovered].size() != + LctState::get_num_subcompartments()) || + (initial_populations[(int)LctState::InfectionState::Dead].size() != + LctState::get_num_subcompartments())) { + mio::log_error("The length of at least one vector in initial_populations does not match the related number of " + "subcompartments."); + return 1; + } + + // Transfer the initial values in initial_populations to the vector init. + Eigen::VectorXd init = Eigen::VectorXd::Zero(LctState::Count); + init[(int)LctState::get_first_index()] = + initial_populations[(int)LctState::InfectionState::Susceptible][0]; + for (unsigned int i = 0; i < LctState::get_num_subcompartments(); i++) { + init[(int)LctState::get_first_index() + i] = + initial_populations[(int)LctState::InfectionState::Exposed][i]; + } + for (unsigned int i = 0; i < LctState::get_num_subcompartments(); + i++) { + init[(int)LctState::get_first_index() + i] = + initial_populations[(int)LctState::InfectionState::InfectedNoSymptoms][i]; + } + for (unsigned int i = 0; i < LctState::get_num_subcompartments(); i++) { + init[(int)LctState::get_first_index() + i] = + initial_populations[(int)LctState::InfectionState::InfectedSymptoms][i]; + } + for (unsigned int i = 0; i < LctState::get_num_subcompartments(); i++) { + init[(int)LctState::get_first_index() + i] = + initial_populations[(int)LctState::InfectionState::InfectedSevere][i]; + } + for (unsigned int i = 0; i < LctState::get_num_subcompartments(); i++) { + init[(int)LctState::get_first_index() + i] = + initial_populations[(int)LctState::InfectionState::InfectedCritical][i]; + } + init[(int)LctState::get_first_index()] = + initial_populations[(int)LctState::InfectionState::Recovered][0]; + init[(int)LctState::get_first_index()] = + initial_populations[(int)LctState::InfectionState::Dead][0]; // Initialize model. - mio::lsecir::Model model(std::move(init), infection_state); + Model model(std::move(init)); // Set Parameters. model.parameters.get() = 3.2; diff --git a/cpp/memilio/CMakeLists.txt b/cpp/memilio/CMakeLists.txt index be65a209c4..0bb0f81a88 100644 --- a/cpp/memilio/CMakeLists.txt +++ b/cpp/memilio/CMakeLists.txt @@ -16,6 +16,7 @@ add_library(memilio epidemiology/damping_sampling.cpp epidemiology/dynamic_npis.h epidemiology/dynamic_npis.cpp + epidemiology/lct_infection_state.h geography/regions.h geography/regions.cpp epidemiology/simulation_day.h diff --git a/cpp/memilio/epidemiology/lct_infection_state.h b/cpp/memilio/epidemiology/lct_infection_state.h new file mode 100644 index 0000000000..709e470465 --- /dev/null +++ b/cpp/memilio/epidemiology/lct_infection_state.h @@ -0,0 +1,86 @@ +/* +* Copyright (C) 2020-2024 MEmilio +* +* Authors: Lena Ploetzke +* +* Contact: Martin J. Kuehn +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#ifndef MIO_EPI_LCT_INFECTION_STATE_H +#define MIO_EPI_LCT_INFECTION_STATE_H + +#include + +namespace mio +{ +/** + * @brief Provides the functionality to be able to work with subcompartments in an LCT model. + * + * @tparam InfectionStates An enum class that defines the basic infection states. + * @tparam Ns Number of subcompartments for each infection state defined in InfectionState. + * The number of given template arguments must be equal to the entry Count from InfectionState. + */ +template +class LctInfectionState +{ +public: + using InfectionState = InfectionStates; + static_assert((unsigned int)InfectionState::Count == sizeof...(Ns), + "The number of integers provided as template parameters must be " + "the same as the entry Count of InfectionState."); + + static_assert(((Ns > 0) && ...), "The number of subcompartments must be at least 1."); + + /** + * @brief Gets the number of subcompartments in an infection state. + * + * @tparam State: Infection state for which the number of subcompartments should be returned. + * @return Number of subcompartments for State. + */ + template + static constexpr unsigned int get_num_subcompartments() + { + static_assert(State < InfectionState::Count, "State must be a a valid InfectionState."); + return m_subcompartment_numbers[(int)State]; + } + + /** + * @brief Gets the index of the first subcompartment of an infection state. + * + * In a simulation, the number of individuals in the subcompartments are stored in vectors. + * Accordingly, the index of the first subcompartment of State in such a vector is returned. + * @tparam State: Infection state for which the index should be returned. + * @return Index of the first subcompartment for a vector with one entry per subcompartment. + */ + template + static constexpr unsigned int get_first_index() + { + static_assert(State < InfectionState::Count, "State must be a a valid InfectionState."); + unsigned int index = 0; + for (int i = 0; i < (int)(State); i++) { + index = index + m_subcompartment_numbers[i]; + } + return index; + } + + static constexpr unsigned int Count{(... + Ns)}; + +private: + static constexpr const std::array m_subcompartment_numbers{ + Ns...}; ///< Vector which defines the number of subcompartments for each infection state of InfectionState. +}; + +} // namespace mio + +#endif \ No newline at end of file diff --git a/cpp/models/lct_secir/CMakeLists.txt b/cpp/models/lct_secir/CMakeLists.txt index e93fbe8bed..7c0aac951b 100644 --- a/cpp/models/lct_secir/CMakeLists.txt +++ b/cpp/models/lct_secir/CMakeLists.txt @@ -2,8 +2,7 @@ add_library(lct_secir infection_state.h model.h model.cpp - simulation.h - simulation.cpp + simulation.h parameters.h ) target_link_libraries(lct_secir PUBLIC memilio) diff --git a/cpp/models/lct_secir/README.md b/cpp/models/lct_secir/README.md index f8780f8ba4..25478c7b29 100644 --- a/cpp/models/lct_secir/README.md +++ b/cpp/models/lct_secir/README.md @@ -35,11 +35,11 @@ Below is an overview of the model architecture and its compartments. | $\xi_{I_{Sy}}$ | `RiskOfInfectionFromSymptomatic` | Proportion of infected people with symptomps who are not isolated. | | $N$ | `m_N0` | Total population. | | $D$ | `D` | Number of death people. | -| $n_E$ | Defined via `InfectionState` | Number of subcompartments of the Exposed compartment. | -| $n_{NS}$ | Defined via `InfectionState` | Number of subcompartments of the InfectedNoSymptoms compartment. | -| $n_{Sy}$ | Defined via `InfectionState` | Number of subcompartments of the InfectedSymptoms compartment. | -| $n_{Sev}$ | Defined via `InfectionState` | Number of subcompartments of the InfectedSevere compartment.| -| $n_{Cr}$ | Defined via `InfectionState` | Number of subcompartments of the InfectedCritical compartment. | +| $n_E$ | `NumExposed` | Number of subcompartments of the Exposed compartment. | +| $n_{NS}$ | `NumInfectedNoSymptoms` | Number of subcompartments of the InfectedNoSymptoms compartment. | +| $n_{Sy}$ | `NumInfectedSymptoms` | Number of subcompartments of the InfectedSymptoms compartment. | +| $n_{Sev}$ |`NumInfectedSevere` | Number of subcompartments of the InfectedSevere compartment.| +| $n_{Cr}$ | `NumInfectedCritical` | Number of subcompartments of the InfectedCritical compartment. | | $T_E$ | `TimeExposed` | Average time in days an individual stays in the Exposed compartment. | | $T_{I_{NS}}$ | `TimeInfectedNoSymptoms` | Average time in days an individual stays in the InfectedNoSymptoms compartment. | | $T_{I_{Sy}}$ | `TimeInfectedSymptoms` | Average time in days an individual stays in the InfectedSymptoms compartment. | diff --git a/cpp/models/lct_secir/infection_state.h b/cpp/models/lct_secir/infection_state.h index 02cdc8d6ef..383087e634 100644 --- a/cpp/models/lct_secir/infection_state.h +++ b/cpp/models/lct_secir/infection_state.h @@ -21,19 +21,16 @@ #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 + * @brief The InfectionState enum describes the basic + * categories for the infection state of persons. */ -enum class InfectionStateBase +enum class InfectionState { Susceptible = 0, Exposed = 1, @@ -46,186 +43,6 @@ enum class InfectionStateBase 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 diff --git a/cpp/models/lct_secir/model.cpp b/cpp/models/lct_secir/model.cpp index 13c606a24d..e804ebc447 100644 --- a/cpp/models/lct_secir/model.cpp +++ b/cpp/models/lct_secir/model.cpp @@ -19,176 +19,11 @@ */ #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 index 6f0f4e1023..0c2f67a769 100644 --- a/cpp/models/lct_secir/model.h +++ b/cpp/models/lct_secir/model.h @@ -25,31 +25,66 @@ #include "lct_secir/infection_state.h" #include "memilio/config.h" #include "memilio/utils/time_series.h" +#include "memilio/utils/logging.h" #include "memilio/math/eigen.h" -#include +#include "memilio/epidemiology/lct_infection_state.h" namespace mio { namespace lsecir { + +/** + * @brief Class that defines an LCT-SECIR model. + * + * @tparam NumExposed The number of subcompartents used for the Exposed compartment. + * @tparam NumInfectedNoSymptoms The number of subcompartents used for the InfectedNoSymptoms compartment. + * @tparam NumInfectedSymptoms The number of subcompartents used for the InfectedSymptoms compartment. + * @tparam NumInfectedSevere The number of subcompartents used for the InfectedSevere compartment. + * @tparam NumInfectedCritical The number of subcompartents used for the InfectedCritical compartment. + */ +template class Model { public: + using LctState = + LctInfectionState; ///< This class specifies the number of subcompartments. + /** * @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()); + Model(Eigen::VectorXd init, Parameters&& parameters_init = Parameters()) + : parameters{parameters_init} + , m_initial_values{std::move(init)} + { + m_N0 = m_initial_values.sum(); + } /** * @brief Checks constraints of the model inclusive check for model parameters. */ - bool check_constraints() const; + bool check_constraints() const + { + if (!(LctState::Count == m_initial_values.size())) { + log_error("Size of the initial values does not match subcompartments."); + return true; + } + for (unsigned int i = 0; i < LctState::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(); + } /** * @brief Evaulates the right-hand-side f of the LCT dydt = f(y, t). @@ -61,20 +96,181 @@ class 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; + void 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(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) + .sum(); + // Calculate sum of all subcompartments for InfectedSymptoms. + I = y.segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) + .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[LctState::template get_first_index()]) * 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 (Eigen::Index i = 0; i < LctState::template get_num_subcompartments(); 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 = LctState::template get_num_subcompartments() * + (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 (Eigen::Index i = 0; i < LctState::template get_num_subcompartments(); + i++) { + dummy = LctState::template get_num_subcompartments() * + (1 / parameters.get()) * + y[LctState::template get_first_index() + i]; + dydt[LctState::template get_first_index() + i] = + dydt[LctState::template get_first_index() + i] - dummy; + dydt[LctState::template get_first_index() + i + 1] = dummy; + } + + // I' + // Flow from last (sub-) compartment of C must be split between I_1 and R. + dydt[LctState::template get_first_index()] = + dydt[LctState::template get_first_index()] * + parameters.get(); + dydt[LctState::template get_first_index()] = + dydt[LctState::template get_first_index()] * + (1 - parameters.get()); + + for (Eigen::Index i = 0; i < LctState::template get_num_subcompartments(); + i++) { + dummy = LctState::template get_num_subcompartments() * + (1 / parameters.get()) * + y[LctState::template get_first_index() + i]; + dydt[LctState::template get_first_index() + i] = + dydt[LctState::template get_first_index() + i] - dummy; + dydt[LctState::template get_first_index() + i + 1] = dummy; + } + + // H' + dydt[LctState::template get_first_index()] = + dydt[LctState::template get_first_index()] + + dydt[LctState::template get_first_index()] * + (1 - parameters.get()); + dydt[LctState::template get_first_index()] = + dydt[LctState::template get_first_index()] * + parameters.get(); + for (Eigen::Index i = 0; i < LctState::template get_num_subcompartments(); + i++) { + dummy = LctState::template get_num_subcompartments() * + (1 / parameters.get()) * + y[LctState::template get_first_index() + i]; + dydt[LctState::template get_first_index() + i] = + dydt[LctState::template get_first_index() + i] - dummy; + dydt[LctState::template get_first_index() + i + 1] = dummy; + } + + // U' + dydt[LctState::template get_first_index()] = + dydt[LctState::template get_first_index()] + + dydt[LctState::template get_first_index()] * + (1 - parameters.get()); + dydt[LctState::template get_first_index()] = + dydt[LctState::template get_first_index()] * + parameters.get(); + for (Eigen::Index i = 0; i < LctState::template get_num_subcompartments() - 1; + i++) { + dummy = LctState::template get_num_subcompartments() * + (1 / parameters.get()) * + y[LctState::template get_first_index() + i]; + dydt[LctState::template get_first_index() + i] = + dydt[LctState::template get_first_index() + i] - dummy; + dydt[LctState::template get_first_index() + 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 = LctState::template get_num_subcompartments() * + (1 / parameters.get()) * + y[LctState::template get_first_index() - 1]; + dydt[LctState::template get_first_index() - 1] = + dydt[LctState::template get_first_index() - 1] - dummy; + dydt[LctState::template get_first_index()] = + dydt[LctState::template get_first_index()] + + (1 - parameters.get()) * dummy; + dydt[LctState::template get_first_index()] = parameters.get() * dummy; + } /** - * @brief Cumulates a simulation result with subcompartments to produce a result that divides the population only into the infection states defined in InfectionStateBase. + * @brief Cumulates a simulation result with subcompartments to produce a result that divides the population only into the infection states defined in InfectionState. * * If the model is used for simulation, we will get a result in form of a TimeSeries with infection states divided in subcompartments. - * Function transforms this TimeSeries in another TimeSeries with just the Basic infectionState. + * The function calculates a TimeSeries without subcompartmens from another TimeSeries with subcompartments. * 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; + TimeSeries calculate_populations(const TimeSeries& result) const + { + if (!(LctState::Count == result.get_num_elements())) { + log_error("Result does not match infectionState of the Model."); + TimeSeries populations((int)InfectionState::Count); + Eigen::VectorXd wrong_size = Eigen::VectorXd::Constant((int)InfectionState::Count, -1); + populations.add_time_point(-1, wrong_size); + return populations; + } + TimeSeries populations((int)InfectionState::Count); + Eigen::VectorXd dummy((int)InfectionState::Count); + for (Eigen::Index i = 0; i < result.get_num_time_points(); ++i) { + // Use segment of vector of the result with subcompartments of InfectionState with index j and sum up values of subcompartments. + dummy[(int)InfectionState::Susceptible] = result[i][0]; + dummy[(int)InfectionState::Exposed] = + result[i] + .segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) + .sum(); + dummy[(int)InfectionState::InfectedNoSymptoms] = + result[i] + .segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) + .sum(); + dummy[(int)InfectionState::InfectedSymptoms] = + result[i] + .segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) + .sum(); + dummy[(int)InfectionState::InfectedSevere] = + result[i] + .segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) + .sum(); + dummy[(int)InfectionState::InfectedCritical] = + result[i] + .segment(LctState::template get_first_index(), + LctState::template get_num_subcompartments()) + .sum(); + dummy[(int)InfectionState::Recovered] = + result[i][LctState::template get_first_index()]; + dummy[(int)InfectionState::Dead] = result[i][LctState::template get_first_index()]; + + populations.add_time_point(result.get_time(i), dummy); + } + + return populations; + } /** * @brief Returns the initial values for the model. @@ -88,7 +284,6 @@ class Model } Parameters parameters{}; ///< Parameters of the model. - InfectionState infectionState; ///< InfectionState specifies number of subcompartments. private: Eigen::VectorXd m_initial_values; ///< Initial values of the model. diff --git a/cpp/models/lct_secir/simulation.cpp b/cpp/models/lct_secir/simulation.cpp deleted file mode 100644 index a2651f1480..0000000000 --- a/cpp/models/lct_secir/simulation.cpp +++ /dev/null @@ -1,58 +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/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 index b66b1e4cb7..2e3369aeb6 100644 --- a/cpp/models/lct_secir/simulation.h +++ b/cpp/models/lct_secir/simulation.h @@ -27,6 +27,7 @@ #include "memilio/utils/metaprogramming.h" #include "memilio/math/stepper_wrapper.h" #include "memilio/math/eigen.h" +#include "boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp" namespace mio { @@ -35,6 +36,7 @@ namespace lsecir /** * @brief A class for the simulation of a LCT model. */ +template class Simulation { public: @@ -47,8 +49,15 @@ class Simulation * @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); - + Simulation(Model const& model, ScalarType t0 = 0., ScalarType dt = 0.1) + : 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) + { + } /** * @brief Set the core integrator used in the simulation. * @@ -168,8 +177,18 @@ class Simulation * 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. */ +template > TimeSeries simulate(ScalarType t0, ScalarType tmax, ScalarType dt, Model const& model, - std::shared_ptr integrator = nullptr); + std::shared_ptr integrator = nullptr) +{ + model.check_constraints(); + Sim sim(model, t0, dt); + if (integrator) { + sim.set_integrator(integrator); + } + sim.advance(tmax); + return sim.get_result(); +} } // namespace lsecir } // namespace mio diff --git a/cpp/tests/test_lct_secir.cpp b/cpp/tests/test_lct_secir.cpp index d38e7ec3bf..32263cb8f8 100644 --- a/cpp/tests/test_lct_secir.cpp +++ b/cpp/tests/test_lct_secir.cpp @@ -30,22 +30,22 @@ #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) { + using Model = mio::lsecir::Model<1, 1, 1, 1, 1>; ScalarType t0 = 0; ScalarType tmax = 1; ScalarType dt = 0.1; - Eigen::VectorXd init = Eigen::VectorXd::Constant((int)mio::lsecir::InfectionStateBase::Count, 15); + Eigen::VectorXd init = Eigen::VectorXd::Constant((int)Model::LctState::InfectionState::Count, 15); init[0] = 200; init[3] = 50; init[5] = 30; - mio::lsecir::Model model(init); + Model model(init); mio::TimeSeries result = mio::lsecir::simulate(t0, tmax, dt, model); EXPECT_NEAR(result.get_last_time(), tmax, 1e-10); @@ -59,21 +59,22 @@ TEST(TestLCTSecir, simulateDefault) with the result of the equivalent ODE SECIR model. */ TEST(TestLCTSecir, compareWithOdeSecir) { + using Model = mio::lsecir::Model<1, 1, 1, 1, 1>; 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); + Eigen::VectorXd init = Eigen::VectorXd::Constant((int)Model::LctState::InfectionState::Count, 15); init[0] = 200; init[3] = 50; init[5] = 30; // Define LCT model. - mio::lsecir::Model model_lct(init); + 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() = 3.2; + model_lct.parameters.get() = 2; model_lct.parameters.get() = 5.8; model_lct.parameters.get() = 9.5; model_lct.parameters.get() = 7.1; @@ -102,21 +103,21 @@ TEST(TestLCTSecir, compareWithOdeSecir) 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)]; + init[Eigen::Index(Model::LctState::InfectionState::Exposed)]; model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] = - init[Eigen::Index(mio::lsecir::InfectionStateBase::InfectedNoSymptoms)]; + init[Eigen::Index(Model::LctState::InfectionState::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)]; + init[Eigen::Index(Model::LctState::InfectionState::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)]; + init[Eigen::Index(Model::LctState::InfectionState::InfectedSevere)]; model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}] = - init[Eigen::Index(mio::lsecir::InfectionStateBase::InfectedCritical)]; + init[Eigen::Index(Model::LctState::InfectionState::InfectedCritical)]; model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Recovered}] = - init[Eigen::Index(mio::lsecir::InfectionStateBase::Recovered)]; + init[Eigen::Index(Model::LctState::InfectionState::Recovered)]; model_ode.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Dead}] = - init[Eigen::Index(mio::lsecir::InfectionStateBase::Dead)]; + init[Eigen::Index(Model::LctState::InfectionState::Dead)]; model_ode.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}, init.sum()); @@ -157,23 +158,23 @@ TEST(TestLCTSecir, compareWithOdeSecir) 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], + ASSERT_NEAR(result_lct[i][(int)Model::LctState::InfectionState::Susceptible], result_ode[i][(int)mio::osecir::InfectionState::Susceptible], 1e-5); - ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::Exposed], + ASSERT_NEAR(result_lct[i][(int)Model::LctState::InfectionState::Exposed], result_ode[i][(int)mio::osecir::InfectionState::Exposed], 1e-5); - ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::InfectedNoSymptoms], + ASSERT_NEAR(result_lct[i][(int)Model::LctState::InfectionState::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], + ASSERT_NEAR(result_lct[i][(int)Model::LctState::InfectionState::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], + ASSERT_NEAR(result_lct[i][(int)Model::LctState::InfectionState::InfectedCritical], result_ode[i][(int)mio::osecir::InfectionState::InfectedCritical], 1e-5); - ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::InfectedSevere], + ASSERT_NEAR(result_lct[i][(int)Model::LctState::InfectionState::InfectedSevere], result_ode[i][(int)mio::osecir::InfectionState::InfectedSevere], 1e-5); - ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::Recovered], + ASSERT_NEAR(result_lct[i][(int)Model::LctState::InfectionState::Recovered], result_ode[i][(int)mio::osecir::InfectionState::Recovered], 1e-5); - ASSERT_NEAR(result_lct[i][(int)mio::lsecir::InfectionStateBase::Dead], + ASSERT_NEAR(result_lct[i][(int)Model::LctState::InfectionState::Dead], result_ode[i][(int)mio::osecir::InfectionState::Dead], 1e-5); } } @@ -182,38 +183,32 @@ TEST(TestLCTSecir, compareWithOdeSecir) 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); + // Chose more than one subcompartment for all compartments except S, R, D so that the function is correct for all selections. + using Model = mio::lsecir::Model<2, 3, 2, 2, 2>; + using LctState = Model::LctState; // 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); + Eigen::VectorXd init(LctState::Count); + init[LctState::get_first_index()] = 750; + init[LctState::get_first_index()] = 30; + init[LctState::get_first_index() + 1] = 20; + init[LctState::get_first_index()] = 20; + init[LctState::get_first_index() + 1] = 10; + init[LctState::get_first_index() + 2] = 10; + init[LctState::get_first_index()] = 30; + init[LctState::get_first_index() + 1] = 20; + init[LctState::get_first_index()] = 40; + init[LctState::get_first_index() + 1] = 10; + init[LctState::get_first_index()] = 10; + init[LctState::get_first_index() + 1] = 20; + init[LctState::get_first_index()] = 20; + init[LctState::get_first_index()] = 10; + + Model model(std::move(init)); // Set parameters. - model.parameters.set(2 * 4.2 - 5.2); - model.parameters.get() = 2 * (5.2 - 4.2); + model.parameters.set(3.2); + model.parameters.get() = 2; model.parameters.get() = 5.8; model.parameters.get() = 9.5; model.parameters.get() = 7.1; @@ -233,7 +228,7 @@ TEST(TestLCTSecir, testEvalRightHandSide) 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(); + unsigned int num_subcompartments = LctState::Count; Eigen::VectorXd dydt(num_subcompartments); model.eval_right_hand_side(model.get_initial_values(), 0, dydt); @@ -241,7 +236,7 @@ TEST(TestLCTSecir, testEvalRightHandSide) 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++) { + for (unsigned int i = 0; i < num_subcompartments; i++) { ASSERT_NEAR(compare[i], dydt[i], 1e-3); } } @@ -249,38 +244,35 @@ TEST(TestLCTSecir, testEvalRightHandSide) // Model setup to compare result with a previous output. class ModelTestLCTSecir : public testing::Test { +public: + using Model = mio::lsecir::Model<2, 3, 1, 1, 5>; + using LctState = Model::LctState; + 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; + // Define initial distribution of the population in the subcompartments. + Eigen::VectorXd init(LctState::Count); + init[LctState::get_first_index()] = 750; + init[LctState::get_first_index()] = 30; + init[LctState::get_first_index() + 1] = 20; + init[LctState::get_first_index()] = 20; + init[LctState::get_first_index() + 1] = 10; + init[LctState::get_first_index() + 2] = 10; + init[LctState::get_first_index()] = 50; + init[LctState::get_first_index()] = 50; + init[LctState::get_first_index()] = 10; + init[LctState::get_first_index() + 1] = 10; + init[LctState::get_first_index() + 2] = 5; + init[LctState::get_first_index() + 3] = 3; + init[LctState::get_first_index() + 4] = 2; + init[LctState::get_first_index()] = 20; + init[LctState::get_first_index()] = 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 = new Model(std::move(init)); + model->parameters.get() = 3.2; + model->parameters.get() = 2; model->parameters.get() = 5.8; model->parameters.get() = 9.5; model->parameters.get() = 7.1; @@ -304,7 +296,7 @@ class ModelTestLCTSecir : public testing::Test } public: - mio::lsecir::Model* model = nullptr; + Model* model = nullptr; }; // Test compares a simulation with the result of a previous run stored in a .csv file. @@ -327,7 +319,7 @@ TEST_F(ModelTestLCTSecir, compareWithPreviousRun) } } - // Compare base compartments. + // Compare InfectionState compartments. mio::TimeSeries population = model->calculate_populations(result); auto compare_population = load_test_data_csv("lct-secir-compartments-compare.csv"); @@ -347,8 +339,8 @@ 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); + mio::TimeSeries init_wrong_size((int)LctState::InfectionState::Count); + Eigen::VectorXd vec_wrong_size = Eigen::VectorXd::Ones((int)LctState::InfectionState::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); @@ -360,61 +352,12 @@ TEST_F(ModelTestLCTSecir, testCalculatePopWrongSize) mio::set_log_level(mio::LogLevel::warn); } -// Check constraints of InfectionState and Parameters. +// Check constraints of Parameters and Model. 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; @@ -434,7 +377,7 @@ TEST(TestLCTSecir, testConstraints) parameters_lct.get() = 0.1; // Check improper TimeExposed. - constraint_check = parameters_lct.check_constraints(); + bool constraint_check = parameters_lct.check_constraints(); EXPECT_TRUE(constraint_check); parameters_lct.get() = 3.1; @@ -515,53 +458,24 @@ TEST(TestLCTSecir, testConstraints) EXPECT_FALSE(constraint_check); // Check for model. + using Model = mio::lsecir::Model<1, 1, 1, 1, 1>; + using LctState = Model::LctState; + // 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)); + Model model1(std::move(Eigen::VectorXd::Ones((int)LctState::Count - 1)), 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)); + Model model2(std::move(Eigen::VectorXd::Constant((int)LctState::Count, -1)), 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)); + Model model3(std::move(Eigen::VectorXd::Constant((int)LctState::Count, 1))); 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