687 create minimal hybrid graph model

cpp/CMakeLists.txt

@@ -116,6 +116,7 @@ add_subdirectory(memilio)
 
 if(MEMILIO_BUILD_MODELS)
     add_subdirectory(models/abm)
+    add_subdirectory(models/graph_abm)
     add_subdirectory(models/ode_secir)
     add_subdirectory(models/ode_secirvvs)
     add_subdirectory(models/ide_secir)

cpp/examples/CMakeLists.txt

@@ -35,6 +35,10 @@ target_link_libraries(abm_minimal_example PRIVATE memilio abm)
 add_executable(abm_history_example abm_history_object.cpp)
 target_link_libraries(abm_history_example PRIVATE memilio abm)
 
+add_executable(graph_abm_example graph_abm.cpp)
+target_link_libraries(graph_abm_example PRIVATE memilio graph_abm)
+target_link_libraries(graph_abm_example PRIVATE memilio abm)
+
 add_executable(twitter_migration_example twitter_migration.cpp)
 target_link_libraries(twitter_migration_example PRIVATE memilio ode_secir)
 

cpp/examples/abm_minimal.cpp

@@ -17,8 +17,8 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
-#include "abm/abm.h"
-#include "abm/household.h"
+#include "models/abm/abm.h"
+#include "models/abm/household.h"
 #include <cstdio>
 #include "abm/world.h"
 #include "memilio/io/io.h"

cpp/examples/graph_abm.cpp

@@ -0,0 +1,230 @@
+/*
+* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC)
+*
+* Authors: Khoa Nguyen
+*
+* 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 "models/abm/abm.h"
+#include "graph_abm/graph_world.h"
+#include "models/abm/household.h"
+#include "memilio/mobility/graph.h"
+#include "memilio/mobility/metapopulation_mobility_instant.h"
+#include <cstdio>
+
+int main()
+{
+    // Set global infection parameters (similar to infection parameters in SECIR model) and initialize the world
+    mio::abm::GlobalInfectionParameters infection_params;
+
+    // Set same infection parameter for all age groups. For example, the incubation period is 4 days.
+    infection_params.get<mio::abm::IncubationPeriod>() = 4.;
+
+    // Create the worlds with infection parameters.
+    auto world1 = mio::graph_abm::GraphWorld(infection_params, 0);
+    auto world2 = mio::graph_abm::GraphWorld(infection_params, 1);
+
+    //add households for world 1
+    auto child = mio::abm::HouseholdMember(); // A child is 50/50% 0-4 or 5-14.
+    child.set_age_weight(mio::abm::AgeGroup::Age0to4, 1);
+    child.set_age_weight(mio::abm::AgeGroup::Age5to14, 1);
+
+    auto parent = mio::abm::HouseholdMember(); // A parent is 50/50% 15-34 or 35-59.
+    parent.set_age_weight(mio::abm::AgeGroup::Age15to34, 1);
+    parent.set_age_weight(mio::abm::AgeGroup::Age35to59, 1);
+
+    auto two_person_household_group1 = mio::abm::HouseholdGroup();
+    auto two_person_household_group2 = mio::abm::HouseholdGroup();
+    auto two_person_family           = mio::abm::Household();
+    two_person_family.add_members(child, 1);
+    two_person_family.add_members(parent, 1);
+    auto two_person_adults = mio::abm::Household();
+    two_person_adults.add_members(parent, 2);
+    two_person_household_group1.add_households(two_person_family, 2);
+    two_person_household_group1.add_households(two_person_adults, 1);
+    two_person_household_group2.add_households(two_person_family, 1);
+    two_person_household_group2.add_households(two_person_adults, 3);
+
+    auto three_person_household_group1 = mio::abm::HouseholdGroup();
+    auto three_person_household_group2 = mio::abm::HouseholdGroup();
+    auto three_person_family           = mio::abm::Household();
+    three_person_family.add_members(child, 1);
+    three_person_family.add_members(parent, 2);
+    three_person_household_group1.add_households(three_person_family, 2);
+    three_person_household_group2.add_households(three_person_family, 1);
+
+    //add households to world 1
+    add_household_group_to_world(world1, two_person_household_group1);
+    add_household_group_to_world(world1, three_person_household_group1);
+    //add households to world 2
+    add_household_group_to_world(world2, two_person_household_group2);
+    add_household_group_to_world(world2, three_person_household_group2);
+
+    //add internal locations for both worlds
+    //world1
+    std::vector<mio::abm::LocationId> events_w1;
+    std::vector<mio::abm::LocationId> works_w1;
+    //add two social events
+    auto event1_w1 = world1.add_location(mio::abm::LocationType::SocialEvent);
+    world1.get_individualized_location(event1_w1).get_infection_parameters().set<mio::abm::MaximumContacts>(5);
+    auto event2_w1 = world1.add_location(mio::abm::LocationType::SocialEvent);
+    world1.get_individualized_location(event2_w1).get_infection_parameters().set<mio::abm::MaximumContacts>(10);
+    events_w1.push_back(event1_w1);
+    events_w1.push_back(event2_w1);
+    // Add hospital and ICU with 5 maximum contacs.
+    auto hospital = world1.add_location(mio::abm::LocationType::Hospital);
+    world1.get_individualized_location(hospital).get_infection_parameters().set<mio::abm::MaximumContacts>(5);
+    auto icu = world1.add_location(mio::abm::LocationType::ICU);
+    world1.get_individualized_location(icu).get_infection_parameters().set<mio::abm::MaximumContacts>(5);
+    // Add one supermarket, maximum constacts are assumed to be 20.
+    auto shop_w1 = world1.add_location(mio::abm::LocationType::BasicsShop);
+    world1.get_individualized_location(shop_w1).get_infection_parameters().set<mio::abm::MaximumContacts>(20);
+    // At every school, the maximum contacts are 20.
+    auto school = world1.add_location(mio::abm::LocationType::School);
+    world1.get_individualized_location(school).get_infection_parameters().set<mio::abm::MaximumContacts>(20);
+    // At every workplace, maximum contacts are 10.
+    auto work1_w1 = world1.add_location(mio::abm::LocationType::Work);
+    world1.get_individualized_location(work1_w1).get_infection_parameters().set<mio::abm::MaximumContacts>(10);
+    auto work2_w1 = world1.add_location(mio::abm::LocationType::Work);
+    world1.get_individualized_location(work2_w1).get_infection_parameters().set<mio::abm::MaximumContacts>(10);
+    works_w1.push_back(work1_w1);
+    works_w1.push_back(work2_w1);
+
+    //world2
+    std::vector<mio::abm::LocationId> events_w2;
+    std::vector<mio::abm::LocationId> works_w2;
+    //add one social event
+    auto event_w2 = world2.add_location(mio::abm::LocationType::SocialEvent);
+    world2.get_individualized_location(event_w2).get_infection_parameters().set<mio::abm::MaximumContacts>(5);
+    events_w2.push_back(event_w2);
+    // Add one supermarket, maximum constacts are assumed to be 20.
+    auto shop_w2 = world2.add_location(mio::abm::LocationType::BasicsShop);
+    world2.get_individualized_location(shop_w2).get_infection_parameters().set<mio::abm::MaximumContacts>(20);
+    // At every workplace, maximum contacts are 10.
+    auto work_w2 = world2.add_location(mio::abm::LocationType::Work);
+    world2.get_individualized_location(work_w2).get_infection_parameters().set<mio::abm::MaximumContacts>(10);
+    works_w2.push_back(work_w2);
+
+    //add external locations
+    //world 1
+    works_w1.push_back(work_w2);
+
+    //world 2
+    //TODO: MaximumContacts for external locations setzen
+    events_w2.push_back(event2_w1);
+    works_w2.push_back(work2_w1);
+
+    auto start_date = mio::abm::TimePoint(0);
+
+    //assign infection states
+    //world 1
+    auto persons_w1 = world1.get_persons();
+    for (auto& person : persons_w1) {
+        mio::abm::InfectionState infection_state =
+            (mio::abm::InfectionState)(rand() % ((uint32_t)mio::abm::InfectionState::Count - 1));
+        if (infection_state != mio::abm::InfectionState::Susceptible)
+            person.add_new_infection(mio::abm::Infection(mio::abm::VirusVariant::Wildtype, person.get_age(),
+                                                         world1.get_global_infection_parameters(), start_date,
+                                                         infection_state));
+    }
+    //world 2
+    auto persons_w2 = world2.get_persons();
+    for (auto& person : persons_w2) {
+        mio::abm::InfectionState infection_state =
+            (mio::abm::InfectionState)(rand() % ((uint32_t)mio::abm::InfectionState::Count - 1));
+        if (infection_state != mio::abm::InfectionState::Susceptible)
+            person.add_new_infection(mio::abm::Infection(mio::abm::VirusVariant::Wildtype, person.get_age(),
+                                                         world2.get_global_infection_parameters(), start_date,
+                                                         infection_state));
+    }
+
+    //assign locations
+    //word 1
+    for (auto& person : persons_w1) {
+        //choose one of the events with same probability and assign to person
+        auto event =
+            events_w1[mio::DiscreteDistribution<size_t>::get_instance()(std::vector<double>(events_w1.size(), 1.0))];
+        person.set_assigned_location(event);
+        //assign shop
+        person.set_assigned_location(shop_w1);
+        //assign hospital and ICU
+        person.set_assigned_location(hospital);
+        person.set_assigned_location(icu);
+        //assign work and school dependent on person's age
+        if (person.get_age() == mio::abm::AgeGroup::Age5to14) {
+            person.set_assigned_location(school);
+        }
+        if (person.get_age() == mio::abm::AgeGroup::Age15to34 || person.get_age() == mio::abm::AgeGroup::Age35to59) {
+            auto work =
+                works_w1[mio::DiscreteDistribution<size_t>::get_instance()(std::vector<double>(works_w1.size(), 1.0))];
+            person.set_assigned_location(work);
+        }
+    }
+    //world 2
+    for (auto& person : persons_w2) {
+        //choose one of the events with same probability and assign to person
+        auto event =
+            events_w2[mio::DiscreteDistribution<size_t>::get_instance()(std::vector<double>(events_w2.size(), 1.0))];
+        person.set_assigned_location(event);
+        //assign shop
+        person.set_assigned_location(shop_w2);
+        //assign hospital and ICU
+        person.set_assigned_location(hospital);
+        person.set_assigned_location(icu);
+        // Assign work and school dependent on Person's AgeGroup.
+        if (person.get_age() == mio::abm::AgeGroup::Age5to14) {
+            person.set_assigned_location(school);
+        }
+        if (person.get_age() == mio::abm::AgeGroup::Age15to34 || person.get_age() == mio::abm::AgeGroup::Age35to59) {
+            auto work =
+                works_w2[mio::DiscreteDistribution<size_t>::get_instance()(std::vector<double>(works_w2.size(), 1.0))];
+            person.set_assigned_location(work);
+        }
+    }
+
+    auto t0   = mio::abm::TimePoint(0);
+    auto dt   = mio::abm::hours(12);
+    auto tmax = mio::abm::TimePoint(0) + mio::abm::days(10);
+
+    mio::Graph<mio::SimulationNode<mio::graph_abm::GraphSimulation>, mio::MobilityEdgeAgents> g;
+    g.add_node(0, t0, std::move(world1));
+    g.add_node(1, t0, std::move(world2));
+    g.add_edge(0, 1);
+    g.add_edge(1, 0);
+
+    auto sim = mio::make_migration_sim(t0, dt, std::move(g));
+
+    sim.advance(tmax, true);
+    for (auto& n : sim.get_graph().nodes()) {
+        std::cout << "node " << n.id << "\n";
+        std::cout << "\n";
+
+        std::cout << "# t S E C I I_s I_c R_C R_I D\n";
+        for (auto i = 0; i < n.property.get_result().get_num_time_points(); ++i) {
+            std::cout << n.property.get_result().get_time(i) << " ";
+            auto v = n.property.get_result().get_value(i);
+            for (auto j = 0; j < v.size(); ++j) {
+                std::cout << v[j] << " ";
+                if (j < v.size() - 1) {
+                    std::cout << " ";
+                }
+            }
+            if (i < n.property.get_result().get_num_time_points() - 1) {
+                std::cout << "\n";
+            }
+        }
+        std::cout << "\n";
+    }
+}

cpp/examples/ode_secir_parameter_study.cpp

@@ -30,9 +30,9 @@
  * @param t0 starting point of simulation
  * @param tmax end point of simulation
  */
-mio::IOResult<void>
-write_single_run_result(const int run,
-                        const mio::Graph<mio::SimulationNode<mio::osecir::Simulation<>>, mio::MigrationEdge>& graph)
+mio::IOResult<void> write_single_run_result(
+    const int run,
+    const mio::Graph<mio::SimulationNode<mio::osecir::Simulation<>>, mio::MigrationEdge>& graph)
 {
     std::string abs_path;
     BOOST_OUTCOME_TRY(created, mio::create_directory("results", abs_path));
@@ -150,7 +150,7 @@ int main()
         auto write_result_status = write_single_run_result(run, graph);
         if (!write_result_status) {
             std::cout << "Error writing result: " << write_result_status.error().formatted_message();
-        }       
+        }
         return 0; //Result handler must return something, but only meaningful when using MPI.
     };
     parameter_study.run(sample_graph, handle_result);

cpp/memilio/CMakeLists.txt

@@ -57,10 +57,16 @@ add_library(memilio
     mobility/metapopulation_mobility_instant.cpp
     mobility/metapopulation_mobility_stochastic.h
     mobility/metapopulation_mobility_stochastic.cpp
+    mobility/metapopulation_mobility_hybrid.h
+    mobility/metapopulation_mobility_hybrid.cpp
     mobility/graph_simulation.h
     mobility/graph_simulation.cpp
     mobility/graph.h
     mobility/graph.cpp
+    mobility/hybrid_graph.h
+    mobility/hybrid_graph.cpp
+    mobility/hybrid_graph_simulation.h
+    mobility/hybrid_graph_simulation.cpp
     utils/visitor.h
     utils/uncertain_value.h
     utils/uncertain_value.cpp

cpp/memilio/compartments/parameter_studies.h

@@ -141,7 +141,8 @@ class ParameterStudy
 #endif
 
         auto run_distribution = distribute_runs(m_num_runs, num_procs);
-        auto start_run_idx = std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0));
+        auto start_run_idx =
+            std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0));
         auto end_run_idx = start_run_idx + run_distribution[size_t(rank)];
 
         std::vector<std::result_of_t<HandleSimulationResultFunction(SimulationGraph, size_t)>> ensemble_result;
@@ -154,7 +155,7 @@ class ParameterStudy
             //- the RNG has been freshly seeded/initialized before this call
             //- the seeds are identical on all MPI processes
             //- the block size of the RNG is sufficiently big to cover one run
-            //  (when in doubt, use a larger block size; fast-forwarding the RNG is cheap and the period length 
+            //  (when in doubt, use a larger block size; fast-forwarding the RNG is cheap and the period length
             //   of the mt19937 RNG is huge)
             mio::thread_local_rng().forward_to_block(run_idx);
 
@@ -296,7 +297,7 @@ class ParameterStudy
 private:
     //sample parameters and create simulation
     template <class SampleGraphFunction>
-    mio::GraphSimulation<SimulationGraph> create_sampled_simulation(SampleGraphFunction sample_graph)
+    mio::GraphSimulation<SimulationGraph, double, double> create_sampled_simulation(SampleGraphFunction sample_graph)
     {
         SimulationGraph sim_graph;
 
@@ -316,7 +317,7 @@ class ParameterStudy
         //evenly distribute runs
         //lower processes do one more run if runs are not evenly distributable
         auto num_runs_local = num_runs / num_procs; //integer division!
-        auto remainder = num_runs % num_procs;
+        auto remainder      = num_runs % num_procs;
 
         std::vector<size_t> run_distribution(num_procs);
         std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1);

cpp/memilio/io/mobility_io.h

@@ -133,7 +133,8 @@ IOResult<void> write_graph(const Graph<Model, MigrationParameters>& graph, const
  * @param read_edges boolean value that decides whether the edges of the graph should also be read in.
  */
 template <class Model>
-IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& directory, int ioflags = IOF_None, bool read_edges = true)
+IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& directory, int ioflags = IOF_None,
+                                                       bool read_edges = true)
 {
     std::string abs_path;
     if (!file_exists(directory, abs_path)) {
@@ -160,7 +161,7 @@ IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& direct
     }
 
     //read edges; nodes must already be available for that)
-    if(read_edges){
+    if (read_edges) {
         for (auto inode = size_t(0); inode < graph.nodes().size(); ++inode) {
             //list of edges
             auto edge_filename = path_join(abs_path, "GraphEdges_node" + std::to_string(inode) + ".json");
@@ -177,7 +178,7 @@ IOResult<Graph<Model, MigrationParameters>> read_graph(const std::string& direct
                 if (end_node_idx >= graph.nodes().size()) {
                     log_error("EndNodeIndex not in range of number of graph nodes.");
                     return failure(StatusCode::OutOfRange,
-                                edge_filename + ", EndNodeIndex not in range of number of graph nodes.");
+                                   edge_filename + ", EndNodeIndex not in range of number of graph nodes.");
                 }
                 BOOST_OUTCOME_TRY(parameters, deserialize_json(e["Parameters"], Tag<MigrationParameters>{}, ioflags));
                 graph.add_edge(start_node_idx, end_node_idx, parameters);