300 python bindings for new vaccination model

cpp/examples/ode_secirvvs.cpp

@@ -66,11 +66,16 @@ int main()
 
     model.parameters.get<mio::osecirvvs::ICUCapacity>()          = 100;
     model.parameters.get<mio::osecirvvs::TestAndTraceCapacity>() = 0.0143;
-    model.parameters.get<mio::osecirvvs::DailyFirstVaccination>().resize(mio::SimulationDay(size_t(1000)));
-    model.parameters.get<mio::osecirvvs::DailyFirstVaccination>().array().setConstant(5);
-    model.parameters.get<mio::osecirvvs::DailyFullVaccination>().resize(mio::SimulationDay(size_t(1000)));
-    model.parameters.get<mio::osecirvvs::DailyFullVaccination>().array().setConstant(3);
-
+    const size_t daily_vaccinations = 10;
+    model.parameters.get<mio::osecirvvs::DailyFirstVaccination>().resize(mio::SimulationDay((size_t)tmax + 1));
+    model.parameters.get<mio::osecirvvs::DailyFullVaccination>().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>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
+            num_vaccinations;
+        model.parameters.get<mio::osecirvvs::DailyFullVaccination>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
+            num_vaccinations;
+    }
     auto& contacts       = model.parameters.get<mio::osecirvvs::ContactPatterns>();
     auto& contact_matrix = contacts.get_cont_freq_mat();
     contact_matrix[0].get_baseline().setConstant(0.5);

cpp/models/ode_secirvvs/parameters_io.h

@@ -1225,6 +1225,24 @@ export_input_data_county_timeseries(std::vector<Model>&& model, const std::strin
 
     return success();
 }
+#else
+template <class Model>
+IOResult<void> export_input_data_county_timeseries(std::vector<Model>&, const std::string&, std::vector<int> const&,
+                                                   Date, const std::vector<double>&, double, int, const std::string&,
+                                                   const std::string&, const std::string&)
+{
+    mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data.");
+    return success();
+}
+
+template <class Model>
+IOResult<void> export_input_data_county_timeseries(std::vector<Model>&&, const std::string&, std::vector<int> const&,
+                                                   Date, const std::vector<double>&, double, int, const std::string&,
+                                                   const std::string&, const std::string&, bool, const std::string&)
+{
+    mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data.");
+    return success();
+}
 
 #endif //MEMILIO_HAS_HDF5
 

pycode/examples/simulation/2020_npis_sarscov2_wildtype_germany.py

@@ -177,8 +177,8 @@ def set_contact_matrices(self, model):
             minimum_file = os.path.join(
                 self.data_dir, "contacts", "minimum_" + location + ".txt")
             contact_matrices[i] = mio.ContactMatrix(
-                mio.secir.read_mobility_plain(baseline_file),
-                mio.secir.read_mobility_plain(minimum_file)
+                mio.read_mobility_plain(baseline_file),
+                mio.read_mobility_plain(minimum_file)
             )
         model.parameters.ContactPatterns.cont_freq_mat = contact_matrices
 
@@ -203,7 +203,7 @@ def set_npis(self, params, end_date):
 
         typ_home = Intervention.Home.value
         typ_school = Intervention.SchoolClosure.value
-        typ_home = Intervention.HomeOffice.value
+        typ_homeoffice = Intervention.HomeOffice.value
         typ_gathering = Intervention.GatheringBanFacilitiesClosure.value
         typ_distance = Intervention.PhysicalDistanceAndMasks.value
         typ_senior = Intervention.SeniorAwareness.value
@@ -231,7 +231,7 @@ def school_closure(t, min, max):
 
         def home_office(t, min, max):
             return damping_helper(
-                t, min, max, lvl_main, typ_home, [loc_work])
+                t, min, max, lvl_main, typ_homeoffice, [loc_work])
 
         def social_events(t, min, max):
             return damping_helper(

pycode/examples/simulation/osecirvvs_simple.py

@@ -0,0 +1,165 @@
+#############################################################################
+# Copyright (C) 2020-2024 MEmilio
+#
+# Authors: Martin J. Kuehn, Maximilian Betz
+#
+# 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.
+#############################################################################
+import argparse
+import os
+from datetime import date, datetime
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+from memilio.simulation import AgeGroup, Damping, SimulationDay
+from memilio.simulation.osecirvvs import InfectionState
+from memilio.simulation.osecirvvs import Model, simulate
+
+
+def run_secirvvs_simulation(show_plot=True):
+    """
+    Runs the c++ SECIRVVS model using a single age group 
+    and plots the results
+    """
+
+    t0 = 0
+    tmax = 30  # number of days to simulate
+    dt = 0.1
+    num_groups = 1
+
+    # Initialize Parameters
+    model = Model(num_groups)
+
+    # set parameters
+    for i in range(num_groups):
+        # Initial number of peaople in each compartment
+        model.populations[AgeGroup(i), InfectionState.ExposedNaive] = 10
+        model.populations[AgeGroup(
+            i), InfectionState.ExposedImprovedImmunity] = 11
+        model.populations[AgeGroup(
+            i), InfectionState.ExposedPartialImmunity] = 12
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedNoSymptomsNaive] = 13
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedNoSymptomsNaiveConfirmed] = 13
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedNoSymptomsPartialImmunity] = 14
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedNoSymptomsPartialImmunityConfirmed] = 14
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedNoSymptomsImprovedImmunity] = 15
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedNoSymptomsImprovedImmunityConfirmed] = 15
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedSymptomsNaive] = 5
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedSymptomsNaiveConfirmed] = 5
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedSymptomsPartialImmunity] = 6
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedSymptomsPartialImmunityConfirmed] = 6
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedSymptomsImprovedImmunity] = 7
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedSymptomsImprovedImmunityConfirmed] = 7
+        model.populations[AgeGroup(i), InfectionState.InfectedSevereNaive] = 8
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedSevereImprovedImmunity] = 1
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedSeverePartialImmunity] = 2
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedCriticalNaive] = 3
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedCriticalPartialImmunity] = 4
+        model.populations[AgeGroup(
+            i), InfectionState.InfectedCriticalImprovedImmunity] = 5
+        model.populations[AgeGroup(
+            i), InfectionState.SusceptibleImprovedImmunity] = 6
+        model.populations[AgeGroup(
+            i), InfectionState.SusceptiblePartialImmunity] = 7
+        model.populations[AgeGroup(i), InfectionState.DeadNaive] = 0
+        model.populations[AgeGroup(i), InfectionState.DeadPartialImmunity] = 0
+        model.populations[AgeGroup(i), InfectionState.DeadImprovedImmunity] = 0
+        model.populations.set_difference_from_group_total_AgeGroup(
+            (AgeGroup(i), InfectionState.SusceptibleNaive), 1000)
+
+    model.parameters.ICUCapacity.value = 100
+    model.parameters.TestAndTraceCapacity.value = 0.0143
+    model.parameters.DailyFirstVaccination.resize_SimulationDay(
+        SimulationDay(tmax + 1))
+    model.parameters.DailyFullVaccination.resize_SimulationDay(
+        SimulationDay(tmax + 1))
+    daily_vaccinations = 10
+    for i, num_vaccinations in enumerate(range(0, daily_vaccinations * (tmax + 1), daily_vaccinations)):
+        model.parameters.DailyFirstVaccination[AgeGroup(
+            0), SimulationDay(i)] = num_vaccinations
+        model.parameters.DailyFullVaccination[AgeGroup(
+            0), SimulationDay(i)] = num_vaccinations
+
+    # contact patterns
+    baseline = np.ones((num_groups, num_groups)) * 0.5
+    np.fill_diagonal(baseline, 5.0)
+    model.parameters.ContactPatterns.cont_freq_mat[0].baseline = baseline
+    model.parameters.ContactPatterns.cont_freq_mat.add_damping(Damping(
+        coeffs=np.ones((num_groups, num_groups)) * 0.3, t=5.0, level=0, type=0))
+
+    # times
+    model.parameters.TimeInfectedSymptoms[AgeGroup(0)] = 7
+    model.parameters.TimeInfectedSevere[AgeGroup(0)] = 6
+    model.parameters.TimeInfectedCritical[AgeGroup(0)] = 7
+
+    # probabilities
+    model.parameters.TransmissionProbabilityOnContact[AgeGroup(0)] = 0.15
+    model.parameters.RelativeTransmissionNoSymptoms[AgeGroup(0)] = 0.5
+    # The precise value between Risk* (situation under control) and MaxRisk* (situation not under control)
+    # depends on incidence and test and trace capacity
+    model.parameters.RiskOfInfectionFromSymptomatic[AgeGroup(0)] = 0.0
+    model.parameters.MaxRiskOfInfectionFromSymptomatic[AgeGroup(0)] = 0.4
+    model.parameters.RecoveredPerInfectedNoSymptoms[AgeGroup(0)] = 0.2
+    model.parameters.SeverePerInfectedSymptoms[AgeGroup(0)] = 0.1
+    model.parameters.CriticalPerSevere[AgeGroup(0)] = 0.1
+    model.parameters.DeathsPerCritical[AgeGroup(0)] = 0.1
+
+    model.parameters.ReducExposedPartialImmunity[AgeGroup(0)] = 0.8
+    model.parameters.ReducExposedImprovedImmunity[AgeGroup(0)] = 0.331
+    model.parameters.ReducInfectedSymptomsPartialImmunity[AgeGroup(0)] = 0.65
+    model.parameters.ReducInfectedSymptomsImprovedImmunity[AgeGroup(0)] = 0.243
+    model.parameters.ReducInfectedSevereCriticalDeadPartialImmunity[AgeGroup(
+        0)] = 0.1
+    model.parameters.ReducInfectedSevereCriticalDeadImprovedImmunity[AgeGroup(
+        0)] = 0.091
+    model.parameters.ReducTimeInfectedMild[AgeGroup(0)] = 0.9
+
+    model.parameters.Seasonality.value = 0.2
+
+    model.apply_constraints()
+
+    # Run Simulation
+    result = simulate(t0, tmax, dt, model)
+
+    # # interpolate results
+    # result = interpolate_simulation_result(result)
+
+    print(result.get_last_value())
+
+
+if __name__ == "__main__":
+    arg_parser = argparse.ArgumentParser(
+        'osecirvvs_simple',
+        description='Simple example demonstrating the setup and simulation of the SECIRVVS model with a single age group.')
+    args = arg_parser.parse_args()
+    run_secirvvs_simulation(**args.__dict__)

pycode/memilio-simulation/CMakeLists.txt

@@ -33,13 +33,6 @@ endif()
 add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/../../cpp ${CMAKE_CURRENT_BINARY_DIR}/cpp EXCLUDE_FROM_ALL)
 
 # build python extensions
-pybind11_add_module(_simulation_secir MODULE
-    memilio/simulation/secir.cpp
-)
-target_link_libraries(_simulation_secir PRIVATE memilio ode_secir)
-target_include_directories(_simulation_secir PRIVATE memilio/simulation)
-install(TARGETS _simulation_secir LIBRARY DESTINATION memilio)
-
 pybind11_add_module(_simulation_abm MODULE
     memilio/simulation/abm.cpp
 )
@@ -62,17 +55,30 @@ target_link_libraries(_simulation PRIVATE memilio)
 target_include_directories(_simulation PRIVATE memilio/simulation)
 install(TARGETS _simulation LIBRARY DESTINATION memilio)
 
+pybind11_add_module(_simulation_osir MODULE
+    memilio/simulation/osir.cpp
+)
+target_link_libraries(_simulation_osir PRIVATE memilio ode_sir)
+target_include_directories(_simulation_osir PRIVATE memilio/simulation)
+install(TARGETS _simulation_osir LIBRARY DESTINATION memilio)
+
 pybind11_add_module(_simulation_oseir MODULE
     memilio/simulation/oseir.cpp
 )
 target_link_libraries(_simulation_oseir PRIVATE memilio ode_seir)
 target_include_directories(_simulation_oseir PRIVATE memilio/simulation)
 install(TARGETS _simulation_oseir LIBRARY DESTINATION memilio)
 
-pybind11_add_module(_simulation_osir MODULE
-    memilio/simulation/osir.cpp
+pybind11_add_module(_simulation_secir MODULE
+    memilio/simulation/secir.cpp
 )
-target_link_libraries(_simulation_osir PRIVATE memilio ode_sir)
-target_include_directories(_simulation_osir PRIVATE memilio/simulation)
-install(TARGETS _simulation_osir LIBRARY DESTINATION memilio)
+target_link_libraries(_simulation_secir PRIVATE memilio ode_secir)
+target_include_directories(_simulation_secir PRIVATE memilio/simulation)
+install(TARGETS _simulation_secir LIBRARY DESTINATION memilio)
 
+pybind11_add_module(_simulation_osecirvvs MODULE
+    memilio/simulation/osecirvvs.cpp
+)
+target_link_libraries(_simulation_osecirvvs PRIVATE memilio ode_secirvvs)
+target_include_directories(_simulation_osecirvvs PRIVATE memilio/simulation)
+install(TARGETS _simulation_osecirvvs LIBRARY DESTINATION memilio)

pycode/memilio-simulation/memilio/simulation/io/mobility_io.h

@@ -34,7 +34,7 @@ void bind_write_graph(pybind11::module_& m)
 {
     m.def("write_graph", [&](const mio::Graph<Model, mio::MigrationParameters>& graph, const std::string& directory) {
         int ioflags = mio::IOF_None;
-        mio::write_graph<Model>(graph, directory, ioflags);
+        auto ioresult = mio::write_graph<Model>(graph, directory, ioflags);
     });
 }
 

pycode/memilio-simulation/memilio/simulation/io/result_io.h

@@ -39,7 +39,7 @@ void bind_save_results(pybind11::module_& m)
               const std::vector<std::vector<Model>>& ensemble_params, const std::vector<int>& county_ids,
               const std::string& result_dir, bool save_single_runs, bool save_percentiles) {
               boost::filesystem::path dir(result_dir);
-              mio::save_results<Model>(ensemble_results, ensemble_params, county_ids, dir, save_single_runs,
+              auto ioresult = mio::save_results<Model>(ensemble_results, ensemble_params, county_ids, dir, save_single_runs,
                                        save_percentiles);
               return NULL;
           });