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

cpp/CMakeLists.txt

@@ -151,6 +151,10 @@ if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/ode_seir)
     add_subdirectory(models/ode_seair)
     add_subdirectory(models/ode_sir)
+    add_subdirectory(models/ode_sir_mobility)
+    # add_subdirectory(models/ode_seir_mobility_massaction)
+    add_subdirectory(models/ode_seir_mobility_improved)
+    add_subdirectory(models/ode_seir_mobility)
     add_subdirectory(models/sde_sir)
     add_subdirectory(models/sde_sirs)
     add_subdirectory(models/sde_seirvv)

cpp/examples/CMakeLists.txt

@@ -30,6 +30,24 @@ 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_sir_mobility_example ode_sir_mobility.cpp)
+target_link_libraries(ode_sir_mobility_example PRIVATE memilio ode_sir_mobility)
+target_compile_options(ode_sir_mobility_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+add_executable(ode_seir_mobility_example ode_seir_mobility.cpp)
+target_link_libraries(ode_seir_mobility_example PRIVATE memilio ode_seir_mobility)
+target_compile_options(ode_seir_mobility_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+if (MEMILIO_ENABLE_OPENMP)
+    add_executable(ode_seir_mobility_timing ode_seir_mobility_timing.cpp)
+    target_link_libraries(ode_seir_mobility_timing PRIVATE memilio ode_seir_mobility_improved)
+    target_compile_options(ode_seir_mobility_timing PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+endif()
+
+add_executable(ode_seir_mobility_example_improved ode_seir_mobility_improved.cpp)
+target_link_libraries(ode_seir_mobility_example_improved PRIVATE memilio ode_seir_mobility_improved)
+target_compile_options(ode_seir_mobility_example_improved 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})
@@ -84,6 +102,16 @@ add_executable(graph_example graph.cpp)
 target_link_libraries(graph_example PRIVATE memilio ode_seir)
 target_compile_options(graph_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(graph_example_extended graph_extended.cpp)
+target_link_libraries(graph_example_extended PRIVATE memilio ode_seir)
+target_compile_options(graph_example_extended PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
+if (MEMILIO_ENABLE_OPENMP)
+    add_executable(graph_timing graph_timing.cpp)
+    target_link_libraries(graph_timing PRIVATE memilio ode_seir)
+    target_compile_options(graph_timing PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+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/graph.cpp

@@ -22,21 +22,24 @@
 #include "ode_seir/parameters.h"
 #include "memilio/mobility/metapopulation_mobility_instant.h"
 #include "memilio/compartments/simulation.h"
+#include "memilio/io/result_io.h"
 
 int main()
 {
     const auto t0   = 0.;
-    const auto tmax = 10.;
+    const auto tmax = 15.;
     const auto dt   = 0.5; //time step of mobility, daily mobility every second step
 
     mio::oseir::Model<> model(1);
 
     // set population
     model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 10000;
 
+    model.parameters.set<mio::oseir::TransmissionProbabilityOnContact<>>(1.);
+
     // set transition times
-    model.parameters.set<mio::oseir::TimeExposed<>>(1);
-    model.parameters.set<mio::oseir::TimeInfected<>>(1);
+    model.parameters.set<mio::oseir::TimeExposed<>>(3.);
+    model.parameters.set<mio::oseir::TimeInfected<>>(5.);
 
     // set contact matrix
     mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::oseir::ContactPatterns<>>().get_cont_freq_mat();
@@ -47,9 +50,9 @@ int main()
     auto model_group2 = model;
 
     //some contact restrictions in group 1
-    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));
+    // 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::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 9990;
@@ -58,12 +61,29 @@ int main()
     mio::Graph<mio::SimulationNode<mio::Simulation<ScalarType, mio::oseir::Model<>>>, mio::MobilityEdge<>> g;
     g.add_node(1001, model_group1, t0);
     g.add_node(1002, model_group2, t0);
-    g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.01));
+    for (auto& node : g.nodes()) {
+        node.property.get_simulation().set_integrator(std::make_shared<mio::EulerIntegratorCore<ScalarType>>());
+    }
+    g.add_edge(0, 1, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.05));
     g.add_edge(1, 0, Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count, 0.01));
 
     auto sim = mio::make_mobility_sim(t0, dt, std::move(g));
 
     sim.advance(tmax);
 
+    auto result_graph = std::move(sim).get_graph();
+    auto result       = mio::interpolate_simulation_result(result_graph);
+
+    std::vector<int> county_ids(result_graph.nodes().size());
+    std::transform(result_graph.nodes().begin(), result_graph.nodes().end(), county_ids.begin(), [](auto& n) {
+        return n.id;
+    });
+
+    auto save_result_status = save_result(result, county_ids, 1, "graph_result.h5");
+
+    for (auto&& node : result_graph.nodes()) {
+        node.property.get_result().print_table();
+    }
+
     return 0;
 }

cpp/examples/graph_extended.cpp

@@ -0,0 +1,271 @@
+
+#include "ode_seir/model.h"
+#include "ode_seir/infection_state.h"
+#include "ode_seir/parameters.h"
+#include "memilio/mobility/metapopulation_mobility_instant.h"
+#include "memilio/compartments/simulation.h"
+#include "memilio/io/result_io.h"
+#include "memilio/io/epi_data.h"
+
+#include <chrono>
+
+/**
+ * indices of contact matrix corresponding to locations where contacts occur.
+ */
+enum class ContactLocation
+{
+    Home = 0,
+    School,
+    Work,
+    Other,
+    Count,
+};
+
+static const std::map<ContactLocation, std::string> contact_locations = {{ContactLocation::Home, "home"},
+                                                                         {ContactLocation::School, "school_pf_eig"},
+                                                                         {ContactLocation::Work, "work"},
+                                                                         {ContactLocation::Other, "other"}};
+
+/**
+ * Set contact matrices.
+ * Reads contact matrices from files in the data directory.
+ * @param data_dir data directory.
+ * @param params Object that the contact matrices will be added to.
+ * @returns any io errors that happen during reading of the files.
+ */
+mio::IOResult<void> set_contact_matrices(const fs::path& data_dir, mio::oseir::Parameters<double>& params,
+                                         bool synthetic_population)
+{
+    if (!synthetic_population) {
+        //TODO: io error handling
+        auto contact_matrices = mio::ContactMatrixGroup(contact_locations.size(), size_t(params.get_num_groups()));
+        for (auto&& contact_location : contact_locations) {
+            BOOST_OUTCOME_TRY(auto&& baseline,
+                              mio::read_mobility_plain(
+                                  (data_dir / "contacts" / ("baseline_" + contact_location.second + ".txt")).string()));
+            BOOST_OUTCOME_TRY(auto&& minimum,
+                              mio::read_mobility_plain(
+                                  (data_dir / "contacts" / ("minimum_" + contact_location.second + ".txt")).string()));
+            contact_matrices[size_t(contact_location.first)].get_baseline() = baseline;
+            contact_matrices[size_t(contact_location.first)].get_minimum()  = minimum;
+        }
+        params.get<mio::oseir::ContactPatterns<double>>() = mio::UncertainContactMatrix<double>(contact_matrices);
+    }
+    else {
+        mio::ContactMatrixGroup& contact_matrix = params.get<mio::oseir::ContactPatterns<>>().get_cont_freq_mat();
+        contact_matrix[0].get_baseline().setConstant(7.95 / (size_t)params.get_num_groups());
+    }
+
+    printf("Setting contact matrices successful.\n");
+    return mio::success();
+}
+
+/**
+ * Set epidemiological parameters of Sars-CoV-2 for a immune-naive
+ * population and wild type variant.
+ * @param params Object that the parameters will be added to.
+ * @returns Currently generates no errors.
+ */
+mio::IOResult<void> set_covid_parameters(mio::oseir::Parameters<double>& params, bool synthetic_population)
+{
+    params.template set<mio::oseir::TimeExposed<>>(3.335);
+    if (!synthetic_population) {
+        params.get<mio::oseir::TimeInfected<>>()[mio::AgeGroup(0)] = 8.0096875;
+        params.get<mio::oseir::TimeInfected<>>()[mio::AgeGroup(1)] = 8.0096875;
+        params.get<mio::oseir::TimeInfected<>>()[mio::AgeGroup(2)] = 8.2182;
+        params.get<mio::oseir::TimeInfected<>>()[mio::AgeGroup(3)] = 8.1158;
+        params.get<mio::oseir::TimeInfected<>>()[mio::AgeGroup(4)] = 8.033;
+        params.get<mio::oseir::TimeInfected<>>()[mio::AgeGroup(5)] = 7.985;
+
+        params.get<mio::oseir::TransmissionProbabilityOnContact<>>()[mio::AgeGroup(0)] = 0.03;
+        params.get<mio::oseir::TransmissionProbabilityOnContact<>>()[mio::AgeGroup(1)] = 0.06;
+        params.get<mio::oseir::TransmissionProbabilityOnContact<>>()[mio::AgeGroup(2)] = 0.06;
+        params.get<mio::oseir::TransmissionProbabilityOnContact<>>()[mio::AgeGroup(3)] = 0.06;
+        params.get<mio::oseir::TransmissionProbabilityOnContact<>>()[mio::AgeGroup(4)] = 0.09;
+        params.get<mio::oseir::TransmissionProbabilityOnContact<>>()[mio::AgeGroup(5)] = 0.175;
+    }
+    else {
+        params.template set<mio::oseir::TimeInfected<>>(8.097612257);
+
+        params.template set<mio::oseir::TransmissionProbabilityOnContact<>>(0.07333);
+    }
+
+    printf("Setting epidemiological parameters successful.\n");
+    return mio::success();
+}
+
+mio::IOResult<std::vector<mio::oseir::Model<double>>>
+set_population_data(const fs::path& data_dir, mio::oseir::Parameters<double>& params, std::vector<int> node_ids)
+{
+    size_t number_regions = node_ids.size();
+
+    std::vector<mio::oseir::Model<double>> nodes(number_regions,
+                                                 mio::oseir::Model(int(size_t(params.get_num_groups()))));
+
+    for (auto& node : nodes) {
+        node.parameters = params;
+    }
+
+    BOOST_OUTCOME_TRY(const auto&& population_data,
+                      mio::read_population_data(
+                          (data_dir / "pydata" / "Germany" / "county_current_population_nrw.json").string(), true));
+
+    std::vector<std::vector<double>> vnum_population(node_ids.size(),
+                                                     std::vector<double>((size_t)params.get_num_groups(), 0.0));
+
+    for (auto&& entry : population_data) {
+        auto it = std::find_if(node_ids.begin(), node_ids.end(), [&entry](auto r) {
+            return r == 0 ||
+                   (entry.county_id && mio::regions::StateId(r) == mio::regions::get_state_id(int(*entry.county_id))) ||
+                   (entry.county_id && mio::regions::CountyId(r) == *entry.county_id) ||
+                   (entry.district_id && mio::regions::DistrictId(r) == *entry.district_id);
+        });
+        if (it != node_ids.end()) {
+            auto region_idx      = size_t(it - node_ids.begin());
+            auto& num_population = vnum_population[region_idx];
+            for (size_t age = 0; age < num_population.size(); age++) {
+                num_population[age] += entry.population[mio::AgeGroup(age)];
+            }
+        }
+    }
+
+    for (size_t region = 0; region < node_ids.size(); region++) {
+        auto num_groups = nodes[region].parameters.get_num_groups();
+        for (auto i = mio::AgeGroup(0); i < num_groups; i++) {
+            nodes[region].populations.template set_difference_from_group_total<mio::AgeGroup>(
+                {i, mio::oseir::InfectionState::Susceptible}, vnum_population[region][size_t(i)]);
+        }
+    }
+    nodes[27].populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] -= 100;
+    nodes[27].populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] += 100;
+
+    return mio::success(nodes);
+}
+
+mio::IOResult<std::vector<mio::oseir::Model<double>>>
+set_synthetic_population_data(mio::oseir::Parameters<double>& params)
+{
+    size_t number_regions = 53;
+
+    std::vector<mio::oseir::Model<double>> nodes(number_regions,
+                                                 mio::oseir::Model(int(size_t(params.get_num_groups()))));
+
+    mio::Populations<double, mio::AgeGroup, mio::oseir::InfectionState> population(
+        {params.get_num_groups(), mio::oseir::InfectionState::Count});
+
+    for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) {
+        population[{i, mio::oseir::InfectionState::Susceptible}] = 1000000;
+    }
+    for (auto& node : nodes) {
+        node.parameters  = params;
+        node.populations = population;
+    }
+    for (auto i = mio::AgeGroup(0); i < params.get_num_groups(); i++) {
+        nodes[0].populations[{i, mio::oseir::InfectionState::Exposed}]     = 100;
+        nodes[0].populations[{i, mio::oseir::InfectionState::Susceptible}] = 999900;
+    }
+
+    return mio::success(nodes);
+}
+
+mio::IOResult<void> run(const fs::path& data_dir, double t0, double tmax, double dt)
+{
+    mio::set_log_level(mio::LogLevel::off);
+    // global parameters
+    bool synthetic_population = false;
+    const int num_age_groups  = 6;
+    if (num_age_groups != 6) {
+        synthetic_population = true;
+    }
+    mio::oseir::Parameters<double> params(num_age_groups);
+
+    BOOST_OUTCOME_TRY(set_covid_parameters(params, synthetic_population));
+
+    // set contact matrix
+    BOOST_OUTCOME_TRY(set_contact_matrices(data_dir, params, synthetic_population));
+
+    // graph of counties with populations and local parameters
+    // and mobility between counties
+    mio::Graph<mio::SimulationNode<mio::Simulation<ScalarType, mio::oseir::Model<>>>, mio::MobilityEdge<>> params_graph;
+
+    BOOST_OUTCOME_TRY(
+        auto&& node_ids,
+        mio::get_node_ids((data_dir / "pydata" / "Germany" / "county_current_population_nrw.json").string(), true,
+                          true));
+
+    if (synthetic_population) {
+        BOOST_OUTCOME_TRY(auto&& nodes, set_synthetic_population_data(params));
+        for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) {
+            params_graph.add_node(node_ids[node_idx], nodes[node_idx]);
+        }
+        printf("Setting synthetic population successful.\n");
+    }
+    else {
+        BOOST_OUTCOME_TRY(auto&& nodes, set_population_data(data_dir, params, node_ids));
+        for (size_t node_idx = 0; node_idx < nodes.size(); ++node_idx) {
+            params_graph.add_node(node_ids[node_idx], nodes[node_idx]);
+        }
+        printf("Setting population from data successful.\n");
+    }
+
+    BOOST_OUTCOME_TRY(auto&& mobility_data_commuter,
+                      mio::read_mobility_plain((data_dir / "mobility" / "commuter_mobility_nrw.txt").string()));
+    if (mobility_data_commuter.rows() != Eigen::Index(params_graph.nodes().size()) ||
+        mobility_data_commuter.cols() != Eigen::Index(params_graph.nodes().size())) {
+        return mio::failure(mio::StatusCode::InvalidValue,
+                            "Mobility matrices do not have the correct size. You may need to run "
+                            "transformMobilitydata.py from pycode memilio epidata package.");
+    }
+
+    for (size_t county_idx_i = 0; county_idx_i < params_graph.nodes().size(); ++county_idx_i) {
+        for (size_t county_idx_j = 0; county_idx_j < params_graph.nodes().size(); ++county_idx_j) {
+            auto& populations = params_graph.nodes()[county_idx_i].property.get_simulation().get_model().populations;
+
+            auto commuter_coeff_ij = mobility_data_commuter(county_idx_i, county_idx_j) / populations.get_total();
+            params_graph.add_edge(
+                county_idx_i, county_idx_j,
+                Eigen::VectorXd::Constant((size_t)mio::oseir::InfectionState::Count * size_t(params.get_num_groups()),
+                                          commuter_coeff_ij));
+        }
+    }
+
+    for (auto& node : params_graph.nodes()) {
+        node.property.get_simulation().set_integrator(std::make_shared<mio::EulerIntegratorCore<ScalarType>>());
+    }
+
+    auto sim = mio::make_mobility_sim(t0, dt, std::move(params_graph));
+
+    printf("Start Simulation\n");
+    auto t1 = std::chrono::high_resolution_clock::now();
+    sim.advance(tmax);
+    auto t2 = std::chrono::high_resolution_clock::now();
+
+    std::chrono::duration<double, std::milli> ms_double = t2 - t1;
+
+    printf("Runtime: %f\n", ms_double.count());
+
+    auto result_graph = std::move(sim).get_graph();
+    auto result       = mio::interpolate_simulation_result(result_graph);
+
+    std::vector<int> county_ids(result_graph.nodes().size());
+    std::transform(result_graph.nodes().begin(), result_graph.nodes().end(), county_ids.begin(), [](auto& n) {
+        return n.id;
+    });
+
+    auto save_result_status = save_result(result, county_ids, num_age_groups, "graph_result_nrw.h5");
+
+    return mio::success();
+}
+
+int main()
+{
+    const auto t0   = 0.;
+    const auto tmax = 50.;
+    const auto dt   = 0.5; //time step of mobility, daily mobility every second step
+
+    const std::string& data_dir = "";
+
+    auto result = run(data_dir, t0, tmax, dt);
+
+    return 0;
+}