680 implement lct model

cpp/CMakeLists.txt

@@ -120,6 +120,7 @@ if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/ode_secir)
     add_subdirectory(models/ode_secirvvs)
     add_subdirectory(models/ide_secir)
+    add_subdirectory(models/lct_secir)
     add_subdirectory(models/ide_seir)
     add_subdirectory(models/ode_seir)
 endif()

cpp/examples/CMakeLists.txt

@@ -58,6 +58,13 @@ 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(compute_parameters compute_parameters.cpp)
+target_link_libraries(compute_parameters PRIVATE memilio)
+
+add_executable(lct_secir_example lct_secir.cpp)
+target_link_libraries(lct_secir_example PRIVATE memilio lct_secir)
+
+
 if(MEMILIO_HAS_JSONCPP)
   add_executable(ode_secir_read_graph_example ode_secir_read_graph.cpp)
   target_link_libraries(ode_secir_read_graph_example PRIVATE memilio ode_secir)
@@ -88,3 +95,18 @@ if(MEMILIO_HAS_HDF5)
   target_link_libraries(ode_secir_save_results_example PRIVATE memilio ode_secir)
   target_compile_options(ode_secir_save_results_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 endif()
+
+if(MEMILIO_HAS_HDF5)
+  add_executable(lct_secir_fictional_scenario lct_secir_fictional_scenario.cpp)
+  target_link_libraries(lct_secir_fictional_scenario PRIVATE memilio lct_secir ide_secir)
+endif()
+
+if(MEMILIO_HAS_HDF5)
+  add_executable(lct_secir_real_scenario lct_secir_real_scenario.cpp)
+  target_link_libraries(lct_secir_real_scenario PRIVATE memilio lct_secir)
+endif()
+
+if(MEMILIO_HAS_HDF5)
+  add_executable(lct_secir_initializations lct_secir_initializations.cpp)
+  target_link_libraries(lct_secir_initializations PRIVATE memilio lct_secir ide_secir)
+endif()

cpp/examples/compute_parameters.cpp

@@ -0,0 +1,209 @@
+/* 
+* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC)
+*
+* Authors: Lena Ploetzke
+*
+* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* 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 <iostream>
+
+int main()
+{
+    /** With this file parameters without age distribution can be calculated to match those specified in the first published paper
+     (https://doi.org/10.1016/j.mbs.2021.108648). 
+
+    For ODE or LCT models, the parameters from the paper are used first to derive just one transitiontime if 2 are specified, eg T_C=\mu_C^R*T_C^R+(1-\mu_C^R)*T_C^I.
+    For an IDE mdel, this step is omitted.
+    For each age group, average values are calculated if lower and upper bounds are given. 
+    For this we assume a uniform distribution so that T=(T_min+T_max)/2.
+    Finally we calculate a weighted average time across the age groups.
+    */
+    bool printResult = true;
+
+    // Age group sizes are calculated using table number 12411-04-02-4-B from www.regionalstatistik.de for the date 31.12.2020.
+    const double age_group_sizes[] = {3969138.0, 7508662, 18921292, 28666166, 18153339, 5936434};
+    const int total                = 83155031.0;
+    const int numagegroups         = 6;
+
+    // Transmission parameters.
+    const double transmissionProbabilityOnContactMin[] = {0.02, 0.05, 0.05, 0.05, 0.08, 0.15};
+    const double transmissionProbabilityOnContactMax[] = {0.04, 0.07, 0.07, 0.07, 0.10, 0.20};
+
+    double transmissionProbabilityOnContact = 0;
+    for (int i = 0; i < numagegroups; i++) {
+        transmissionProbabilityOnContact +=
+            (age_group_sizes[i] / total) * 0.5 *
+            (transmissionProbabilityOnContactMin[i] + transmissionProbabilityOnContactMax[i]);
+    }
+
+    if (printResult) {
+        std::cout << "transmissionProbabilityOnContact: " << transmissionProbabilityOnContact << std::endl;
+    }
+
+    const double relativeTransmissionNoSymptoms = 1;
+    const double riskOfInfectionFromSymptomatic = 0.3;
+
+    if (printResult) {
+        std::cout << "relativeTransmissionNoSymptoms: " << relativeTransmissionNoSymptoms << std::endl;
+        std::cout << "riskOfInfectionFromSymptomatic: " << riskOfInfectionFromSymptomatic << std::endl;
+    }
+
+    // E
+    const double timeExposedMin = 2.67;
+    const double timeExposedMax = 4.00;
+
+    double timeExposed = 0.5 * (timeExposedMin + timeExposedMax);
+
+    if (printResult) {
+        std::cout << "timeExposed: " << timeExposed << std::endl;
+    }
+
+    // Calculate parameters for I first because a value of I is needed for C.
+    // I
+    const double timeInfectedSymptomstoRecoveredMin        = 5.6;
+    const double timeInfectedSymptomstoRecoveredMax        = 8.4;
+    const double timeInfectedSymptomstoInfectedSevereMin[] = {9, 9, 9, 5, 5, 5};
+    const double timeInfectedSymptomstoInfectedSevereMax[] = {12, 12, 12, 7, 7, 7};
+    const double severePerInfectedSymptomsMin[]            = {0.006, 0.006, 0.015, 0.049, 0.15, 0.20};
+    const double severePerInfectedSymptomsMax[]            = {0.009, 0.009, 0.023, 0.074, 0.18, 0.25};
+
+    double timeInfectedSymptomsMindummy;
+    double timeInfectedSymptomsMaxdummy;
+    double severePerInfectedSymptomsdummy;
+    double timeInfectedSymptoms      = 0;
+    double severePerInfectedSymptoms = 0;
+
+    for (int i = 0; i < numagegroups; i++) {
+        severePerInfectedSymptomsdummy = 0.5 * (severePerInfectedSymptomsMin[i] + severePerInfectedSymptomsMax[i]);
+        timeInfectedSymptomsMindummy   = (1 - severePerInfectedSymptomsdummy) * timeInfectedSymptomstoRecoveredMin +
+                                       severePerInfectedSymptomsdummy * timeInfectedSymptomstoInfectedSevereMin[i];
+        timeInfectedSymptomsMaxdummy = (1 - severePerInfectedSymptomsdummy) * timeInfectedSymptomstoRecoveredMax +
+                                       severePerInfectedSymptomsdummy * timeInfectedSymptomstoInfectedSevereMax[i];
+
+        timeInfectedSymptoms +=
+            (age_group_sizes[i] / total) * 0.5 * (timeInfectedSymptomsMindummy + timeInfectedSymptomsMaxdummy);
+        severePerInfectedSymptoms += (age_group_sizes[i] / total) * severePerInfectedSymptomsdummy;
+    }
+
+    // Calculation for an IDE model.
+    const double timeInfectedSymptomstoRecovered =
+        0.5 * (timeInfectedSymptomstoRecoveredMin + timeInfectedSymptomstoRecoveredMax);
+    double timeInfectedSymptomstoInfectedSevere = 0;
+
+    for (int i = 0; i < numagegroups; i++) {
+        severePerInfectedSymptomsdummy = 0.5 * (severePerInfectedSymptomsMin[i] + severePerInfectedSymptomsMax[i]);
+        timeInfectedSymptomstoInfectedSevere +=
+            severePerInfectedSymptomsdummy * (age_group_sizes[i] / total) * 0.5 *
+            (timeInfectedSymptomstoInfectedSevereMin[i] + timeInfectedSymptomstoInfectedSevereMax[i]);
+    }
+    timeInfectedSymptomstoInfectedSevere = timeInfectedSymptomstoInfectedSevere / severePerInfectedSymptoms;
+    // C
+    const double timeInfectedNoSymptomstoInfectedSymptoms = 5.2 - timeExposed;
+    const double timeInfectedNoSymptomstoRecovered =
+        timeInfectedNoSymptomstoInfectedSymptoms +
+        0.5 * (timeInfectedSymptomstoRecoveredMin + timeInfectedSymptomstoRecoveredMax);
+    const double recoveredPerInfectedNoSymptomsMin[] = {0.2, 0.2, 0.15, 0.15, 0.15, 0.15};
+    const double recoveredPerInfectedNoSymptomsMax[] = {0.3, 0.3, 0.25, 0.25, 0.25, 0.25};
+
+    double timeInfectedNoSymptoms         = 0;
+    double recoveredPerInfectedNoSymptoms = 0;
+
+    for (int i = 0; i < numagegroups; i++) {
+        recoveredPerInfectedNoSymptoms += (age_group_sizes[i] / total) * 0.5 *
+                                          (recoveredPerInfectedNoSymptomsMin[i] + recoveredPerInfectedNoSymptomsMax[i]);
+    }
+
+    timeInfectedNoSymptoms = recoveredPerInfectedNoSymptoms * timeInfectedNoSymptomstoRecovered +
+                             (1 - recoveredPerInfectedNoSymptoms) * timeInfectedNoSymptomstoInfectedSymptoms;
+
+    if (printResult) {
+        std::cout << "timeInfectedNoSymptoms: " << timeInfectedNoSymptoms << std::endl;
+        std::cout << "recoveredPerInfectedNoSymptoms: " << recoveredPerInfectedNoSymptoms << std::endl;
+        std::cout << "timeInfectedSymptoms: " << timeInfectedSymptoms << std::endl;
+        std::cout << "severePerInfectedSymptoms: " << severePerInfectedSymptoms << std::endl;
+    }
+
+    // H
+    const double timeInfectedSeveretoRecoveredMin[]      = {4, 4, 5, 7, 9, 13};
+    const double timeInfectedSeveretoRecoveredMax[]      = {6, 6, 7, 9, 11, 17};
+    const double timeInfectedSeveretoInfectedCriticalMin = 3;
+    const double timeInfectedSeveretoInfectedCriticalMax = 7;
+    const double criticalPerSevereMin[]                  = {0.05, 0.05, 0.05, 0.10, 0.25, 0.35};
+    const double criticalPerSevereMax[]                  = {0.10, 0.10, 0.10, 0.20, 0.35, 0.45};
+
+    double timeInfectedSevereMindummy;
+    double timeInfectedSevereMaxdummy;
+    double criticalPerSeveredummy;
+    double timeInfectedSevere = 0;
+    double criticalPerSevere  = 0;
+
+    for (int i = 0; i < numagegroups; i++) {
+        criticalPerSeveredummy     = 0.5 * (criticalPerSevereMin[i] + criticalPerSevereMax[i]);
+        timeInfectedSevereMindummy = (1 - criticalPerSeveredummy) * timeInfectedSeveretoRecoveredMin[i] +
+                                     criticalPerSeveredummy * timeInfectedSeveretoInfectedCriticalMin;
+        timeInfectedSevereMaxdummy = (1 - criticalPerSeveredummy) * timeInfectedSeveretoRecoveredMax[i] +
+                                     criticalPerSeveredummy * timeInfectedSeveretoInfectedCriticalMax;
+
+        timeInfectedSevere +=
+            (age_group_sizes[i] / total) * 0.5 * (timeInfectedSevereMindummy + timeInfectedSevereMaxdummy);
+        criticalPerSevere += (age_group_sizes[i] / total) * criticalPerSeveredummy;
+    }
+
+    if (printResult) {
+        std::cout << "timeInfectedSevere: " << timeInfectedSevere << std::endl;
+        std::cout << "criticalPerSevere: " << criticalPerSevere << std::endl;
+    }
+
+    // U
+    const double timeInfectedCriticaltoRecoveredMin[] = {5, 5, 5, 14, 14, 10};
+    const double timeInfectedCriticaltoRecoveredMax[] = {9, 9, 9, 21, 21, 15};
+    const double timeInfectedCriticaltoDeadMin[]      = {4, 4, 4, 15, 15, 10};
+    const double timeInfectedCriticaltoDeadMax[]      = {8, 8, 8, 18, 18, 12};
+    const double deathsPerCriticalMin[]               = {0.00, 0.00, 0.10, 0.10, 0.30, 0.5};
+    const double deathsPerCriticalMax[]               = {0.10, 0.10, 0.18, 0.18, 0.50, 0.7};
+
+    double timeInfectedCriticalMindummy;
+    double timeInfectedCriticalMaxdummy;
+    double deathsPerCriticaldummy;
+    double timeInfectedCritical = 0;
+    double deathsPerCritical    = 0;
+
+    for (int i = 0; i < numagegroups; i++) {
+        deathsPerCriticaldummy       = 0.5 * (deathsPerCriticalMin[i] + deathsPerCriticalMax[i]);
+        timeInfectedCriticalMindummy = (1 - deathsPerCriticaldummy) * timeInfectedCriticaltoRecoveredMin[i] +
+                                       deathsPerCriticaldummy * timeInfectedCriticaltoDeadMin[i];
+        timeInfectedCriticalMaxdummy = (1 - deathsPerCriticaldummy) * timeInfectedCriticaltoRecoveredMax[i] +
+                                       deathsPerCriticaldummy * timeInfectedCriticaltoDeadMax[i];
+
+        timeInfectedCritical +=
+            (age_group_sizes[i] / total) * 0.5 * (timeInfectedCriticalMindummy + timeInfectedCriticalMaxdummy);
+        deathsPerCritical += (age_group_sizes[i] / total) * deathsPerCriticaldummy;
+    }
+
+    if (printResult) {
+        std::cout << "timeInfectedCritical: " << timeInfectedCritical << std::endl;
+        std::cout << "deathsPerCritical: " << deathsPerCritical << std::endl;
+    }
+
+    if (printResult) {
+        std::cout << "\nFor IDE model: " << std::endl;
+        std::cout << "timeInfectedNoSymptomstoInfectedSymptoms: " << timeInfectedNoSymptomstoInfectedSymptoms
+                  << std::endl;
+        std::cout << "timeInfectedNoSymptomstoRecovered: " << timeInfectedNoSymptomstoRecovered << std::endl;
+        std::cout << "timeInfectedSymptomstoInfectedSevere: " << timeInfectedSymptomstoInfectedSevere << std::endl;
+        std::cout << "timeInfectedSymptomstoRecovered: " << timeInfectedSymptomstoRecovered << std::endl;
+    }
+}

cpp/examples/ide_secir.cpp

@@ -32,10 +32,10 @@ int main()
 {
     using Vec = mio::TimeSeries<ScalarType>::Vector;
 
-    ScalarType tmax        = 10;
-    ScalarType N           = 10000;
-    ScalarType Dead_before = 12;
-    ScalarType dt          = 1;
+    ScalarType tmax   = 10;
+    ScalarType N      = 10000;
+    ScalarType deaths = 13.10462213;
+    ScalarType dt     = 1;
 
     int num_transitions = (int)mio::isecir::InfectionTransition::Count;
 
@@ -63,7 +63,7 @@ int main()
     }
 
     // Initialize model.
-    mio::isecir::Model model(std::move(init), N, Dead_before);
+    mio::isecir::Model model(std::move(init), N, deaths);
 
     // model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 1000;
     // model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered]   = 0;
@@ -100,4 +100,4 @@ int main()
     sim.print_transitions();
 
     sim.print_compartments();
-}
+}

cpp/examples/lct_secir.cpp

@@ -0,0 +1,91 @@
+/* 
+* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC)
+*
+* Authors: Lena Ploetzke
+*
+* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* 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 simulate using an LCT SECIR model. 
+    Parameters, initial values and subcompartments are not realistic. */
+
+    // 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 infectionState(num_subcompartments);
+
+    ScalarType tmax = 20;
+
+    // Define initial distribution of the population in the subcompartments.
+    Eigen::VectorXd init(infectionState.get_count());
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::Susceptible)]            = 750;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::Exposed)]                = 30;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::Exposed) + 1]            = 20;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms)]     = 20;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms) + 1] = 10;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedNoSymptoms) + 2] = 10;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSymptoms)]       = 50;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedSevere)]         = 50;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical)]       = 10;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 1]   = 10;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 2]   = 5;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 3]   = 3;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::InfectedCritical) + 4]   = 2;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::Recovered)]              = 20;
+    init[infectionState.get_firstindex(mio::lsecir::InfectionStateBase::Dead)]                   = 10;
+
+    // Initialize model.
+    mio::lsecir::Model model(std::move(init), infectionState);
+
+    // Set Parameters.
+    model.parameters.get<mio::lsecir::TimeExposed>()            = 2 * 4.2 - 5.2;
+    model.parameters.get<mio::lsecir::TimeInfectedNoSymptoms>() = 2 * (5.2 - 4.2);
+    model.parameters.get<mio::lsecir::TimeInfectedSymptoms>()   = 5.8;
+    model.parameters.get<mio::lsecir::TimeInfectedSevere>()     = 9.5;
+    // Also possible to change values with setter.
+    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));
+    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 infectionState without subcompartments of the result and print it.
+    mio::TimeSeries<ScalarType> populations = model.calculate_populations(result);
+    mio::lsecir::print_TimeSeries(populations, model.get_heading_CompartmentsBase());
+}