Implement Parallel Spatial-Ordering/Cut Selection

Extends explanation for recursive bisection ids
Cleans up Bisection State
Removes license boilerplate from partitioner config
Sorts Spatially and picks Sources and Sinks
Uses sets for sources and sinks for now; see how large they will get
Runs n cuts in parallel changing the slope and uses the best
Clarifies balance <-> ratio naming

include/partition/dinic_max_flow.hpp

@@ -31,7 +31,11 @@ namespace partition
 class DinicMaxFlow
 {
   public:
-    using PartitionResult = std::vector<bool>;
+    using PartitionResult = struct
+    {
+        std::size_t num_edges;
+        std::vector<bool> flags;
+    };
     using SourceSinkNodes = std::set<NodeID>;
     using LevelGraph = std::unordered_map<NodeID, std::uint32_t>;
     using FlowEdges = std::unordered_set<std::pair<NodeID, NodeID>>;

include/partition/graph_view.hpp

@@ -81,6 +81,9 @@ class GraphView
 
     NodeID GetTarget(const EdgeID eid) const;
 
+    const BisectionNode &GetNode(const NodeID nid) const;
+    const BisectionEdge &GetEdge(const EdgeID eid) const;
+
   private:
     const BisectionGraph &bisection_graph;
     const RecursiveBisectionState &bisection_state;

include/partition/inertial_flow.hpp

@@ -2,6 +2,8 @@
 #define OSRM_PARTITION_INERTIAL_FLOW_HPP_
 
 #include "partition/graph_view.hpp"
+
+#include <unordered_set>
 #include <vector>
 
 namespace osrm
@@ -15,7 +17,24 @@ class InertialFlow
     InertialFlow(const GraphView &view);
 
     std::vector<bool> ComputePartition(const double balance, const double source_sink_rate);
+
   private:
+    // Spatially ordered sources and sink ids.
+    // The node ids refer to nodes in the GraphView.
+    struct SpatialOrder
+    {
+        std::unordered_set<NodeID> sources;
+        std::unordered_set<NodeID> sinks;
+    };
+
+    // Creates a spatial order of n * sources "first" and n * sink "last" node ids.
+    // The slope determines the spatial order for sorting node coordinates.
+    SpatialOrder MakeSpatialOrder(double ratio, double slope) const;
+
+    // Makes n cuts with different spatial orders and returns the best.
+    MinCut bestMinCut(std::size_t n, double ratio) const;
+
+    // The subgraph to partition into two parts.
     const GraphView &view;
 };
 

include/partition/partition_config.hpp

@@ -13,7 +13,7 @@ namespace partition
 
 struct PartitionConfig
 {
-    PartitionConfig() noexcept : requested_num_threads(0) {}
+    PartitionConfig() : requested_num_threads(0) {}
 
     void UseDefaults()
     {

include/partition/recursive_bisection_state.hpp

@@ -41,10 +41,26 @@ namespace partition
 // 
 // bisection-ids: [00,10,01,10,00,11,01,11,00,10]
 
+/* Written out in a recursive tree form:
+
+ ids:     [0,1,2,3,4,5,6,7,8,9]
+ mask:    [0,1,0,1,0,1,0,1,0,1]
+              /          \
+ ids:    [0,2,4,6,8] [1,3,5,7,9]
+ mask:   [0,1,0,1,0] [0,0,1,1,0]
+           /    \      /    \
+ ids:  [0,4,8] [2,6] [1,3,9] [5,7]
+
+ The bisection ids then trace the path (left: 0, right: 1) through the tree:
+
+ ids:     [0,  1,  2,  3,  4,  5,  6,  7,  8,  9 ]
+ path:    [00, 10, 01, 10, 00, 11, 01, 11, 00, 10]
+
+*/
 class RecursiveBisectionState
 {
   public:
-    // the ID in the partition arr
+    // The ID in the partition array
     using BisectionID = std::uint32_t;
     using IDIterator = std::vector<NodeID>::const_iterator;
 
@@ -53,7 +69,8 @@ class RecursiveBisectionState
 
     BisectionID GetBisectionID(const NodeID nid) const;
 
-    // returns the center of the bisection
+    // Bisects the node id array's sub-range based on the partition mask.
+    // Returns: partition point of the bisection: iterator to the second group's first element.
     IDIterator ApplyBisection(const IDIterator begin,
                               const IDIterator end,
                               const std::vector<bool> &partition);

include/partition/reorder_first_last.hpp

@@ -0,0 +1,55 @@
+#ifndef OSRM_REORDER_FIRST_LAST_HPP
+#define OSRM_REORDER_FIRST_LAST_HPP
+
+#include <boost/assert.hpp>
+
+#include <algorithm>
+#include <cstddef>
+#include <iterator>
+
+namespace osrm
+{
+namespace partition
+{
+
+// Reorders the first n elements in the range to satisfy the comparator,
+// and the last n elements to satisfy the comparator with arguments flipped.
+// Note: no guarantees to the element's ordering inside the reordered ranges.
+template <typename RandomIt, typename Comparator>
+void reorderFirstLast(RandomIt first, RandomIt last, std::size_t n, Comparator comp)
+{
+    BOOST_ASSERT_MSG(n <= (last - first) / std::size_t{2}, "overlapping subranges not allowed");
+
+    if (n == 0 or (last - first < 2))
+        return;
+
+    // Reorder first n: guarantees that the predicate holds for the first elements.
+    std::nth_element(first, first + (n - 1), last, comp);
+
+    // Reorder last n: guarantees that the flipped predicate holds for the last k elements.
+    // We reorder from the end backwards up to the end of the already reordered range.
+    // We can not use std::not2, since then e.g. std::less<> would lose its irreflexive
+    // requirements.
+    std::reverse_iterator<RandomIt> rfirst{last}, rlast{first + n};
+
+    const auto flipped = [](auto fn) {
+        return [fn](auto &&lhs, auto &&rhs) {
+            return fn(std::forward<decltype(lhs)>(rhs), std::forward<decltype(rhs)>(lhs));
+        };
+    };
+
+    std::nth_element(rfirst, rfirst + (n - 1), rlast, flipped(comp));
+}
+
+template <typename RandomAccessRange, typename Compare>
+void reorderFirstLast(RandomAccessRange &rng, std::size_t n, Compare comp)
+{
+    using std::begin;
+    using std::end;
+    return reorderFirstLast(begin(rng), end(rng), n, comp);
+}
+
+} // ns partition
+} // ns osrm
+
+#endif

src/partition/graph_view.cpp

@@ -71,5 +71,15 @@ NodeID GraphView::GetTarget(const EdgeID eid) const
     return bisection_graph.GetTarget(eid);
 }
 
+const BisectionNode &GraphView::GetNode(const NodeID nid) const
+{
+    return bisection_graph.GetNode(nid);
+}
+
+const BisectionEdge &GraphView::GetEdge(const EdgeID eid) const
+{
+    return bisection_graph.GetEdge(eid);
+}
+
 } // namespace partition
 } // namespace osrm

src/partition/inertial_flow.cpp

@@ -1,10 +1,18 @@
 #include "partition/inertial_flow.hpp"
+#include "partition/bisection_graph.hpp"
 #include "partition/dinic_max_flow.hpp"
+#include "partition/reorder_first_last.hpp"
 
-#include <cstddef>
-#include <set>
+#include <algorithm>
 #include <cmath>
-#include <bitset>
+#include <cstddef>
+#include <iterator>
+#include <mutex>
+#include <set>
+#include <utility>
+
+#include <tbb/blocked_range.h>
+#include <tbb/parallel_for.h>
 
 namespace osrm
 {
@@ -13,30 +21,102 @@ namespace partition
 
 InertialFlow::InertialFlow(const GraphView &view_) : view(view_) {}
 
-std::vector<bool> InertialFlow::ComputePartition(const double balance, const double source_sink_rate)
+std::vector<bool> InertialFlow::ComputePartition(const double balance,
+                                                 const double source_sink_rate)
 {
-    std::set<NodeID> sources;
-    std::set<NodeID> sinks;
+    auto cut = bestMinCut(10 /* should be taken from outside */, source_sink_rate);
 
-    std::size_t count = std::ceil(source_sink_rate * view.NumberOfNodes());
-
-    auto itr = view.Begin();
-    auto itr_end = view.End();
-    while(count--)
-    {
-        --itr_end;
-        sinks.insert(*itr_end);
-        sources.insert(*itr);
-        ++itr;
-    }
-
-    std::cout << "Running Flow" << std::endl;
-    auto result = DinicMaxFlow()(view,sources,sinks);
     std::cout << "Partition: ";
-    for( auto b : result)
+    for (auto b : cut.flags)
         std::cout << b;
     std::cout << std::endl;
-    return result;
+
+    return cut.flags;
+}
+
+InertialFlow::SpatialOrder InertialFlow::MakeSpatialOrder(const double ratio,
+                                                          const double slope) const
+{
+    struct NodeWithCoordinate
+    {
+        NodeWithCoordinate(NodeID nid_, util::Coordinate coordinate_)
+            : nid{nid_}, coordinate{std::move(coordinate_)}
+        {
+        }
+
+        NodeID nid;
+        util::Coordinate coordinate;
+    };
+
+    using Embedding = std::vector<NodeWithCoordinate>;
+
+    Embedding embedding;
+    embedding.reserve(view.NumberOfNodes());
+
+    std::transform(view.Begin(), view.End(), std::back_inserter(embedding), [&](const auto nid) {
+        return NodeWithCoordinate{nid, view.GetNode(nid).coordinate};
+    });
+
+    const auto project = [slope](const auto &each) {
+        auto lon = static_cast<std::int32_t>(each.coordinate.lon);
+        auto lat = static_cast<std::int32_t>(each.coordinate.lat);
+
+        return slope * lon + (1. - std::fabs(slope)) * lat;
+    };
+
+    const auto spatially = [&](const auto &lhs, const auto &rhs) {
+        return project(lhs) < project(rhs);
+    };
+
+    const std::size_t n = ratio * embedding.size();
+
+    reorderFirstLast(embedding, n, spatially);
+
+    InertialFlow::SpatialOrder order;
+
+    order.sources.reserve(n);
+    order.sinks.reserve(n);
+
+    for (auto it = begin(embedding), last = begin(embedding) + n; it != last; ++it)
+        order.sources.insert(it->nid);
+
+    for (auto it = end(embedding) - n, last = end(embedding); it != last; ++it)
+        order.sinks.insert(it->nid);
+
+    return order;
+}
+
+MinCut InertialFlow::bestMinCut(const std::size_t n, const double ratio) const
+{
+    auto base = MakeSpatialOrder(ratio, -1.);
+    auto best = DinicMaxFlow()(view, base.sources, base.sinks);
+
+    std::mutex lock;
+
+    tbb::blocked_range<std::size_t> range{1, n + 1};
+
+    tbb::parallel_for(range, [&, this](const auto &chunk) {
+        for (auto round = chunk.begin(), end = chunk.end(); round != end; ++round)
+        {
+            const auto slope = -1. + round * (2. / n);
+
+            auto order = this->MakeSpatialOrder(ratio, slope);
+            auto cut = Dinic(view, order.sources, order.sinks);
+            auto cut = DinicMaxFlow()(view, order.sources, order.sinks);
+
+            {
+                std::lock_guard<std::mutex> guard{lock};
+
+                // Swap to keep the destruction of the old object outside of critical section.
+                if (cut.num_edges < best.num_edges)
+                    std::swap(best, cut);
+            }
+
+            // cut gets destroyed here
+        }
+    });
+
+    return best;
 }
 
 } // namespace partition

src/partition/recursive_bisection_state.cpp

@@ -2,9 +2,9 @@
 
 #include <numeric>
 
-//TODO remove
-#include <iostream>
+// TODO remove
 #include <bitset>
+#include <iostream>
 
 namespace osrm
 {
@@ -15,20 +15,20 @@ RecursiveBisectionState::RecursiveBisectionState(const BisectionGraph &bisection
     : bisection_graph(bisection_graph_)
 {
     id_array.resize(bisection_graph.GetNumberOfNodes());
-    std::iota(id_array.begin(), id_array.end(), 0);
-    bisection_ids.resize(bisection_graph.GetNumberOfNodes(), 0);
+    std::iota(id_array.begin(), id_array.end(), NodeID{0});
+    bisection_ids.resize(bisection_graph.GetNumberOfNodes(), BisectionID{0});
 }
 
 RecursiveBisectionState::~RecursiveBisectionState()
 {
     std::cout << "Internal Result\n";
     std::cout << "IDArray:";
-    for( auto id : id_array )
+    for (auto id : id_array)
         std::cout << " " << id;
     std::cout << std::endl;
 
     std::cout << "BisectionIDs:";
-    for( auto id : bisection_ids)
+    for (auto id : bisection_ids)
         std::cout << " " << (std::bitset<4>(id));
 
     std::cout << std::endl;
@@ -36,12 +36,12 @@ RecursiveBisectionState::~RecursiveBisectionState()
 
 const RecursiveBisectionState::IDIterator RecursiveBisectionState::Begin() const
 {
-    return id_array.begin();
+    return id_array.cbegin();
 }
 
 const RecursiveBisectionState::IDIterator RecursiveBisectionState::End() const
 {
-    return id_array.end();
+    return id_array.cend();
 }
 
 RecursiveBisectionState::BisectionID RecursiveBisectionState::GetBisectionID(const NodeID nid) const
@@ -59,10 +59,15 @@ RecursiveBisectionState::IDIterator RecursiveBisectionState::ApplyBisection(
         bisection_ids[*itr] |= partition[std::distance(begin, itr)];
     }
 
-    // keep items with `0` as partition id to the left, move other to the right
-    return std::stable_partition(id_array.begin() + std::distance(id_array.cbegin(), begin),
-                                 id_array.begin() + std::distance(id_array.cbegin(), end),
-                                 [this](const auto nid) { return 0 == (bisection_ids[nid] & 1); });
+    auto first = id_array.begin() + std::distance(id_array.cbegin(), begin);
+    auto last = id_array.begin() + std::distance(id_array.cbegin(), end);
+
+    // Keep items with `0` as partition id to the left, move other to the right
+    auto by_last_bit = [this](const auto nid) {
+        return BisectionID{0} == (bisection_ids[nid] & 1);
+    };
+
+    return std::stable_partition(first, last, by_last_bit);
 }
 
 } // namespace partition