892 implement simple lct model

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)

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})

cpp/examples/lct_secir.cpp

@@ -0,0 +1,92 @@
+/* 
+* Copyright (C) 2020-2024 MEmilio
+*
+* Authors: Lena Ploetzke
+*
+* Contact: Martin J. Kuehn <[email protected]>
+*
+* 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 <vector>
+
+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<int> 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<mio::lsecir::TimeExposed>()            = 3.2;
+    model.parameters.get<mio::lsecir::TimeInfectedNoSymptoms>() = 2;
+    model.parameters.get<mio::lsecir::TimeInfectedSymptoms>()   = 5.8;
+    model.parameters.get<mio::lsecir::TimeInfectedSevere>()     = 9.5;
+    // It is also possible to change values with the set function.
+    model.parameters.set<mio::lsecir::TimeInfectedCritical>(7.1);
+
+    model.parameters.get<mio::lsecir::TransmissionProbabilityOnContact>() = 0.05;
+
+    mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::lsecir::ContactPatterns>();
+    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<mio::lsecir::RelativeTransmissionNoSymptoms>() = 0.7;
+    model.parameters.get<mio::lsecir::RiskOfInfectionFromSymptomatic>() = 0.25;
+    model.parameters.get<mio::lsecir::RecoveredPerInfectedNoSymptoms>() = 0.09;
+    model.parameters.get<mio::lsecir::SeverePerInfectedSymptoms>()      = 0.2;
+    model.parameters.get<mio::lsecir::CriticalPerSevere>()              = 0.25;
+    model.parameters.set<mio::lsecir::DeathsPerCritical>(0.3);
+
+    // Perform a simulation.
+    mio::TimeSeries<ScalarType> 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<ScalarType> population_no_subcompartments = model.calculate_populations(result);
+    population_no_subcompartments.print_table({"S", "E", "C", "I", "H", "U", "R", "D "}, 16, 8);
+}

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
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+target_compile_options(lct_secir PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

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. 
+
+

cpp/models/lct_secir/infection_state.h

@@ -0,0 +1,232 @@
+/* 
+* Copyright (C) 2020-2024 MEmilio
+*
+* Authors: Lena Ploetzke
+*
+* Contact: Martin J. Kuehn <[email protected]>
+*
+* 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 <vector>
+
+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<int> 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<int> subcompartment_numbers)
+    {
+        std::vector<int> 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<int>
+        m_subcompartment_numbers; ///< Vector which defines the number of subcompartments for each infection state of InfectionStateBase.
+    std::vector<int>
+        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