Implement Dinic Algorithm for MaximumFlow/MinimumCut

include/partition/bisection_graph.hpp

@@ -55,32 +55,34 @@ template <typename InputEdge>
 std::vector<BisectionNode> computeNodes(const std::vector<util::Coordinate> &coordinates,
                                         const std::vector<InputEdge> &edges)
 {
-    std::vector<BisectionNode> result(coordinates.size() + 1 /*sentinel*/);
+    std::vector<BisectionNode> result;
+    result.reserve(coordinates.size() + 1 /*sentinel*/);
 
-    // stateful transform, counting node ids and moving the edge itr forward
-    const auto coordinate_to_bisection_node =
-        [ edge_itr = edges.begin(), node_id = 0u, &edges ](
-            const util::Coordinate coordinate) mutable->BisectionNode
-    {
-        const auto edges_of_node = edge_itr;
-
-        // count all edges with this source
+    // find the end of edges that belong to node_id
+    const auto advance_edge_itr = [&edges](const std::size_t node_id, auto edge_itr) {
         while (edge_itr != edges.end() && edge_itr->source == node_id)
             ++edge_itr;
-
-        // go to the next node
-        ++node_id;
-
-        return {static_cast<std::size_t>(std::distance(edges.begin(), edges_of_node)), coordinate};
+        return edge_itr;
     };
 
-    std::transform(
-        begin(coordinates), end(coordinates), begin(result), coordinate_to_bisection_node);
+    // create a bisection node, requires the ID of the node as well as the lower bound to its edges
+    const auto make_bisection_node = [&edges, &coordinates](const std::size_t node_id,
+                                                            const auto edge_itr) -> BisectionNode {
+        return {static_cast<std::size_t>(std::distance(edges.begin(), edge_itr)),
+                coordinates[node_id]};
+    };
+
+    auto edge_itr = edges.begin();
+    for (std::size_t node_id = 0; node_id < coordinates.size(); ++node_id)
+    {
+        result.emplace_back(make_bisection_node(node_id,edge_itr));
+        edge_itr = advance_edge_itr(node_id,edge_itr);
+    }
 
     auto null_island = util::Coordinate(util::FloatLongitude{0.0}, util::FloatLatitude{0.0});
     auto sentinel = BisectionNode{edges.size(), std::move(null_island)};
 
-    result.back() = std::move(sentinel);
+    result.emplace_back(std::move(sentinel));
 
     return result;
 }

include/partition/dinic_max_flow.hpp

@@ -0,0 +1,68 @@
+#ifndef OSRM_PARTITION_DINIC_MAX_FLOW_HPP_
+#define OSRM_PARTITION_DINIC_MAX_FLOW_HPP_
+
+#include "partition/graph_view.hpp"
+
+#include <cstdint>
+#include <functional>
+#include <set>
+#include <unordered_map>
+#include <unordered_set>
+#include <utility>
+#include <vector>
+
+namespace std
+{
+template <> struct hash<std::pair<NodeID, NodeID>>
+{
+    std::size_t operator()(const std::pair<NodeID, NodeID> &flow_edge) const
+    {
+        std::size_t combined = (static_cast<std::size_t>(flow_edge.first) << 32) | flow_edge.second;
+        return std::hash<std::size_t>()(combined);
+    }
+};
+}
+
+namespace osrm
+{
+namespace partition
+{
+
+class DinicMaxFlow
+{
+  public:
+    using PartitionResult = std::vector<bool>;
+    using SourceSinkNodes = std::set<NodeID>;
+    using LevelGraph = std::unordered_map<NodeID, std::uint32_t>;
+    using FlowEdges = std::unordered_set<std::pair<NodeID, NodeID>>;
+
+    PartitionResult operator()(const GraphView &view,
+                               const SourceSinkNodes &sink_nodes,
+                               const SourceSinkNodes &source_nodes) const;
+
+  private:
+    LevelGraph ComputeLevelGraph(const GraphView &view,
+                                 const SourceSinkNodes &source_nodes,
+                                 const FlowEdges &flow) const;
+
+    void AugmentFlow(FlowEdges &flow,
+                     const GraphView &view,
+                     const SourceSinkNodes &source_nodes,
+                     const SourceSinkNodes &sink_nodes,
+                     const LevelGraph &levels) const;
+
+    bool findPath(const NodeID from,
+                  std::vector<NodeID> &path,
+                  const GraphView &view,
+                  const LevelGraph &levels,
+                  const FlowEdges &flow,
+                  const SourceSinkNodes &sink_nodes) const;
+};
+
+} // namespace partition
+} // namespace osrm
+
+// Implementation of Dinics [1] algorithm for max-flow/min-cut.
+// [1] https://www.cs.bgu.ac.il/~dinitz/D70.pdf
+
+#endif // OSRM_PARTITION_DINIC_MAX_FLOW_HPP_

include/partition/graph_view.hpp

@@ -49,7 +49,7 @@ class EdgeIDIterator : public boost::iterator_facade<EdgeIDIterator,
     void decrement() { --position; }
     void advance(difference_type offset) { position += offset; }
     bool equal(const EdgeIDIterator &other) const { return position == other.position; }
-    const reference dereference() const { return position; }
+    reference dereference() const { return position; }
     difference_type distance_to(const EdgeIDIterator &other) const
     {
         return static_cast<difference_type>(other.position - position);
@@ -70,7 +70,7 @@ class GraphView
               const RecursiveBisectionState::IDIterator begin,
               const RecursiveBisectionState::IDIterator end);
 
-    // Number of nodes _in this sub-graph_.
+    // Number of nodes _in this sub-graph.
     std::size_t NumberOfNodes() const;
 
     RecursiveBisectionState::IDIterator Begin() const;
@@ -79,6 +79,8 @@ class GraphView
     EdgeIterator EdgeBegin(const NodeID nid) const;
     EdgeIterator EdgeEnd(const NodeID nid) const;
 
+    NodeID GetTarget(const EdgeID eid) const;
+
   private:
     const BisectionGraph &bisection_graph;
     const RecursiveBisectionState &bisection_state;

src/partition/dinic_max_flow.cpp

@@ -0,0 +1,177 @@
+#include "partition/dinic_max_flow.hpp"
+
+#include <queue>
+
+namespace osrm
+{
+namespace partition
+{
+
+DinicMaxFlow::PartitionResult DinicMaxFlow::operator()(const GraphView &view,
+                                                       const SourceSinkNodes &source_nodes,
+                                                       const SourceSinkNodes &sink_nodes) const
+{
+    FlowEdges flow;
+    do
+    {
+        const auto levels = ComputeLevelGraph(view, source_nodes, flow);
+
+        // check if the sink can be reached from the source
+        const auto separated =
+            std::find_if(sink_nodes.begin(), sink_nodes.end(), [&levels](const auto node) {
+                return levels.find(node) != levels.end();
+            }) == sink_nodes.end();
+
+        // no further elements can be found
+        if (separated)
+        {
+            // all elements within `levels` are on the source side
+            PartitionResult result(view.NumberOfNodes(), true);
+            for (auto itr = view.Begin(); itr != view.End(); ++itr)
+            {
+                if (levels.find(*itr) != levels.end())
+                    result[std::distance(view.Begin(), itr)] = false;
+            }
+            return result;
+        }
+
+        AugmentFlow(flow, view, source_nodes, sink_nodes, levels);
+    } while (true);
+}
+
+DinicMaxFlow::LevelGraph DinicMaxFlow::ComputeLevelGraph(const GraphView &view,
+                                                         const SourceSinkNodes &source_nodes,
+                                                         const FlowEdges &flow) const
+{
+    LevelGraph levels;
+    std::queue<NodeID> level_queue;
+
+    for (const auto node : source_nodes)
+    {
+        levels[node] = 0;
+        level_queue.push(node);
+    }
+
+    // check if there is flow present on an edge
+    const auto has_flow = [&](const NodeID from, const NodeID to) {
+        return flow.find(std::make_pair(from, to)) != flow.end();
+    };
+
+    // perform a relaxation step in the BFS algorithm
+    const auto relax_node = [&](const NodeID node_id, const std::uint32_t level) {
+        for (auto itr = view.EdgeBegin(node_id); itr != view.EdgeEnd(node_id); ++itr)
+        {
+            const auto target = view.GetTarget(*itr);
+
+            // don't relax edges with flow on them
+            if (has_flow(node_id, target))
+                continue;
+
+            // don't go back, only follow edges to new nodes
+            if (levels.find(target) == levels.end())
+            {
+                level_queue.push(target);
+                levels[target] = level;
+            }
+        }
+    };
+
+    // compute the levels of level graph using BFS
+    for (std::uint32_t level = 1; !level_queue.empty(); ++level)
+    {
+        // run through the current level
+        auto steps = level_queue.size();
+        while (steps--)
+        {
+            relax_node(level_queue.front(), level);
+            level_queue.pop();
+        }
+    }
+
+    return levels;
+}
+
+void DinicMaxFlow::AugmentFlow(FlowEdges &flow,
+                               const GraphView &view,
+                               const SourceSinkNodes &source_nodes,
+                               const SourceSinkNodes &sink_nodes,
+                               const LevelGraph &levels) const
+{
+    const auto has_flow = [&](const auto from, const auto to) {
+        return flow.find(std::make_pair(from, to)) != flow.end();
+    };
+
+    // find a path and augment the flow along its edges
+    // do a dfs on the level graph, augment all paths that are found
+    const auto find_and_augment_path = [&](const NodeID source) {
+        std::vector<NodeID> path;
+        auto has_path = findPath(source, path, view, levels, flow, sink_nodes);
+
+        if (has_path)
+        {
+            NodeID last = path.back();
+            path.pop_back();
+            // augment the flow graph with the flow
+            while (!path.empty())
+            {
+                flow.insert(std::make_pair(path.back(), last));
+
+                // update the residual capacities correctly
+                const auto reverse_flow = flow.find(std::make_pair(last, path.back()));
+                if (reverse_flow != flow.end())
+                    flow.erase(reverse_flow);
+
+                last = path.back();
+                path.pop_back();
+            }
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    };
+
+    // find and augment the blocking flow
+    for (auto source : source_nodes)
+    {
+        while (find_and_augment_path(source))
+        { /*augment*/
+        }
+    }
+}
+
+bool DinicMaxFlow::findPath(const NodeID node_id,
+                            std::vector<NodeID> &path,
+                            const GraphView &view,
+                            const LevelGraph &levels,
+                            const FlowEdges &flow,
+                            const SourceSinkNodes &sink_nodes) const
+{
+    path.push_back(node_id);
+    if (sink_nodes.find(node_id) != sink_nodes.end())
+        return true;
+
+    const auto node_level = levels.find(node_id)->second;
+    for (auto itr = view.EdgeBegin(node_id); itr != view.EdgeEnd(node_id); ++itr)
+    {
+        const auto target = view.GetTarget(*itr);
+
+        // don't relax edges with flow on them
+        if (flow.find(std::make_pair(node_id, target)) != flow.end())
+            continue;
+
+        // don't go back, only follow edges to new nodes
+        const auto level = levels.find(target)->second;
+        if( level != node_level + 1 )
+            continue;
+
+        if( findPath(target,path,view,levels,flow,sink_nodes) )
+            return true;
+    }
+    path.pop_back();
+    return false;
+}
+
+} // namespace partition
+} // namespace osrm

src/partition/graph_view.cpp

@@ -34,7 +34,7 @@ GraphView::GraphView(const BisectionGraph &bisection_graph_,
         std::cout << "Node: " << *itr << std::endl;
         for (auto eitr = EdgeBegin(*itr); eitr != EdgeEnd(*itr); ++eitr)
         {
-            std::cout << "\t" << *eitr << " -> " << bisection_graph.GetTarget(*eitr) << std::endl;
+            std::cout << "\t" << *eitr << " -> " << GetTarget(*eitr) << std::endl;
         }
     }
 }
@@ -66,5 +66,10 @@ GraphView::EdgeIterator GraphView::EdgeEnd(const NodeID nid) const
     return boost::make_filter_iterator(predicate, last, last);
 }
 
+NodeID GraphView::GetTarget(const EdgeID eid) const
+{
+    return bisection_graph.GetTarget(eid);
+}
+
 } // namespace partition
 } // namespace osrm

src/partition/inertial_flow.cpp

@@ -1,4 +1,10 @@
 #include "partition/inertial_flow.hpp"
+#include "partition/dinic_max_flow.hpp"
+
+#include <cstddef>
+#include <set>
+#include <cmath>
+#include <bitset>
 
 namespace osrm
 {
@@ -9,13 +15,28 @@ InertialFlow::InertialFlow(const GraphView &view_) : view(view_) {}
 
 std::vector<bool> InertialFlow::ComputePartition(const double balance, const double source_sink_rate)
 {
-    std::vector<bool> partition(view.NumberOfNodes());
-    std::size_t i = 0;
-    for( auto itr = partition.begin(); itr != partition.end(); ++itr )
+    std::set<NodeID> sources;
+    std::set<NodeID> sinks;
+
+    std::size_t count = std::ceil(source_sink_rate * view.NumberOfNodes());
+
+    auto itr = view.Begin();
+    auto itr_end = view.End();
+    while(count--)
     {
-        *itr = (i++ % 2) != 0;
+        --itr_end;
+        sinks.insert(*itr_end);
+        sources.insert(*itr);
+        ++itr;
     }
-    return partition;
+
+    std::cout << "Running Flow" << std::endl;
+    auto result = DinicMaxFlow()(view,sources,sinks);
+    std::cout << "Partition: ";
+    for( auto b : result)
+        std::cout << b;
+    std::cout << std::endl;
+    return result;
 }
 
 } // namespace partition