1003 implement model in which the mobility is integrated into the ODEs

cpp/CMakeLists.txt

@@ -161,6 +161,8 @@ if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/ode_seir)
     add_subdirectory(models/ode_seair)
     add_subdirectory(models/ode_sir)
+    add_subdirectory(models/ode_metapop)
+    add_subdirectory(models/ode_metapop_wang)
     add_subdirectory(models/sde_sir)
     add_subdirectory(models/sde_sirs)
     add_subdirectory(models/sde_seirvv)

cpp/examples/CMakeLists.txt

@@ -30,6 +30,41 @@ add_executable(sde_sir_example sde_sir.cpp)
 target_link_libraries(sde_sir_example PRIVATE memilio sde_sir)
 target_compile_options(sde_sir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(ode_metapop_wang_nrw examples_thesis/ode_metapop_wang_nrw.cpp)
+target_link_libraries(ode_metapop_wang_nrw PRIVATE memilio ode_metapop_wang)
+target_compile_options(ode_metapop_wang_nrw PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+if (MEMILIO_ENABLE_OPENMP)
+    add_executable(ode_metapop_timing examples_thesis/ode_metapop_timing.cpp)
+    target_link_libraries(ode_metapop_timing PRIVATE memilio ode_metapop likwid stdc++)
+    target_compile_options(ode_metapop_timing PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+    target_compile_definitions(ode_metapop_timing PRIVATE "-DLIKWID_PERFMON")
+endif()
+
+add_executable(ode_metapop_steps examples_thesis/ode_metapop_steps.cpp)
+target_link_libraries(ode_metapop_steps PRIVATE memilio ode_metapop)
+target_compile_options(ode_metapop_steps PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(ode_metapop_seir_example ode_seir_metapop.cpp)
+target_link_libraries(ode_metapop_seir_example PRIVATE memilio ode_metapop)
+target_compile_options(ode_metapop_seir_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(basic_reproduction_number_modela examples_thesis/basic_reproduction_number_modela.cpp)
+target_link_libraries(basic_reproduction_number_modela PRIVATE memilio ode_seir)
+target_compile_options(basic_reproduction_number_modela PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(basic_reproduction_number_modelb examples_thesis/basic_reproduction_number_modelb.cpp)
+target_link_libraries(basic_reproduction_number_modelb PRIVATE memilio ode_metapop_wang)
+target_compile_options(basic_reproduction_number_modelb PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(basic_reproduction_number_modelc examples_thesis/basic_reproduction_number_modelc.cpp)
+target_link_libraries(basic_reproduction_number_modelc PRIVATE memilio ode_metapop)
+target_compile_options(basic_reproduction_number_modelc PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(ode_metapop_nrw examples_thesis/ode_metapop_nrw.cpp)
+target_link_libraries(ode_metapop_nrw PRIVATE memilio ode_metapop)
+target_compile_options(ode_metapop_nrw PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
 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})
@@ -92,6 +127,21 @@ add_executable(ode_secir_graph_example ode_secir_graph.cpp)
 target_link_libraries(ode_secir_graph_example PRIVATE memilio ode_secir)
 target_compile_options(ode_secir_graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(graph_nrw examples_thesis/graph_nrw.cpp)
+target_link_libraries(graph_nrw PRIVATE memilio ode_seir)
+target_compile_options(graph_nrw PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(graph_steps examples_thesis/graph_steps.cpp)
+target_link_libraries(graph_steps PRIVATE memilio ode_seir)
+target_compile_options(graph_steps PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+if (MEMILIO_ENABLE_OPENMP)
+    add_executable(graph_timing examples_thesis/graph_timing.cpp)
+    target_link_libraries(graph_timing PRIVATE memilio ode_seir likwid stdc++)
+    target_compile_options(graph_timing PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+    target_compile_definitions(graph_timing PRIVATE "-DLIKWID_PERFMON")
+endif()
+
 add_executable(graph_stochastic_mobility_example graph_stochastic_mobility.cpp)
 target_link_libraries(graph_stochastic_mobility_example PRIVATE memilio ode_secir)
 target_compile_options(graph_stochastic_mobility_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

cpp/examples/examples_thesis/basic_reproduction_number_modela.cpp

@@ -0,0 +1,64 @@
+#include "models/ode_seir/model.h"
+
+#include "memilio/math/euler.h"
+#include "memilio/compartments/simulation.h"
+#include "memilio/utils/custom_index_array.h"
+
+Eigen::MatrixXd get_contact_matrix()
+{
+    Eigen::MatrixXd contact_matrix_eigen(6, 6);
+    contact_matrix_eigen << 3.9547, 1.1002, 2.9472, 2.05, 0.3733, 0.0445, 0.3327, 3.5892, 1.236, 1.9208, 0.2681, 0.0161,
+        0.246, 0.7124, 5.6518, 3.2939, 0.2043, 0.0109, 0.1742, 0.8897, 3.3124, 4.5406, 0.4262, 0.0214, 0.0458, 0.1939,
+        0.5782, 1.3825, 1.473, 0.0704, 0.1083, 0.1448, 0.4728, 0.9767, 0.6266, 0.1724;
+
+    return contact_matrix_eigen;
+}
+
+const ScalarType TimeExposed[]                      = {3.335, 3.335, 3.335, 3.335, 3.335, 3.335};
+const ScalarType TimeInfected[]                     = {8.0096875, 8.0096875, 8.2182, 8.1158, 8.033, 7.985};
+const ScalarType TransmissionProbabilityOnContact[] = {0.03, 0.06, 0.06, 0.06, 0.09, 0.175};
+
+void calculate_basic_reproduction_number(size_t number_regions, ScalarType tmax)
+{
+    mio::set_log_level(mio::LogLevel::off);
+    ScalarType t0                = 0.;
+    ScalarType dt                = 0.1;
+    ScalarType number_age_groups = 6;
+
+    mio::oseir::Model<ScalarType> model(number_age_groups);
+    auto& population = model.populations;
+
+    for (size_t j = 0; j < number_age_groups; j++) {
+
+        population[{mio::AgeGroup(j), mio::oseir::InfectionState::Susceptible}] = number_regions * 10000;
+    }
+
+    mio::ContactMatrixGroup& contact_matrix =
+        model.parameters.template get<mio::oseir::ContactPatterns<>>().get_cont_freq_mat();
+    contact_matrix[0].get_baseline() = get_contact_matrix();
+
+    for (size_t j = 0; j < number_age_groups; j++) {
+        model.parameters.template get<mio::oseir::TimeExposed<>>()[mio::AgeGroup(j)]  = TimeExposed[j];
+        model.parameters.template get<mio::oseir::TimeInfected<>>()[mio::AgeGroup(j)] = TimeInfected[j];
+        model.parameters.template get<mio::oseir::TransmissionProbabilityOnContact<>>()[mio::AgeGroup(j)] =
+            TransmissionProbabilityOnContact[j];
+    }
+
+    std::shared_ptr<mio::IntegratorCore<ScalarType>> integrator = std::make_shared<mio::EulerIntegratorCore<>>();
+
+    auto result = simulate(t0, tmax, dt, model, integrator);
+
+    auto basic_reproduction_number = model.get_reproduction_number(t0, result).value();
+    std::cout << "\"Model A\": " << basic_reproduction_number << ", " << std::endl;
+}
+
+int main()
+{
+    const ScalarType tmax = 1.;
+    size_t num_regions    = 1;
+
+    std::cout << "{ \"Regions\": " << num_regions << ", " << std::endl;
+
+    calculate_basic_reproduction_number(num_regions, tmax);
+    return 0;
+}

cpp/examples/examples_thesis/basic_reproduction_number_modelb.cpp

@@ -0,0 +1,78 @@
+
+#include "models/ode_metapop_wang/model.h"
+
+#include "memilio/math/euler.h"
+#include "memilio/compartments/simulation.h"
+#include "memilio/utils/custom_index_array.h"
+
+Eigen::MatrixXd get_contact_matrix()
+{
+    Eigen::MatrixXd contact_matrix_eigen(6, 6);
+    contact_matrix_eigen << 3.9547, 1.1002, 2.9472, 2.05, 0.3733, 0.0445, 0.3327, 3.5892, 1.236, 1.9208, 0.2681, 0.0161,
+        0.246, 0.7124, 5.6518, 3.2939, 0.2043, 0.0109, 0.1742, 0.8897, 3.3124, 4.5406, 0.4262, 0.0214, 0.0458, 0.1939,
+        0.5782, 1.3825, 1.473, 0.0704, 0.1083, 0.1448, 0.4728, 0.9767, 0.6266, 0.1724;
+
+    return contact_matrix_eigen;
+}
+
+const ScalarType TimeExposed[]                      = {3.335, 3.335, 3.335, 3.335, 3.335, 3.335};
+const ScalarType TimeInfected[]                     = {8.0096875, 8.0096875, 8.2182, 8.1158, 8.033, 7.985};
+const ScalarType TransmissionProbabilityOnContact[] = {0.03, 0.06, 0.06, 0.06, 0.09, 0.175};
+
+void calculate_basic_reproduction_number(size_t number_regions, ScalarType tmax)
+{
+    mio::set_log_level(mio::LogLevel::off);
+    ScalarType t0                = 0.;
+    ScalarType dt                = 0.1;
+    ScalarType number_age_groups = 6;
+
+    mio::oseirmetapopwang::Model<ScalarType> model(number_regions, number_age_groups);
+    auto& population = model.populations;
+
+    for (size_t j = 0; j < number_age_groups; j++) {
+        for (size_t i = 0; i < number_regions; i++) {
+            population[{mio::oseirmetapopwang::Region(i), mio::AgeGroup(j),
+                        mio::oseirmetapopwang::InfectionState::Susceptible}] = 10000;
+        }
+    }
+
+    double fraction_commuter = 1. / (2 * number_regions);
+    Eigen::MatrixXd mobility_data_commuter =
+        Eigen::MatrixXd::Constant(number_regions, number_regions, fraction_commuter) -
+        fraction_commuter *
+            Eigen::MatrixXd::Identity(number_regions, number_regions); // Ensure that the diagonal is zero
+    for (size_t county_idx_i = 0; county_idx_i < number_regions; ++county_idx_i) {
+        mobility_data_commuter(county_idx_i, county_idx_i) = 1 - mobility_data_commuter.rowwise().sum()(county_idx_i);
+    }
+    model.parameters.template get<mio::oseirmetapopwang::CommutingStrengths<>>().get_cont_freq_mat()[0].get_baseline() =
+        mobility_data_commuter;
+
+    mio::ContactMatrixGroup& contact_matrix =
+        model.parameters.template get<mio::oseirmetapopwang::ContactPatterns<>>().get_cont_freq_mat();
+    contact_matrix[0].get_baseline() = get_contact_matrix();
+
+    for (size_t j = 0; j < number_age_groups; j++) {
+        model.parameters.template get<mio::oseirmetapopwang::TimeExposed<>>()[mio::AgeGroup(j)]  = TimeExposed[j];
+        model.parameters.template get<mio::oseirmetapopwang::TimeInfected<>>()[mio::AgeGroup(j)] = TimeInfected[j];
+        model.parameters.template get<mio::oseirmetapopwang::TransmissionProbabilityOnContact<>>()[mio::AgeGroup(j)] =
+            TransmissionProbabilityOnContact[j];
+    }
+
+    std::shared_ptr<mio::IntegratorCore<ScalarType>> integrator = std::make_shared<mio::EulerIntegratorCore<>>();
+
+    auto result = simulate(t0, tmax, dt, model, integrator);
+
+    auto basic_reproduction_number = model.get_reproduction_number(t0, result).value();
+    std::cout << "\"Model B\": " << basic_reproduction_number << "}" << std::endl;
+}
+
+int main()
+{
+    const ScalarType tmax = 1.;
+    size_t num_regions    = 1;
+
+    std::cout << "{ \"Regions\": " << num_regions << ", " << std::endl;
+
+    calculate_basic_reproduction_number(num_regions, tmax);
+    return 0;
+}

cpp/examples/examples_thesis/basic_reproduction_number_modelc.cpp

@@ -0,0 +1,98 @@
+#include "models/ode_metapop/model.h"
+
+#include "memilio/math/euler.h"
+#include "memilio/compartments/simulation.h"
+#include "memilio/utils/custom_index_array.h"
+
+Eigen::MatrixXd get_contact_matrix()
+{
+    Eigen::MatrixXd contact_matrix_eigen(6, 6);
+    contact_matrix_eigen << 3.9547, 1.1002, 2.9472, 2.05, 0.3733, 0.0445, 0.3327, 3.5892, 1.236, 1.9208, 0.2681, 0.0161,
+        0.246, 0.7124, 5.6518, 3.2939, 0.2043, 0.0109, 0.1742, 0.8897, 3.3124, 4.5406, 0.4262, 0.0214, 0.0458, 0.1939,
+        0.5782, 1.3825, 1.473, 0.0704, 0.1083, 0.1448, 0.4728, 0.9767, 0.6266, 0.1724;
+
+    return contact_matrix_eigen;
+}
+
+const ScalarType TimeExposed[]                      = {3.335, 3.335, 3.335, 3.335, 3.335, 3.335};
+const ScalarType TimeInfected[]                     = {8.0096875, 8.0096875, 8.2182, 8.1158, 8.033, 7.985};
+const ScalarType TransmissionProbabilityOnContact[] = {0.03, 0.06, 0.06, 0.06, 0.09, 0.175};
+
+void calculate_basic_reproduction_number(size_t number_regions, ScalarType tmax)
+{
+    mio::set_log_level(mio::LogLevel::off);
+    ScalarType t0                = 0.;
+    ScalarType dt                = 0.1;
+    ScalarType number_age_groups = 6;
+
+    mio::oseirmetapop::Model<ScalarType> model(number_regions, number_age_groups);
+    auto& population = model.populations;
+
+    for (size_t j = 0; j < number_age_groups; j++) {
+        for (size_t i = 0; i < number_regions; i++) {
+            population[{mio::oseirmetapop::Region(i), mio::AgeGroup(j),
+                        mio::oseirmetapop::InfectionState::Susceptible}] = 10000;
+        }
+    }
+
+    double fraction_commuter = 1. / (2 * number_regions);
+    Eigen::MatrixXd mobility_data_commuter =
+        Eigen::MatrixXd::Constant(number_regions, number_regions, fraction_commuter) -
+        fraction_commuter *
+            Eigen::MatrixXd::Identity(number_regions, number_regions); // Ensure that the diagonal is zero
+    for (size_t county_idx_i = 0; county_idx_i < number_regions; ++county_idx_i) {
+        mobility_data_commuter(county_idx_i, county_idx_i) = 1 - mobility_data_commuter.rowwise().sum()(county_idx_i);
+    }
+    model.parameters.template get<mio::oseirmetapop::CommutingStrengths<>>().get_cont_freq_mat()[0].get_baseline() =
+        mobility_data_commuter;
+
+    mio::ContactMatrixGroup& contact_matrix =
+        model.parameters.template get<mio::oseirmetapop::ContactPatterns<>>().get_cont_freq_mat();
+    contact_matrix[0].get_baseline() = get_contact_matrix();
+
+    for (size_t j = 0; j < number_age_groups; j++) {
+        model.parameters.template get<mio::oseirmetapop::TimeExposed<>>()[mio::AgeGroup(j)]  = TimeExposed[j];
+        model.parameters.template get<mio::oseirmetapop::TimeInfected<>>()[mio::AgeGroup(j)] = TimeInfected[j];
+        model.parameters.template get<mio::oseirmetapop::TransmissionProbabilityOnContact<>>()[mio::AgeGroup(j)] =
+            TransmissionProbabilityOnContact[j];
+    }
+
+    mio::ContactMatrixGroup& commuting_strengths =
+        model.parameters.template get<mio::oseirmetapop::CommutingStrengths<>>().get_cont_freq_mat();
+
+    auto& population_after_commuting = model.m_population_after_commuting;
+    for (size_t region_n = 0; region_n < number_regions; ++region_n) {
+        for (size_t age = 0; age < number_age_groups; ++age) {
+            double population_n = 0;
+            for (size_t state = 0; state < (size_t)mio::oseirmetapop::InfectionState::Count; state++) {
+                population_n += population[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age),
+                                            mio::oseirmetapop::InfectionState(state)}];
+            }
+            population_after_commuting[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age)}] += population_n;
+            for (size_t region_m = 0; region_m < number_regions; ++region_m) {
+                population_after_commuting[{mio::oseirmetapop::Region(region_n), mio::AgeGroup(age)}] -=
+                    commuting_strengths[0].get_baseline()(region_n, region_m) * population_n;
+                population_after_commuting[{mio::oseirmetapop::Region(region_m), mio::AgeGroup(age)}] +=
+                    commuting_strengths[0].get_baseline()(region_n, region_m) * population_n;
+            }
+        }
+    }
+
+    std::shared_ptr<mio::IntegratorCore<ScalarType>> integrator = std::make_shared<mio::EulerIntegratorCore<>>();
+
+    auto result = simulate(t0, tmax, dt, model, integrator);
+
+    auto basic_reproduction_number = model.get_reproduction_number(t0, result).value();
+    std::cout << "\"Metapopulation\": " << basic_reproduction_number << "}" << std::endl;
+}
+
+int main()
+{
+    const ScalarType tmax = 1.;
+    size_t num_regions    = 1;
+
+    std::cout << "{ \"Regions\": " << num_regions << ", " << std::endl;
+
+    calculate_basic_reproduction_number(num_regions, tmax);
+    return 0;
+}