323 SECIRTS model with waning immunity

cpp/CMakeLists.txt

@@ -144,6 +144,7 @@ add_subdirectory(memilio)
 if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/abm)
     add_subdirectory(models/ode_secir)
+    add_subdirectory(models/ode_secirts)
     add_subdirectory(models/ode_secirvvs)
     add_subdirectory(models/lct_secir)
     add_subdirectory(models/ide_secir)

cpp/benchmarks/flow_simulation_ode_secirvvs.h

@@ -486,15 +486,17 @@ class Simulation : public Base
             double first_vacc;
             double full_vacc;
             if (t_idx == SimulationDay(0)) {
-                first_vacc = params.template get<osecirvvs::DailyFirstVaccination<ScalarType>>()[{(AgeGroup)i, t_idx}];
-                full_vacc  = params.template get<osecirvvs::DailyFullVaccination<ScalarType>>()[{(AgeGroup)i, t_idx}];
+                first_vacc =
+                    params.template get<osecirvvs::DailyPartialVaccinations<ScalarType>>()[{(AgeGroup)i, t_idx}];
+                full_vacc = params.template get<osecirvvs::DailyFullVaccinations<ScalarType>>()[{(AgeGroup)i, t_idx}];
             }
             else {
-                first_vacc = params.template get<osecirvvs::DailyFirstVaccination<ScalarType>>()[{(AgeGroup)i, t_idx}] -
-                             params.template get<osecirvvs::DailyFirstVaccination<ScalarType>>()[{
-                                 (AgeGroup)i, t_idx - SimulationDay(1)}];
-                full_vacc = params.template get<osecirvvs::DailyFullVaccination<ScalarType>>()[{(AgeGroup)i, t_idx}] -
-                            params.template get<osecirvvs::DailyFullVaccination<ScalarType>>()[{
+                first_vacc =
+                    params.template get<osecirvvs::DailyPartialVaccinations<ScalarType>>()[{(AgeGroup)i, t_idx}] -
+                    params.template get<osecirvvs::DailyPartialVaccinations<ScalarType>>()[{(AgeGroup)i,
+                                                                                            t_idx - SimulationDay(1)}];
+                full_vacc = params.template get<osecirvvs::DailyFullVaccinations<ScalarType>>()[{(AgeGroup)i, t_idx}] -
+                            params.template get<osecirvvs::DailyFullVaccinations<ScalarType>>()[{
                                 (AgeGroup)i, t_idx - SimulationDay(1)}];
             }
 
@@ -665,10 +667,11 @@ void setup_model(Model& model)
 
     model.parameters.template get<osecirvvs::ICUCapacity<ScalarType>>()          = 100;
     model.parameters.template get<osecirvvs::TestAndTraceCapacity<ScalarType>>() = 0.0143;
-    model.parameters.template get<osecirvvs::DailyFirstVaccination<ScalarType>>().resize(SimulationDay(size_t(1000)));
-    model.parameters.template get<osecirvvs::DailyFirstVaccination<ScalarType>>().array().setConstant(5);
-    model.parameters.template get<osecirvvs::DailyFullVaccination<ScalarType>>().resize(SimulationDay(size_t(1000)));
-    model.parameters.template get<osecirvvs::DailyFullVaccination<ScalarType>>().array().setConstant(3);
+    model.parameters.template get<osecirvvs::DailyPartialVaccinations<ScalarType>>().resize(
+        SimulationDay(size_t(1000)));
+    model.parameters.template get<osecirvvs::DailyPartialVaccinations<ScalarType>>().array().setConstant(5);
+    model.parameters.template get<osecirvvs::DailyFullVaccinations<ScalarType>>().resize(SimulationDay(size_t(1000)));
+    model.parameters.template get<osecirvvs::DailyFullVaccinations<ScalarType>>().array().setConstant(3);
 
     auto& contacts       = model.parameters.template get<osecirvvs::ContactPatterns<ScalarType>>();
     auto& contact_matrix = contacts.get_cont_freq_mat();

cpp/benchmarks/graph_simulation.cpp

@@ -69,10 +69,10 @@ mio::osecirvvs::Model<ScalarType> create_model(size_t num_agegroups, const Scala
     const size_t vacc_full                                                   = 5;
     model.parameters.get<mio::osecirvvs::ICUCapacity<ScalarType>>()          = 100;
     model.parameters.get<mio::osecirvvs::TestAndTraceCapacity<ScalarType>>() = 0.0143;
-    model.parameters.get<mio::osecirvvs::DailyFirstVaccination<ScalarType>>().resize(mio::SimulationDay(tmax));
-    model.parameters.get<mio::osecirvvs::DailyFirstVaccination<ScalarType>>().array().setConstant(vacc_first);
-    model.parameters.get<mio::osecirvvs::DailyFullVaccination<ScalarType>>().resize(mio::SimulationDay(tmax));
-    model.parameters.get<mio::osecirvvs::DailyFullVaccination<ScalarType>>().array().setConstant(vacc_full);
+    model.parameters.get<mio::osecirvvs::DailyPartialVaccinations<ScalarType>>().resize(mio::SimulationDay(tmax));
+    model.parameters.get<mio::osecirvvs::DailyPartialVaccinations<ScalarType>>().array().setConstant(vacc_first);
+    model.parameters.get<mio::osecirvvs::DailyFullVaccinations<ScalarType>>().resize(mio::SimulationDay(tmax));
+    model.parameters.get<mio::osecirvvs::DailyFullVaccinations<ScalarType>>().array().setConstant(vacc_full);
 
     auto& contacts       = model.parameters.get<mio::osecirvvs::ContactPatterns<ScalarType>>();
     auto& contact_matrix = contacts.get_cont_freq_mat();

cpp/examples/CMakeLists.txt

@@ -76,6 +76,10 @@ add_executable(ode_secirvvs_example ode_secirvvs.cpp)
 target_link_libraries(ode_secirvvs_example PRIVATE memilio ode_secirvvs)
 target_compile_options(ode_secirvvs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(ode_secirts_example ode_secirts.cpp)
+target_link_libraries(ode_secirts_example PRIVATE memilio ode_secirts)
+target_compile_options(ode_secirts_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
 add_executable(ode_secir_ageres_example ode_secir_ageres.cpp)
 target_link_libraries(ode_secir_ageres_example PRIVATE memilio ode_secir)
 target_compile_options(ode_secir_ageres_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

cpp/examples/ode_secirts.cpp

@@ -0,0 +1,154 @@
+/* 
+* Copyright (C) 2020-2024 MEmilio
+*
+* Authors: Henrik Zunker
+*
+* 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 "ode_secirts/analyze_result.h"
+#include "ode_secirts/model.h"
+#include "ode_secirts/parameters.h"
+#include "memilio/compartments/simulation.h"
+#include "memilio/utils/logging.h"
+
+int main()
+{
+    // This example demonstrates how to simulate a SECIRTS model.
+    // The SECIRTS model is an extension of the SECIRVVS model that includes waning and  temporary immunity.
+    // After the simulation, the aggregated size of the temporary immunity states are printed.
+    mio::set_log_level(mio::LogLevel::debug);
+
+    double t0   = 0;
+    double tmax = 100;
+    double dt   = 0.1;
+
+    mio::log_info("Simulating SECIRTS; t={} ... {} with dt = {}.", t0, tmax, dt);
+
+    mio::osecirts::Model<double> model(3);
+    auto nb_groups = model.parameters.get_num_groups();
+
+    for (mio::AgeGroup i = 0; i < nb_groups; i++) {
+        // population
+        model.populations[{i, mio::osecirts::InfectionState::ExposedNaive}]                                = 20;
+        model.populations[{i, mio::osecirts::InfectionState::ExposedImprovedImmunity}]                     = 20;
+        model.populations[{i, mio::osecirts::InfectionState::ExposedPartialImmunity}]                      = 20;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsNaive}]                     = 30;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsNaiveConfirmed}]            = 0;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsPartialImmunity}]           = 30;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}]  = 0;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsImprovedImmunity}]          = 30;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] = 0;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsNaive}]                       = 40;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsNaiveConfirmed}]              = 0;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsPartialImmunity}]             = 40;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsPartialImmunityConfirmed}]    = 0;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsImprovedImmunity}]            = 40;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsImprovedImmunityConfirmed}]   = 0;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedSevereNaive}]                         = 30;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedSevereImprovedImmunity}]              = 30;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedSeverePartialImmunity}]               = 30;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedCriticalNaive}]                       = 20;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedCriticalPartialImmunity}]             = 20;
+        model.populations[{i, mio::osecirts::InfectionState::InfectedCriticalImprovedImmunity}]            = 20;
+        model.populations[{i, mio::osecirts::InfectionState::SusceptibleNaive}]                            = 1000;
+        model.populations[{i, mio::osecirts::InfectionState::SusceptiblePartialImmunity}]                  = 1200;
+        model.populations[{i, mio::osecirts::InfectionState::SusceptibleImprovedImmunity}]                 = 1000;
+        model.populations[{i, mio::osecirts::InfectionState::TemporaryImmunPartialImmunity}]               = 60;
+        model.populations[{i, mio::osecirts::InfectionState::TemporaryImmunImprovedImmunity}]              = 70;
+        model.populations[{i, mio::osecirts::InfectionState::DeadNaive}]                                   = 0;
+        model.populations[{i, mio::osecirts::InfectionState::DeadPartialImmunity}]                         = 0;
+        model.populations[{i, mio::osecirts::InfectionState::DeadImprovedImmunity}]                        = 0;
+
+        // parameters
+        //times
+        model.parameters.get<mio::osecirts::TimeExposed<double>>()[i]                = 3.33;
+        model.parameters.get<mio::osecirts::TimeInfectedNoSymptoms<double>>()[i]     = 1.87;
+        model.parameters.get<mio::osecirts::TimeInfectedSymptoms<double>>()[i]       = 7;
+        model.parameters.get<mio::osecirts::TimeInfectedSevere<double>>()[i]         = 6;
+        model.parameters.get<mio::osecirts::TimeInfectedCritical<double>>()[i]       = 7;
+        model.parameters.get<mio::osecirts::TimeTemporaryImmunityPI<double>>()[i]    = 60;
+        model.parameters.get<mio::osecirts::TimeTemporaryImmunityII<double>>()[i]    = 60;
+        model.parameters.get<mio::osecirts::TimeWaningPartialImmunity<double>>()[i]  = 180;
+        model.parameters.get<mio::osecirts::TimeWaningImprovedImmunity<double>>()[i] = 180;
+
+        //probabilities
+        model.parameters.get<mio::osecirts::TransmissionProbabilityOnContact<double>>()[i]  = 0.15;
+        model.parameters.get<mio::osecirts::RelativeTransmissionNoSymptoms<double>>()[i]    = 0.5;
+        model.parameters.get<mio::osecirts::RiskOfInfectionFromSymptomatic<double>>()[i]    = 0.0;
+        model.parameters.get<mio::osecirts::MaxRiskOfInfectionFromSymptomatic<double>>()[i] = 0.4;
+        model.parameters.get<mio::osecirts::RecoveredPerInfectedNoSymptoms<double>>()[i]    = 0.2;
+        model.parameters.get<mio::osecirts::SeverePerInfectedSymptoms<double>>()[i]         = 0.1;
+        model.parameters.get<mio::osecirts::CriticalPerSevere<double>>()[i]                 = 0.1;
+        model.parameters.get<mio::osecirts::DeathsPerCritical<double>>()[i]                 = 0.1;
+
+        model.parameters.get<mio::osecirts::ReducExposedPartialImmunity<double>>()[i]                     = 0.8;
+        model.parameters.get<mio::osecirts::ReducExposedImprovedImmunity<double>>()[i]                    = 0.331;
+        model.parameters.get<mio::osecirts::ReducInfectedSymptomsPartialImmunity<double>>()[i]            = 0.65;
+        model.parameters.get<mio::osecirts::ReducInfectedSymptomsImprovedImmunity<double>>()[i]           = 0.243;
+        model.parameters.get<mio::osecirts::ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[i]  = 0.1;
+        model.parameters.get<mio::osecirts::ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[i] = 0.091;
+        model.parameters.get<mio::osecirts::ReducTimeInfectedMild<double>>()[i]                           = 0.9;
+    }
+
+    model.parameters.get<mio::osecirts::ICUCapacity<double>>()          = 100;
+    model.parameters.get<mio::osecirts::TestAndTraceCapacity<double>>() = 0.0143;
+    const size_t daily_vaccinations                                     = 10;
+    const size_t num_days                                               = 300;
+    model.parameters.get<mio::osecirts::DailyPartialVaccinations<double>>().resize(mio::SimulationDay(num_days));
+    model.parameters.get<mio::osecirts::DailyFullVaccinations<double>>().resize(mio::SimulationDay(num_days));
+    model.parameters.get<mio::osecirts::DailyBoosterVaccinations<double>>().resize(mio::SimulationDay(num_days));
+    for (size_t i = 0; i < num_days; ++i) {
+        for (mio::AgeGroup j = 0; j < nb_groups; ++j) {
+            auto num_vaccinations = static_cast<double>(i * daily_vaccinations);
+            model.parameters.get<mio::osecirts::DailyPartialVaccinations<double>>()[{j, mio::SimulationDay(i)}] =
+                num_vaccinations;
+            model.parameters.get<mio::osecirts::DailyFullVaccinations<double>>()[{j, mio::SimulationDay(i)}] =
+                num_vaccinations;
+            model.parameters.get<mio::osecirts::DailyBoosterVaccinations<double>>()[{j, mio::SimulationDay(i)}] =
+                num_vaccinations;
+        }
+    }
+
+    mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::osecirts::ContactPatterns<double>>();
+    const double cont_freq                  = 10;
+    const double fact                       = 1.0 / (double)(size_t)nb_groups;
+    contact_matrix[0] =
+        mio::ContactMatrix(Eigen::MatrixXd::Constant((size_t)nb_groups, (size_t)nb_groups, fact * cont_freq));
+    contact_matrix.add_damping(Eigen::MatrixXd::Constant((size_t)nb_groups, (size_t)nb_groups, 0.7),
+                               mio::SimulationTime(30.));
+
+    model.parameters.get<mio::osecirts::Seasonality<double>>() = 0.2;
+
+    model.apply_constraints();
+
+    mio::TimeSeries<double> result = simulate(t0, tmax, dt, model);
+
+    bool print_to_terminal = true;
+
+    if (print_to_terminal) {
+        auto result_interpolated = mio::interpolate_simulation_result(result);
+        for (auto t_indx = 0; t_indx < result_interpolated.get_num_time_points(); t_indx++) {
+            double timm_pi = 0.0;
+            double timm_ii = 0.0;
+            for (mio::AgeGroup i = 0; i < nb_groups; i++) {
+                timm_pi += result_interpolated.get_value(t_indx)[model.populations.get_flat_index(
+                    {i, mio::osecirts::InfectionState::TemporaryImmunPartialImmunity})];
+                timm_ii += result_interpolated.get_value(t_indx)[model.populations.get_flat_index(
+                    {i, mio::osecirts::InfectionState::TemporaryImmunImprovedImmunity})];
+            }
+            printf("t=%i, timm_pi=%f, timm_ii=%f\n", int(result_interpolated.get_time(t_indx)), timm_pi, timm_ii);
+        }
+    }
+}

cpp/examples/ode_secirvvs.cpp

@@ -66,16 +66,17 @@ int main()
 
     model.parameters.get<mio::osecirvvs::ICUCapacity<double>>()          = 100;
     model.parameters.get<mio::osecirvvs::TestAndTraceCapacity<double>>() = 0.0143;
-    const size_t daily_vaccinations = 10;
-    model.parameters.get<mio::osecirvvs::DailyFirstVaccination<double>>().resize(mio::SimulationDay((size_t)tmax + 1));
-    model.parameters.get<mio::osecirvvs::DailyFullVaccination<double>>().resize(mio::SimulationDay((size_t)tmax + 1));
+    const size_t daily_vaccinations                                      = 10;
+    model.parameters.get<mio::osecirvvs::DailyPartialVaccinations<double>>().resize(
+        mio::SimulationDay((size_t)tmax + 1));
+    model.parameters.get<mio::osecirvvs::DailyFullVaccinations<double>>().resize(mio::SimulationDay((size_t)tmax + 1));
     for (size_t i = 0; i < tmax + 1; ++i) {
         auto num_vaccinations = static_cast<double>(i * daily_vaccinations);
         model.parameters
-            .get<mio::osecirvvs::DailyFirstVaccination<double>>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
+            .get<mio::osecirvvs::DailyPartialVaccinations<double>>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
             num_vaccinations;
         model.parameters
-            .get<mio::osecirvvs::DailyFullVaccination<double>>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
+            .get<mio::osecirvvs::DailyFullVaccinations<double>>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
             num_vaccinations;
     }
     auto& contacts       = model.parameters.get<mio::osecirvvs::ContactPatterns<double>>();