1410 Accelerate add_edge for large graphs

cpp/memilio/mobility/graph.h

@@ -20,15 +20,17 @@
 #ifndef GRAPH_H
 #define GRAPH_H
 
-#include <functional>
 #include "memilio/utils/stl_util.h"
 #include "memilio/epidemiology/age_group.h"
 #include "memilio/utils/date.h"
 #include "memilio/utils/uncertain_value.h"
 #include "memilio/utils/parameter_distributions.h"
 #include "memilio/epidemiology/damping.h"
 #include "memilio/geography/regions.h"
+#include <algorithm>
+#include <functional>
 #include <iostream>
+#include <ranges>
 
 #include "boost/filesystem.hpp"
 
@@ -152,30 +154,34 @@ class Graph
     using NodeProperty = NodePropertyT;
     using EdgeProperty = EdgePropertyT;
 
+    Graph() = default;
+    Graph(std::vector<Node<NodePropertyT>>&& nodes, std::vector<Edge<EdgePropertyT>>&& edges)
+        : m_nodes(std::move(nodes))
+        , m_edges(std::move(edges))
+    {
+    }
+
     /**
      * @brief add a node to the graph. property of the node is constructed from arguments.
      */
     template <class... Args>
-    Node<NodePropertyT>& add_node(int id, Args&&... args)
+    void add_node(int id, Args&&... args)
     {
         m_nodes.emplace_back(id, std::forward<Args>(args)...);
-        return m_nodes.back();
     }
 
     /**
      * @brief add an edge to the graph. property of the edge is constructed from arguments.
      */
     template <class... Args>
-    Edge<EdgePropertyT>& add_edge(size_t start_node_idx, size_t end_node_idx, Args&&... args)
+    void add_edge(size_t start_node_idx, size_t end_node_idx, Args&&... args)
     {
         assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx);
-        return *insert_sorted_replace(m_edges,
-                                      Edge<EdgePropertyT>(start_node_idx, end_node_idx, std::forward<Args>(args)...),
-                                      [](auto&& e1, auto&& e2) {
-                                          return e1.start_node_idx == e2.start_node_idx
-                                                     ? e1.end_node_idx < e2.end_node_idx
-                                                     : e1.start_node_idx < e2.start_node_idx;
-                                      });
+        insert_sorted_replace(m_edges, Edge<EdgePropertyT>(start_node_idx, end_node_idx, std::forward<Args>(args)...),
+                              [](auto&& e1, auto&& e2) {
+                                  return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx
+                                                                                : e1.start_node_idx < e2.start_node_idx;
+                              });
     }
 
     /**
@@ -436,6 +442,112 @@ void print_graph(std::ostream& os, const Graph& g)
     }
 }
 
+/**
+ * @brief A builder class for constructing graphs.
+ *
+ * This class provides a interface for adding nodes and edges to a graph. It allows for efficient construction of large 
+ * graphs by reserving space for nodes and edges in advance. The build method finalizes the graph by sorting edges and 
+ * optionally removing duplicates.
+ * The advantage over the :ref add_edge function of the Graph class is that edges are only sorted once during the build 
+ * process, improving performance when adding many edges.
+ * 
+ * @tparam NodePropertyT Type of the node property.
+ * @tparam EdgePropertyT Type of the edge property.
+ */
+template <class NodePropertyT, class EdgePropertyT>
+class GraphBuilder
+{
+public:
+    using NodeProperty = NodePropertyT;
+    using EdgeProperty = EdgePropertyT;
+
+    GraphBuilder() = default;
+    GraphBuilder(const size_t num_nodes, const size_t num_edges)
+    {
+        m_nodes.reserve(num_nodes);
+        m_edges.reserve(num_edges);
+    }
+
+    /**
+     * @brief Add a node to the GraphBuilder.
+     *
+     * The property of the node is constructed from arguments.
+     * @param id Id for the node.
+     * @tparam args Additional arguments for node construction.
+     */
+    template <class... Args>
+    void add_node(int id, Args&&... args)
+    {
+        m_nodes.emplace_back(id, std::forward<Args>(args)...);
+    }
+
+    /**
+     * @brief Add an edge to the GraphBuilder.
+     * 
+     * @param start_node_idx Id of start node
+     * @param end_node_idx Id of end node
+     * @tparam args Additional arguments for edge construction
+     */
+    template <class... Args>
+    void add_edge(size_t start_node_idx, size_t end_node_idx, Args&&... args)
+    {
+        assert(m_nodes.size() > start_node_idx && m_nodes.size() > end_node_idx);
+        m_edges.emplace_back(start_node_idx, end_node_idx, std::forward<Args>(args)...);
+    }
+
+    /**
+     * @brief Build the graph from the added nodes and edges.
+     * 
+     * Sorts the edges and optionally removes duplicate edges (same start and end node indices).
+     * @param make_unique If true, duplicate edges are removed. The first added edge is kept!
+     * @return Graph<NodePropertyT, EdgePropertyT> The constructed graph.
+     * @tparam NodeProperty The type of the node property.
+     * @tparam EdgeProperty The type of the edge property.
+     */
+    Graph<NodeProperty, EdgeProperty> build(bool make_unique = false)
+    {
+        sort_edges();
+        if (make_unique) {
+            remove_duplicate_edges();
+        }
+        Graph<NodeProperty, EdgeProperty> graph(std::move(m_nodes), std::move(m_edges));
+        return graph;
+    }
+
+private:
+    /**
+     * @brief Sort the edge vector of a graph.
+     */
+    void sort_edges()
+    {
+        std::stable_sort(m_edges.begin(), m_edges.end(), [](auto&& e1, auto&& e2) {
+            return e1.start_node_idx == e2.start_node_idx ? e1.end_node_idx < e2.end_node_idx
+                                                          : e1.start_node_idx < e2.start_node_idx;
+        });
+    }
+
+    /**
+     * @brief Remove duplicate edges from a sorted edge vector.
+     * 
+     * Copies all the unique edges to a new vector and replaces the original edge vector with it. Unique means that
+     * the start and end node indices are unique. Other edge properties are not checked and may get lost. Only the first 
+     * edge in the vector is kept.
+     */
+    void remove_duplicate_edges()
+    {
+        std::vector<Edge<EdgePropertyT>> unique_edges;
+        unique_edges.reserve(m_edges.size());
+        std::ranges::unique_copy(m_edges, std::back_inserter(unique_edges), [](auto&& e1, auto&& e2) {
+            return e1.start_node_idx == e2.start_node_idx && e1.end_node_idx == e2.end_node_idx;
+        });
+        m_edges = std::move(unique_edges);
+    }
+
+private:
+    std::vector<Node<NodePropertyT>> m_nodes;
+    std::vector<Edge<EdgePropertyT>> m_edges;
+};
+
 } // namespace mio
 
 #endif //GRAPH_H

cpp/memilio/utils/stl_util.h

@@ -65,23 +65,24 @@ inline std::ostream& set_ostream_format(std::ostream& out, size_t width, size_t
  * @return iterator to inserted or replaced item in vec
  */
 template <typename T, typename Pred>
-typename std::vector<T>::iterator insert_sorted_replace(std::vector<T>& vec, T const& item, Pred pred)
+void insert_sorted_replace(std::vector<T>& vec, T const& item, Pred pred)
 {
     auto bounds = std::equal_range(begin(vec), end(vec), item, pred);
     auto lb     = bounds.first;
     auto ub     = bounds.second;
     assert(ub - lb <= 1); //input vector contains at most one item that is equal to the new item
     if (ub - lb == 1) {
         *lb = item;
-        return lb;
+        return;
     }
     else {
-        return vec.insert(lb, item);
+        vec.insert(lb, item);
+        return;
     }
 }
 
 template <typename T>
-typename std::vector<T>::iterator insert_sorted_replace(std::vector<T>& vec, T const& item)
+void insert_sorted_replace(std::vector<T>& vec, T const& item)
 {
     return insert_sorted_replace(vec, item, std::less<T>());
 }

cpp/tests/test_graph.cpp

@@ -325,7 +325,7 @@ TEST(TestGraph, set_edges_saving_edges)
     EXPECT_EQ(indices_edge1, indices_save_edges);
 }
 
-TEST(TestGraph, ot_edges)
+TEST(TestGraph, out_edges)
 {
     mio::Graph<int, int> g;
     g.add_node(0);
@@ -344,6 +344,42 @@ TEST(TestGraph, ot_edges)
     EXPECT_THAT(g.out_edges(1), testing::ElementsAreArray(v1));
 }
 
+TEST(TestGraphBuilder, Build)
+{
+    mio::GraphBuilder<int, int> builder(3, 3);
+    builder.add_node(0, 100);
+    builder.add_node(1, 100);
+    builder.add_node(2, 100);
+    builder.add_edge(0, 1, 100);
+    builder.add_edge(2, 1, 100);
+    builder.add_edge(1, 2, 100);
+
+    auto g = builder.build();
+
+    EXPECT_EQ(g.nodes().size(), 3);
+    EXPECT_EQ(g.edges().size(), 3);
+}
+
+TEST(TestGraphBuilder, Build_unique)
+{
+    mio::GraphBuilder<int, int> builder;
+    builder.add_node(0, 100);
+    builder.add_node(1, 100);
+    builder.add_node(2, 100);
+    builder.add_edge(1, 2, 100);
+    builder.add_edge(0, 1, 100);
+    builder.add_edge(2, 1, 100);
+    builder.add_edge(1, 2, 200);
+
+    auto g = builder.build(true);
+
+    EXPECT_EQ(g.nodes().size(), 3);
+    EXPECT_EQ(g.edges().size(), 3);
+    for (const auto& e : g.edges()) {
+        EXPECT_EQ(e.property, 100);
+    }
+}
+
 namespace
 {
 

cpp/tests/test_stl_util.cpp

@@ -109,22 +109,6 @@ TEST(TestInsertSortedReplace, normal)
     EXPECT_THAT(v, testing::ElementsAre(1, 2, 5, 6, 7));
 }
 
-TEST(TestInsertSortedReplace, returnsValidIterator)
-{
-    std::vector<int> v;
-    int x;
-
-    //There is no GTEST_NO_DEATH macro so we just let the test crash.
-    //If this test crashes, the function does not return a valid iterator.
-    //Dereferencing an invalid iterator is undefined behavior so the test
-    //may behave unexpectedly (pass, fail, or something else) if the iterator is invalid.
-    x = *mio::insert_sorted_replace(v, 5);
-    x = *mio::insert_sorted_replace(v, 1);
-    x = *mio::insert_sorted_replace(v, 4);
-    x = *mio::insert_sorted_replace(v, 7);
-    ASSERT_EQ(x, 7);
-}
-
 TEST(TestInsertSortedReplace, reverse)
 {
     std::vector<int> v = {5};

docs/source/cpp/graph_metapop.rst

@@ -177,6 +177,36 @@ The following steps detail how to configure and execute a graph simulation:
         graph.add_edge(0, 1, std::move(transition_rates));
         graph.add_edge(1, 0, std::move(transition_rates));
 
+.. dropdown:: :fa:`gears` Working with large graphs
+
+    When working with very large graphs, i.e. starting from a few thousand edges, it will be faster to not use the standard ``add_edge`` function.
+    For this case, we provide a ``GraphBuilder``. There you can add all edges without any checks and the edges will be sorted when the graph is generated:
+
+    .. code-block:: cpp
+
+        mio::GraphBuilder<mio::SimulationNode<mio::Simulation<double, mio::osecir::Model<double>>>, mio::MobilityEdgeStochastic> builder;
+        builder.add_node(1001, model_group1, t0);
+        builder.add_node(1002, model_group2, to);
+        builder.add_edge(0, 1, std::move(transition_rates));
+        builder.add_edge(1, 0, std::move(transition_rates));
+        auto graph = builder.build();
+
+
+    Usually, there should be no duplicate edges. If this is not certain, the ``GraphBuilder`` can also remove duplicates, based on the start and end node. 
+    The parameters in the edge will not be compared. In this case it will only keep the first edge that was inserted:
+
+    .. code-block:: cpp
+
+        mio::GraphBuilder<Int, Int> builder;
+        builder.add_node(1001, 100);
+        builder.add_node(1002, 100);
+        builder.add_edge(0, 1, 100);
+        builder.add_edge(1, 0, 100);
+        builder.add_edge(0, 1, 200);
+        auto graph = builder.build(true);
+        // graph contains the edges (0, 1, 100) and (1, 0, 100)
+
+
 5. **Initialize and Advance the Mobility Simulation:**
 
    With the graph constructed, initialize the simulation with the starting time and time step. Then, advance the simulation until the final time :math:`t_{max}`.

pycode/memilio-simulation/memilio/simulation/bindings/mobility/metapopulation_mobility_instant.h

@@ -37,8 +37,8 @@ void bind_MobilityGraph(pybind11::module_& m, std::string const& name)
         .def(pybind11::init<>())
         .def(
             "add_node",
-            [](G& self, int id, const typename Simulation::Model& p, double t0, double dt) -> auto& {
-                return self.add_node(id, p, t0, dt);
+            [](G& self, int id, const typename Simulation::Model& p, double t0, double dt) -> void {
+                self.add_node(id, p, t0, dt);
             },
             pybind11::arg("id"), pybind11::arg("model"), pybind11::arg("t0") = 0.0, pybind11::arg("dt") = 0.1,
             pybind11::return_value_policy::reference_internal)