905 add age group resolution to seir and sir model

cpp/examples/CMakeLists.txt

@@ -18,6 +18,10 @@ add_executable(ode_seir_example ode_seir.cpp)
 target_link_libraries(ode_seir_example PRIVATE memilio ode_seir)
 target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(ode_seir_ageres_example ode_seir_ageres.cpp)
+target_link_libraries(ode_seir_ageres_example PRIVATE memilio ode_seir)
+target_compile_options(ode_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
 add_executable(ode_sir_example ode_sir.cpp)
 target_link_libraries(ode_sir_example PRIVATE memilio ode_sir)
 target_compile_options(ode_sir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
@@ -30,6 +34,10 @@ add_executable(sde_sirs_example sde_sirs.cpp)
 target_link_libraries(sde_sirs_example PRIVATE memilio sde_sirs)
 target_compile_options(sde_sirs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(ode_sir_ageres_example ode_sir_ageres.cpp)
+target_link_libraries(ode_sir_ageres_example PRIVATE memilio ode_sir)
+target_compile_options(ode_sir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
 add_executable(seir_flows_example ode_seir_flows.cpp)
 target_link_libraries(seir_flows_example PRIVATE memilio ode_seir)
 target_compile_options(seir_flows_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

cpp/examples/graph.cpp

@@ -31,20 +31,30 @@ int main()
     const auto tmax = 10.;
     const auto dt   = 0.5; //time step of migration, daily migration every second step
 
-    mio::oseir::Model model;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Susceptible)}] = 10000;
+    mio::oseir::Model model(1);
+
+    // set population
+    model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 10000;
+
+    // set transition times
     model.parameters.set<mio::oseir::TimeExposed>(1);
-    model.parameters.get<mio::oseir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
     model.parameters.set<mio::oseir::TimeInfected>(1);
 
+    // set contact matrix
+    mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::oseir::ContactPatterns>().get_cont_freq_mat();
+    contact_matrix[0].get_baseline().setConstant(2.7);
+
     //two mostly identical groups
     auto model_group1 = model;
     auto model_group2 = model;
+
     //some contact restrictions in group 1
-    model_group1.parameters.get<mio::oseir::ContactPatterns>().add_damping(0.5, mio::SimulationTime(5));
+    mio::ContactMatrixGroup& contact_matrix1 = model_group1.parameters.get<mio::oseir::ContactPatterns>().get_cont_freq_mat();
+    contact_matrix1[0].add_damping(0.5, mio::SimulationTime(5));
+
     //infection starts in group 1
-    model_group1.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Susceptible)}] = 9990;
-    model_group1.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}]     = 10;
+    model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 9990;
+    model_group1.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}]     = 10;
 
     mio::Graph<mio::SimulationNode<mio::Simulation<mio::oseir::Model>>, mio::MigrationEdge> g;
     g.add_node(1001, model_group1, t0);

cpp/examples/ode_seir.cpp

@@ -23,6 +23,8 @@
 #include "memilio/compartments/simulation.h"
 #include "memilio/utils/logging.h"
 
+#include "memilio/utils/time_series.h"
+
 int main()
 {
     mio::set_log_level(mio::LogLevel::debug);
@@ -33,26 +35,25 @@ int main()
 
     mio::log_info("Simulating ODE SEIR; t={} ... {} with dt = {}.", t0, tmax, dt);
 
-    mio::oseir::Model model;
+    mio::oseir::Model model(1);
 
-    double total_population                                                                            = 1061000;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}]   = 10000;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}]  = 1000;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}] = 1000;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Susceptible)}] =
-        total_population -
-        model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}] -
-        model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}] -
-        model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}];
+    double total_population                                                      = 10000;
+    model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}]   = 100;
+    model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}]  = 100;
+    model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100;
+    model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] =
+        total_population - model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] -
+        model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}] -
+        model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}];
 
-    model.parameters.set<mio::oseir::TransmissionProbabilityOnContact>(0.04);
     model.parameters.set<mio::oseir::TimeExposed>(5.2);
-    model.parameters.set<mio::oseir::TimeInfected>(2);
-
-    model.parameters.get<mio::oseir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
+    model.parameters.set<mio::oseir::TimeInfected>(6);
+    model.parameters.set<mio::oseir::TransmissionProbabilityOnContact>(0.1);
 
-    // contacts increase by 100% after 12.5 days
-    model.parameters.get<mio::oseir::ContactPatterns>().add_damping(-1., mio::SimulationTime(12.5));
+    mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::oseir::ContactPatterns>();
+    contact_matrix[0].get_baseline().setConstant(2.7);
+    contact_matrix[0].add_damping(0.7, mio::SimulationTime(30.));
+
 
     model.check_constraints();
 

cpp/examples/ode_seir_ageres.cpp

@@ -0,0 +1,62 @@
+#include "ode_seir/model.h"
+#include "ode_seir/infection_state.h"
+#include "ode_seir/parameters.h"
+#include "memilio/compartments/simulation.h"
+#include "memilio/utils/logging.h"
+
+int main()
+{
+    mio::set_log_level(mio::LogLevel::debug);
+
+    double t0   = 0;
+    double tmax = 50;
+    double dt   = 0.001;
+
+    mio::log_info("Simulating SEIR; t={} ... {} with dt = {}.", t0, tmax, dt);
+
+    double cont_freq = 10;
+
+    double nb_total_t0 = 10000, nb_exp_t0 = 100, nb_inf_t0 = 50, nb_rec_t0 = 10;
+    const size_t num_groups = 3;
+
+    mio::oseir::Model model(num_groups);
+    double fact = 1.0 / num_groups;
+
+    auto& params = model.parameters;
+
+    for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); i++) {
+        model.populations[{i, mio::oseir::InfectionState::Exposed}]   = fact * nb_exp_t0;
+        model.populations[{i, mio::oseir::InfectionState::Infected}]  = fact * nb_inf_t0;
+        model.populations[{i, mio::oseir::InfectionState::Recovered}] = fact * nb_rec_t0;
+        model.populations.set_difference_from_group_total<mio::AgeGroup>({i, mio::oseir::InfectionState::Susceptible},
+                                                                         fact * nb_total_t0);
+
+        model.parameters.get<mio::oseir::TimeExposed>()[i]                      = 5.2;
+        model.parameters.get<mio::oseir::TimeInfected>()[i]                     = 6;
+        model.parameters.get<mio::oseir::TransmissionProbabilityOnContact>()[i] = 0.04;
+    }
+
+    mio::ContactMatrixGroup& contact_matrix = params.get<mio::oseir::ContactPatterns>();
+    contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_groups, num_groups, fact * cont_freq));
+    contact_matrix.add_damping(Eigen::MatrixXd::Constant(num_groups, num_groups, 0.7), mio::SimulationTime(30.));
+
+    model.apply_constraints();
+
+    mio::TimeSeries<double> seir = simulate(t0, tmax, dt, model);
+
+    std::vector<std::string> vars = {"S", "E", "I", "R"};
+    printf("Number of time points :%d\n", static_cast<int>(seir.get_num_time_points()));
+    printf("People in\n");
+
+    for (size_t k = 0; k < (size_t)mio::oseir::InfectionState::Count; k++) {
+        double dummy = 0;
+
+        for (size_t i = 0; i < (size_t)params.get_num_groups(); i++) {
+            printf("\t %s[%d]: %.0f", vars[k].c_str(), (int)i,
+                   seir.get_last_value()[k + (size_t)mio::oseir::InfectionState::Count * (int)i]);
+            dummy += seir.get_last_value()[k + (size_t)mio::oseir::InfectionState::Count * (int)i];
+        }
+
+        printf("\t %s_otal: %.0f\n", vars[k].c_str(), dummy);
+    }
+}

cpp/examples/ode_seir_flows.cpp

@@ -38,25 +38,24 @@ int main()
 
     mio::log_info("Simulating SEIR; t={} ... {} with dt = {}.", t0, tmax, dt);
 
-    mio::oseir::Model model;
+    mio::oseir::Model model(1);
+
+    constexpr double total_population                                            = 10000;
+    model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}]   = 100;
+    model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Infected}]  = 100;
+    model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Recovered}] = 100;
+    model.populations.set_difference_from_group_total<mio::AgeGroup>(
+        {mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}, total_population);
 
-    double total_population                                                                            = 10000;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}]   = 100;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}]  = 100;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}] = 100;
-    model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Susceptible)}] =
-        total_population -
-        model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Exposed)}] -
-        model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Infected)}] -
-        model.populations[{mio::Index<mio::oseir::InfectionState>(mio::oseir::InfectionState::Recovered)}];
-    // suscetible now set with every other update
-    // params.nb_sus_t0   = params.nb_total_t0 - params.nb_exp_t0 - params.nb_inf_t0 - params.nb_rec_t0;
     model.parameters.set<mio::oseir::TimeExposed>(5.2);
     model.parameters.set<mio::oseir::TimeInfected>(6);
     model.parameters.set<mio::oseir::TransmissionProbabilityOnContact>(0.04);
-    model.parameters.get<mio::oseir::ContactPatterns>().get_baseline()(0, 0) = 10;
+
+    mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::oseir::ContactPatterns>().get_cont_freq_mat();
+    contact_matrix[0].get_baseline().setConstant(10);
 
     model.check_constraints();
+
     auto seir = simulate_flows(t0, tmax, dt, model);
 
     printf("Compartments: \n");

cpp/examples/ode_sir.cpp

@@ -24,6 +24,8 @@
 #include "ode_sir/infection_state.h"
 #include "ode_sir/model.h"
 #include "ode_sir/parameters.h"
+#include <fstream>
+#include <stdio.h>
 
 int main()
 {
@@ -37,24 +39,23 @@ int main()
 
     mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt);
 
-    mio::osir::Model model;
+    mio::osir::Model model(1);
 
-    model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Infected)}]  = 1000;
-    model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Recovered)}] = 1000;
-    model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Susceptible)}] =
-        total_population -
-        model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Infected)}] -
-        model.populations[{mio::Index<mio::osir::InfectionState>(mio::osir::InfectionState::Recovered)}];
+    model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Infected}]  = 1000;
+    model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Recovered}] = 1000;
+    model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Susceptible}] =
+        total_population - model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Infected}] -
+        model.populations[{mio::AgeGroup(0), mio::osir::InfectionState::Recovered}];
     model.parameters.set<mio::osir::TimeInfected>(2);
-    model.parameters.set<mio::osir::TransmissionProbabilityOnContact>(0.04);
-    model.parameters.get<mio::osir::ContactPatterns>().get_baseline()(0, 0) = 2.7;
-    model.parameters.get<mio::osir::ContactPatterns>().add_damping(0.6, mio::SimulationTime(12.5));
-
-    auto integrator = std::make_shared<mio::EulerIntegratorCore>();
+    model.parameters.set<mio::osir::TransmissionProbabilityOnContact>(0.5);
 
+    mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::osir::ContactPatterns>().get_cont_freq_mat();
+    contact_matrix[0].get_baseline().setConstant(2.7);
+    contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5));
     model.check_constraints();
 
-    auto sir = simulate(t0, tmax, dt, model, integrator);
+    auto integrator = std::make_shared<mio::EulerIntegratorCore>();
+    auto sir        = simulate(t0, tmax, dt, model, integrator);
 
     bool print_to_terminal = true;
 

cpp/examples/ode_sir_ageres.cpp

@@ -0,0 +1,61 @@
+#include "memilio/epidemiology/age_group.h"
+#include "memilio/utils/time_series.h"
+#include "memilio/utils/logging.h"
+#include "memilio/compartments/simulation.h"
+#include "ode_sir/model.h"
+
+int main()
+{
+    mio::set_log_level(mio::LogLevel::debug);
+
+    double t0   = 0;
+    double tmax = 50;
+    double dt   = 0.1;
+
+    mio::log_info("Simulating SIR; t={} ... {} with dt = {}.", t0, tmax, dt);
+
+    double cont_freq = 10; // see Polymod study
+
+    double nb_total_t0 = 10000, nb_inf_t0 = 50, nb_rec_t0 = 10;
+
+    const size_t num_groups = 3;
+    mio::osir::Model model(num_groups);
+
+    double fact = 1.0 / num_groups;
+
+    auto& params = model.parameters;
+
+    for (auto i = mio::AgeGroup(0); i < mio::AgeGroup(num_groups); i++) {
+        model.populations[{i, mio::osir::InfectionState::Infected}]  = fact * nb_inf_t0;
+        model.populations[{i, mio::osir::InfectionState::Recovered}] = fact * nb_rec_t0;
+        model.populations.set_difference_from_group_total<mio::AgeGroup>({i, mio::osir::InfectionState::Susceptible},
+                                                                         fact * nb_total_t0);
+
+        model.parameters.get<mio::osir::TimeInfected>()[i]                     = 2.0;
+        model.parameters.get<mio::osir::TransmissionProbabilityOnContact>()[i] = 0.3;
+    }
+
+    mio::ContactMatrixGroup& contact_matrix = params.get<mio::osir::ContactPatterns>();
+    contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(num_groups, num_groups, fact * cont_freq));
+    contact_matrix.add_damping(Eigen::MatrixXd::Constant(num_groups, num_groups, 0.7), mio::SimulationTime(30.));
+
+    model.apply_constraints();
+
+    mio::TimeSeries<double> sir = simulate(t0, tmax, dt, model);
+
+    std::vector<std::string> vars = {"S", "I", "R"};
+    printf("Number of time points :%d\n", static_cast<int>(sir.get_num_time_points()));
+    printf("People in\n");
+
+    for (size_t k = 0; k < (size_t)mio::osir::InfectionState::Count; k++) {
+        double dummy = 0;
+
+        for (size_t i = 0; i < (size_t)params.get_num_groups(); i++) {
+            printf("\t %s[%d]: %.0f", vars[k].c_str(), (int)i,
+                   sir.get_last_value()[k + (size_t)mio::osir::InfectionState::Count * (int)i]);
+            dummy += sir.get_last_value()[k + (size_t)mio::osir::InfectionState::Count * (int)i];
+        }
+
+        printf("\t %s_otal: %.0f\n", vars[k].c_str(), dummy);
+    }
+}