Rename namespace partition to partitioner

Rename module partition to partitioner.
This cultivates naming used in existing modules like extractor,
customizer, etc. - noun vs verb (word partition is both though).

src/partitioner/annotated_partition.cpp

@@ -0,0 +1,387 @@
+#include "partitioner/annotated_partition.hpp"
+
+#include <algorithm>
+#include <climits> // for CHAR_BIT
+#include <cstddef>
+#include <cstdint>
+#include <iostream>
+#include <limits>
+#include <map>
+#include <numeric>
+#include <queue>
+#include <string>
+#include <unordered_map>
+
+#include "util/timing_util.hpp"
+
+namespace osrm
+{
+namespace partitioner
+{
+
+namespace
+{
+
+// the shift value needed to access the most significant bit of the bisection ID
+const constexpr auto SHIFT_TO_MSB_BISECTION_ID = sizeof(BisectionID) * CHAR_BIT - 1;
+// an invalid ID for a cell
+const constexpr std::uint32_t INVALID_CELLID = std::numeric_limits<std::uint32_t>::max();
+
+auto masked(const BisectionID id, const std::int32_t level)
+{
+    // special treatment for negative level
+    if (level == -1)
+        return 0u;
+
+    // 0.01.1 with 1 starting at the level+1_th most significant bit (level = 0 -> 01..1)
+    const auto cut_below_level = (1 << (SHIFT_TO_MSB_BISECTION_ID - level)) - 1;
+    const auto mask = std::numeric_limits<BisectionID>::max() ^ cut_below_level;
+    return id & mask;
+}
+
+// create a comparator for a given level
+auto makeCompare(const std::uint32_t level)
+{
+    return [level](const AnnotatedPartition::SizedID lhs, const AnnotatedPartition::SizedID rhs) {
+        return masked(lhs.id, level) < masked(rhs.id, level);
+    };
+}
+
+// build a tree of cells from the IDs present:
+auto leftChild(const BisectionID id_prefix, const std::int32_t /*level*/) { return id_prefix; }
+
+// given the prefix 10.... on level 1 (second level), the the right child would be
+// 101.... on level 2
+auto rightChild(const BisectionID id_prefix, const std::int32_t level)
+{
+    return id_prefix | (1 << (SHIFT_TO_MSB_BISECTION_ID - (level + 1)));
+}
+
+// get the range of all children
+auto getChildrenRange(const std::vector<AnnotatedPartition::SizedID> &implicit_tree,
+                      const BisectionID id_prefix,
+                      const std::int32_t level)
+{
+    AnnotatedPartition::SizedID id = {id_prefix, 0};
+
+    // find all elements of the same prefix as id_prefi
+    auto range =
+        std::equal_range(implicit_tree.begin(), implicit_tree.end(), id, makeCompare(level));
+
+    // don't ever return our sentinel element as included
+    if (range.second == implicit_tree.end())
+        --range.second;
+
+    return range;
+}
+
+auto getCellSize(const std::vector<AnnotatedPartition::SizedID> &implicit_tree,
+                 const BisectionID id_prefix,
+                 const std::uint32_t level)
+{
+    auto range = getChildrenRange(implicit_tree, id_prefix, level);
+    return range.second->count - range.first->count;
+}
+
+bool hasChildren(const std::vector<AnnotatedPartition::SizedID> &implicit_tree,
+                 const BisectionID id_prefix,
+                 const std::uint32_t level)
+{
+    auto range = getChildrenRange(implicit_tree, id_prefix, level);
+    return std::distance(range.first, range.second) > 1;
+}
+
+} // namespace
+
+AnnotatedPartition::AnnotatedPartition(const BisectionGraph &graph,
+                                       const std::vector<BisectionID> &bisection_ids)
+{
+    // create a sorted vector of bisection ids that exist in the network
+    std::vector<SizedID> implicit_tree = [&]() {
+        std::map<BisectionID, SizedID> existing_ids;
+
+        // insert an ID into the sized_id set or increase the count if the element should be already
+        // present in the set of known ids
+        const auto insert_or_augment = [&existing_ids](const BisectionID id) {
+            SizedID sized_id = {id, 1};
+            auto maybe_existing_id = existing_ids.find(id);
+            if (maybe_existing_id == existing_ids.end())
+                existing_ids[id] = sized_id;
+            else
+                maybe_existing_id->second.count++;
+        };
+        std::for_each(bisection_ids.begin(), bisection_ids.end(), insert_or_augment);
+
+        std::vector<SizedID> result;
+        result.resize(existing_ids.size() + 1);
+        std::transform(existing_ids.begin(),
+                       existing_ids.end(),
+                       result.begin(),
+                       [](const auto &pair) { return pair.second; });
+
+        // sentinel
+        result.back() = {std::numeric_limits<BisectionID>::max(), 0};
+
+        return result;
+    }();
+
+    // calculate a prefix sum over all sorted IDs, this allows to get the size of any partition in
+    // the array/level based on the prefix and lower bound on prefixes.
+    // e.g 00,01,10,11 allow to search for (0) (1) to find (00) and (10) as lower bounds. The
+    // difference in count is the size of all cells in the left part of the partition.
+    std::transform(implicit_tree.begin(),
+                   implicit_tree.end(),
+                   implicit_tree.begin(),
+                   [sum = std::size_t{0}](SizedID id) mutable {
+                       const auto new_sum = sum + id.count;
+                       id.count = sum;
+                       sum = new_sum;
+                       return id;
+                   });
+
+    PrintBisection(implicit_tree, graph, bisection_ids);
+    SearchLevels(implicit_tree, graph, bisection_ids);
+}
+
+void AnnotatedPartition::PrintBisection(const std::vector<SizedID> &implicit_tree,
+                                        const BisectionGraph &graph,
+                                        const std::vector<BisectionID> &bisection_ids) const
+{
+    // print some statistics on the bisection tree
+    std::queue<BisectionID> id_queue;
+    id_queue.push(0);
+
+    const auto add_child = [&id_queue, &implicit_tree](const BisectionID prefix,
+                                                       const std::uint32_t level) {
+        const auto child_range = getChildrenRange(implicit_tree, prefix, level);
+        if (std::distance(child_range.first, child_range.second) > 1)
+            id_queue.push(prefix);
+    };
+
+    std::vector<std::pair<BisectionID, std::int32_t>> current_level;
+    for (std::int32_t level = -1; !id_queue.empty(); ++level)
+    {
+        auto level_size = id_queue.size();
+        current_level.clear();
+        while (level_size--)
+        {
+            const auto prefix = id_queue.front();
+            id_queue.pop();
+            if (level == -1 || hasChildren(implicit_tree, prefix, level))
+            {
+                current_level.push_back(
+                    std::pair<BisectionID, std::uint32_t>(leftChild(prefix, level), level + 1));
+                current_level.push_back(
+                    std::pair<BisectionID, std::uint32_t>(rightChild(prefix, level), level + 1));
+            }
+            add_child(leftChild(prefix, level), level);
+            add_child(rightChild(prefix, level), level);
+        }
+        if (!current_level.empty())
+        {
+            const auto cell_ids = ComputeCellIDs(current_level, graph, bisection_ids);
+            const auto stats = AnalyseLevel(graph, cell_ids);
+            stats.logMachinereadable(std::cout, "bisection", level, level == -1);
+        }
+    }
+}
+
+void AnnotatedPartition::SearchLevels(const std::vector<SizedID> &implicit_tree,
+                                      const BisectionGraph &graph,
+                                      const std::vector<BisectionID> &bisection_ids) const
+{
+    std::vector<std::pair<BisectionID, std::int32_t>> current_level;
+
+    // start searching with level 0 at prefix 0
+    current_level.push_back({static_cast<BisectionID>(0), -1});
+    std::int32_t level = -1;
+
+    const auto print_level = [&]() {
+        if (current_level.empty())
+            return;
+
+        const auto cell_ids = ComputeCellIDs(current_level, graph, bisection_ids);
+        const auto stats = AnalyseLevel(graph, cell_ids);
+        stats.logMachinereadable(std::cout, "dfs-balanced", level, level == -1);
+        ++level;
+    };
+
+    std::size_t max_size = 0.5 * graph.NumberOfNodes();
+    std::queue<std::pair<BisectionID, std::int32_t>> id_queue;
+    while (!current_level.empty())
+    {
+        std::size_t total_size = 0;
+        std::size_t count = 0;
+        for (auto element : current_level)
+        {
+            // don't relax final cells
+            if (element.second == -1 || hasChildren(implicit_tree, element.first, element.second))
+            {
+                total_size += getCellSize(
+                    implicit_tree, leftChild(element.first, element.second), element.second + 1);
+                id_queue.push(std::pair<BisectionID, std::uint32_t>(
+                    leftChild(element.first, element.second), element.second + 1));
+
+                total_size += getCellSize(
+                    implicit_tree, rightChild(element.first, element.second), element.second + 1);
+                id_queue.push(std::pair<BisectionID, std::uint32_t>(
+                    rightChild(element.first, element.second), element.second + 1));
+                count += 2;
+            }
+        }
+        auto avg_size = (total_size / static_cast<double>(count));
+
+        current_level.clear();
+
+        const auto relax = [&id_queue, &implicit_tree, avg_size, &current_level](
+            const std::pair<BisectionID, std::uint32_t> &element) {
+            const auto size = getCellSize(implicit_tree, element.first, element.second);
+            if (!hasChildren(implicit_tree, element.first, element.second))
+            {
+                current_level.push_back(element);
+            }
+            else
+            {
+                const auto left = leftChild(element.first, element.second);
+                const auto right = rightChild(element.first, element.second);
+
+                const auto get_penalty = [avg_size](const auto size) {
+                    return std::abs(size - avg_size);
+                };
+
+                if (get_penalty(size) <
+                    0.5 * (get_penalty(getCellSize(implicit_tree, left, element.second + 1)) +
+                           get_penalty(getCellSize(implicit_tree, right, element.second + 1))))
+                {
+                    current_level.push_back(element);
+                }
+                else
+                {
+                    id_queue.push(std::pair<BisectionID, std::uint32_t>(left, element.second + 1));
+                    id_queue.push(std::pair<BisectionID, std::uint32_t>(right, element.second + 1));
+                }
+            }
+        };
+
+        while (!id_queue.empty())
+        {
+            relax(id_queue.front());
+            id_queue.pop();
+        }
+        print_level();
+        max_size *= 0.5;
+    }
+}
+
+AnnotatedPartition::LevelMetrics
+AnnotatedPartition::AnalyseLevel(const BisectionGraph &graph,
+                                 const std::vector<std::uint32_t> &cell_ids) const
+{
+    std::unordered_map<std::uint32_t, std::size_t> cell_sizes;
+    std::unordered_map<std::uint32_t, std::size_t> border_nodes;
+    std::unordered_map<std::uint32_t, std::size_t> border_arcs;
+
+    // compute basic metrics of the level
+    std::size_t border_nodes_total = 0;
+    std::size_t border_arcs_total = 0;
+    std::size_t contained_nodes = 0;
+
+    // only border nodes on the lowest level can be border nodes in general
+    for (const auto &node : graph.Nodes())
+    {
+        const auto cell_id = cell_ids[node.original_id];
+        if (cell_id == INVALID_CELLID)
+            continue;
+
+        ++contained_nodes;
+
+        const auto edge_range = graph.Edges(node);
+        const auto border_arcs_at_node = std::count_if(
+            edge_range.begin(), edge_range.end(), [&cell_id, &cell_ids, &graph](const auto &edge) {
+                const auto target_cell_id = cell_ids[graph.Node(edge.target).original_id];
+                return target_cell_id != cell_id;
+            });
+
+        cell_sizes[cell_id]++;
+        border_arcs[cell_id] += border_arcs_at_node;
+        border_arcs_total += border_arcs_at_node;
+        if (border_arcs_at_node)
+        {
+            border_nodes[cell_id]++;
+            ++border_nodes_total;
+        }
+    }
+
+    const auto by_size = [](const std::pair<std::uint32_t, std::size_t> &lhs,
+                            const std::pair<std::uint32_t, std::size_t> &rhs) {
+        return lhs.second < rhs.second;
+    };
+    const auto max_nodes =
+        border_nodes.empty()
+            ? 0
+            : std::max_element(border_nodes.begin(), border_nodes.end(), by_size)->second;
+    const auto max_arcs =
+        border_arcs.empty()
+            ? 0
+            : std::max_element(border_arcs.begin(), border_arcs.end(), by_size)->second;
+
+    const auto squarded_size = [](const std::size_t accumulated,
+                                  const std::pair<std::uint32_t, std::size_t> &element) {
+        return accumulated + element.second * element.second;
+    };
+
+    const auto memory =
+        4 * std::accumulate(border_arcs.begin(), border_nodes.end(), std::size_t(0), squarded_size);
+
+    std::vector<std::size_t> cell_sizes_vec;
+    cell_sizes_vec.resize(cell_sizes.size());
+    std::transform(cell_sizes.begin(),
+                   cell_sizes.end(),
+                   cell_sizes_vec.begin(),
+                   [](const auto &pair) { return pair.second; });
+    return {border_nodes_total,
+            border_arcs_total,
+            contained_nodes,
+            border_nodes.size(),
+            max_nodes,
+            max_arcs,
+            memory,
+            std::move(cell_sizes_vec)};
+}
+
+std::vector<std::uint32_t>
+AnnotatedPartition::ComputeCellIDs(std::vector<std::pair<BisectionID, std::int32_t>> &prefixes,
+                                   const BisectionGraph &graph,
+                                   const std::vector<BisectionID> &bisection_ids) const
+{
+    std::vector<std::uint32_t> cell_ids(graph.NumberOfNodes(), INVALID_CELLID);
+
+    std::sort(prefixes.begin(), prefixes.end(), [](const auto lhs, const auto rhs) {
+        return lhs.first < rhs.first;
+    });
+
+    for (const auto &node : graph.Nodes())
+    {
+        // find the cell_id of node in the current levels
+        const auto id = bisection_ids[node.original_id];
+
+        const auto is_prefixed_by = [id](const auto &prefix) {
+            return masked(id, prefix.second) == prefix.first;
+        };
+
+        const auto prefix = std::lower_bound(
+            prefixes.begin(), prefixes.end(), id, [&](const auto prefix, const BisectionID id) {
+                return prefix.first < masked(id, prefix.second);
+            });
+
+        if (prefix == prefixes.end())
+            continue;
+
+        if (is_prefixed_by(*prefix))
+            cell_ids[node.original_id] = std::distance(prefixes.begin(), prefix);
+    }
+    return cell_ids;
+}
+
+} // namespace partitioner
+} // namespace osrm

src/partitioner/bisection_graph_view.cpp

@@ -0,0 +1,56 @@
+#include "partitioner/bisection_graph_view.hpp"
+
+#include <iostream>
+#include <iterator>
+
+#include <boost/assert.hpp>
+
+namespace osrm
+{
+namespace partitioner
+{
+
+BisectionGraphView::BisectionGraphView(const BisectionGraph &bisection_graph_)
+    : BisectionGraphView(bisection_graph_, bisection_graph_.CBegin(), bisection_graph_.CEnd())
+{
+}
+
+BisectionGraphView::BisectionGraphView(const BisectionGraph &bisection_graph_,
+                                       const BisectionGraph::ConstNodeIterator begin_,
+                                       const BisectionGraph::ConstNodeIterator end_)
+    : bisection_graph(bisection_graph_), begin(begin_), end(end_)
+{
+}
+
+BisectionGraphView::BisectionGraphView(const BisectionGraphView &other_view,
+                                       const BisectionGraph::ConstNodeIterator begin_,
+                                       const BisectionGraph::ConstNodeIterator end_)
+    : BisectionGraphView(other_view.bisection_graph, begin_, end_)
+{
+}
+
+std::size_t BisectionGraphView::NumberOfNodes() const { return std::distance(begin, end); }
+
+NodeID BisectionGraphView::GetID(const NodeT &node) const
+{
+    const auto node_id = static_cast<NodeID>(&node - &(*begin));
+    BOOST_ASSERT(node_id < NumberOfNodes());
+    return node_id;
+}
+
+BisectionGraph::ConstNodeIterator BisectionGraphView::Begin() const { return begin; }
+
+BisectionGraph::ConstNodeIterator BisectionGraphView::End() const { return end; }
+
+const BisectionGraphView::NodeT &BisectionGraphView::Node(const NodeID nid) const
+{
+    return *(begin + nid);
+}
+
+const BisectionGraphView::EdgeT &BisectionGraphView::Edge(const EdgeID eid) const
+{
+    return bisection_graph.Edge(eid);
+}
+
+} // namespace partitioner
+} // namespace osrm

src/partitioner/bisection_to_partition.cpp

@@ -0,0 +1,144 @@
+#include "partitioner/bisection_to_partition.hpp"
+
+namespace osrm
+{
+namespace partitioner
+{
+
+namespace
+{
+struct CellBisection
+{
+    std::uint32_t begin;
+    std::uint32_t end;
+    std::uint8_t bit;
+    bool tabu; // we will not attempt to split this cell anymore
+};
+static constexpr std::size_t NUM_BISECTION_BITS = sizeof(BisectionID) * CHAR_BIT;
+
+std::vector<std::uint32_t> getLargeCells(const std::size_t max_cell_size,
+                                         const std::vector<CellBisection> &cells)
+{
+    std::vector<std::uint32_t> large_cells;
+
+    for (auto index = 0u; index < cells.size(); ++index)
+    {
+        if (!cells[index].tabu && cells[index].end - cells[index].begin > max_cell_size)
+            large_cells.push_back(index);
+    }
+
+    return large_cells;
+}
+
+Partition cellsToPartition(const std::vector<CellBisection> &cells,
+                           const std::vector<std::uint32_t> &permutation)
+{
+    Partition partition(permutation.size(), INVALID_CELL_ID);
+    CellID cell_id = 0;
+    for (const auto &cell : cells)
+    {
+        std::for_each(permutation.begin() + cell.begin,
+                      permutation.begin() + cell.end,
+                      [&partition, cell_id](const auto node_id) { partition[node_id] = cell_id; });
+        cell_id++;
+    }
+    BOOST_ASSERT(std::find(partition.begin(), partition.end(), INVALID_CELL_ID) == partition.end());
+
+    return partition;
+}
+
+void partitionLevel(const std::vector<BisectionID> &node_to_bisection_id,
+                    std::size_t max_cell_size,
+                    std::vector<std::uint32_t> &permutation,
+                    std::vector<CellBisection> &cells)
+{
+    for (auto large_cells = getLargeCells(max_cell_size, cells); large_cells.size() > 0;
+         large_cells = getLargeCells(max_cell_size, cells))
+    {
+        for (const auto cell_index : large_cells)
+        {
+            auto &cell = cells[cell_index];
+            BOOST_ASSERT(cell.bit < NUM_BISECTION_BITS);
+
+            // Go over all nodes and sum up the bits to determine at which position the first one
+            // bit is
+            BisectionID sum =
+                std::accumulate(permutation.begin() + cell.begin,
+                                permutation.begin() + cell.end,
+                                BisectionID{0},
+                                [&node_to_bisection_id](const BisectionID lhs, const NodeID rhs) {
+                                    return lhs | node_to_bisection_id[rhs];
+                                });
+            // masks all bit strictly higher then cell.bit
+            BOOST_ASSERT(sizeof(unsigned long long) * CHAR_BIT > sizeof(BisectionID) * CHAR_BIT);
+            const BisectionID mask = (1ULL << (cell.bit + 1)) - 1;
+            BOOST_ASSERT(mask == 0 || util::msb(mask) == cell.bit);
+            const auto masked_sum = sum & mask;
+            // we can't split the cell anymore, but it also doesn't conform to the max size
+            // constraint
+            // -> we need to remove it from the optimization
+            if (masked_sum == 0)
+            {
+                cell.tabu = true;
+                continue;
+            }
+            const auto bit = util::msb(masked_sum);
+            // determines if an bisection ID is on the left side of the partition
+            const BisectionID is_left_mask = 1ULL << bit;
+            BOOST_ASSERT(util::msb(is_left_mask) == bit);
+
+            std::uint32_t middle =
+                std::partition(permutation.begin() + cell.begin,
+                               permutation.begin() + cell.end,
+                               [is_left_mask, &node_to_bisection_id](const auto node_id) {
+                                   return node_to_bisection_id[node_id] & is_left_mask;
+                               }) -
+                permutation.begin();
+
+            if (bit > 0)
+                cell.bit = bit - 1;
+            else
+                cell.tabu = true;
+            if (middle != cell.begin && middle != cell.end)
+            {
+                auto old_end = cell.end;
+                cell.end = middle;
+                cells.push_back(
+                    CellBisection{middle, old_end, static_cast<std::uint8_t>(cell.bit), cell.tabu});
+            }
+        }
+    }
+}
+}
+
+// Implements a greedy algorithm that split cells using the bisection until a target cell size is
+// reached
+std::tuple<std::vector<Partition>, std::vector<std::uint32_t>>
+bisectionToPartition(const std::vector<BisectionID> &node_to_bisection_id,
+                     const std::vector<std::size_t> &max_cell_sizes)
+{
+    std::vector<std::uint32_t> permutation(node_to_bisection_id.size());
+    std::iota(permutation.begin(), permutation.end(), 0);
+
+    std::vector<CellBisection> cells;
+    cells.push_back(CellBisection{
+        0, static_cast<std::uint32_t>(node_to_bisection_id.size()), NUM_BISECTION_BITS - 1, false});
+
+    std::vector<Partition> partitions(max_cell_sizes.size());
+    std::vector<std::uint32_t> num_cells(max_cell_sizes.size());
+
+    int level_idx = max_cell_sizes.size() - 1;
+    for (auto max_cell_size : boost::adaptors::reverse(max_cell_sizes))
+    {
+        BOOST_ASSERT(level_idx >= 0);
+        partitionLevel(node_to_bisection_id, max_cell_size, permutation, cells);
+
+        partitions[level_idx] = cellsToPartition(cells, permutation);
+        num_cells[level_idx] = cells.size();
+        level_idx--;
+    }
+
+    return std::make_tuple(std::move(partitions), std::move(num_cells));
+}
+}
+}

src/partitioner/dinic_max_flow.cpp

@@ -0,0 +1,314 @@
+#include "partitioner/dinic_max_flow.hpp"
+#include "util/integer_range.hpp"
+
+#include <algorithm>
+#include <limits>
+#include <numeric>
+#include <queue>
+#include <set>
+#include <stack>
+
+namespace osrm
+{
+namespace partitioner
+{
+
+namespace
+{
+
+const auto constexpr INVALID_LEVEL = std::numeric_limits<DinicMaxFlow::Level>::max();
+
+auto makeHasNeighborNotInCheck(const DinicMaxFlow::SourceSinkNodes &set,
+                               const BisectionGraphView &view)
+{
+    return [&](const NodeID nid) {
+        const auto is_not_contained = [&set](const BisectionEdge &edge) {
+            return set.count(edge.target) == 0;
+        };
+        return view.EndEdges(nid) !=
+               std::find_if(view.BeginEdges(nid), view.EndEdges(nid), is_not_contained);
+    };
+}
+
+} // end namespace
+
+DinicMaxFlow::MinCut DinicMaxFlow::operator()(const BisectionGraphView &view,
+                                              const SourceSinkNodes &source_nodes,
+                                              const SourceSinkNodes &sink_nodes) const
+{
+    BOOST_ASSERT(Validate(view, source_nodes, sink_nodes));
+    // for the inertial flow algorithm, we use quite a large set of nodes as source/sink nodes. Only
+    // a few of them can be part of the process, since they are grouped together. A standard
+    // parameterisation would be 25% sink/source nodes. This already includes 50% of the graph. By
+    // only focussing on a small set on the outside of the source/sink blob, we can save quite some
+    // overhead in initialisation/search cost.
+
+    std::vector<NodeID> border_source_nodes;
+    border_source_nodes.reserve(0.01 * source_nodes.size());
+
+    std::copy_if(source_nodes.begin(),
+                 source_nodes.end(),
+                 std::back_inserter(border_source_nodes),
+                 makeHasNeighborNotInCheck(source_nodes, view));
+
+    std::vector<NodeID> border_sink_nodes;
+    border_sink_nodes.reserve(0.01 * sink_nodes.size());
+    std::copy_if(sink_nodes.begin(),
+                 sink_nodes.end(),
+                 std::back_inserter(border_sink_nodes),
+                 makeHasNeighborNotInCheck(sink_nodes, view));
+
+    // edges in current flow that have capacity
+    // The graph (V,E) contains undirected edges for all (u,v) \in V x V. We describe the flow as a
+    // set of vertices (s,t) with flow set to `true`. Since flow can be either from `s` to `t` or
+    // from `t` to `s`, we can remove `(s,t)` from the flow, if we send flow back the first time,
+    // and insert `(t,s)` only if we send flow again.
+
+    // allocate storage for the flow
+    FlowEdges flow(view.NumberOfNodes());
+    std::size_t flow_value = 0;
+    do
+    {
+        auto levels = ComputeLevelGraph(view, border_source_nodes, source_nodes, sink_nodes, flow);
+
+        // check if the sink can be reached from the source, it's enough to check the border
+        const auto separated = std::find_if(border_sink_nodes.begin(),
+                                            border_sink_nodes.end(),
+                                            [&levels, &view](const auto node) {
+                                                return levels[node] != INVALID_LEVEL;
+                                            }) == border_sink_nodes.end();
+
+        if (!separated)
+        {
+            flow_value += BlockingFlow(flow, levels, view, source_nodes, border_sink_nodes);
+        }
+        else
+        {
+            // mark levels for all sources to not confuse make-cut (due to the border nodes
+            // heuristic)
+            for (auto s : source_nodes)
+                levels[s] = 0;
+            const auto cut = MakeCut(view, levels, flow_value);
+            return cut;
+        }
+    } while (true);
+}
+
+DinicMaxFlow::MinCut DinicMaxFlow::MakeCut(const BisectionGraphView &view,
+                                           const LevelGraph &levels,
+                                           const std::size_t flow_value) const
+{
+    const auto is_valid_level = [](const Level level) { return level != INVALID_LEVEL; };
+
+    // all elements within `levels` are on the source side
+    // This part should opt to find the most balanced cut, which is not necessarily the case right
+    // now. There is potential for optimisation here.
+    std::vector<bool> result(view.NumberOfNodes());
+
+    BOOST_ASSERT(view.NumberOfNodes() == levels.size());
+    std::size_t source_side_count = std::count_if(levels.begin(), levels.end(), is_valid_level);
+    std::transform(levels.begin(), levels.end(), result.begin(), is_valid_level);
+
+    return {source_side_count, flow_value, std::move(result)};
+}
+
+DinicMaxFlow::LevelGraph
+DinicMaxFlow::ComputeLevelGraph(const BisectionGraphView &view,
+                                const std::vector<NodeID> &border_source_nodes,
+                                const SourceSinkNodes &source_nodes,
+                                const SourceSinkNodes &sink_nodes,
+                                const FlowEdges &flow) const
+{
+    LevelGraph levels(view.NumberOfNodes(), INVALID_LEVEL);
+    std::queue<NodeID> level_queue;
+
+    // set the front of the source nodes to zero and add them to the BFS queue. In addition, set all
+    // neighbors to zero as well (which allows direct usage of the levels to see what we visited,
+    // and still don't go back into the hughe set of sources)
+    for (const auto node_id : border_source_nodes)
+    {
+        levels[node_id] = 0;
+        level_queue.push(node_id);
+        for (const auto &edge : view.Edges(node_id))
+            if (source_nodes.count(edge.target))
+                levels[edge.target] = 0;
+    }
+    // check if there is flow present on an edge
+    const auto has_flow = [&](const NodeID from, const NodeID to) {
+        return flow[from].find(to) != flow[from].end();
+    };
+
+    // perform a relaxation step in the BFS algorithm
+    const auto relax_node = [&](const NodeID node_id) {
+        // don't relax sink nodes
+        if (sink_nodes.count(node_id))
+            return;
+
+        const auto level = levels[node_id] + 1;
+        for (const auto &edge : view.Edges(node_id))
+        {
+            const auto target = edge.target;
+            // 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[target] > level)
+            {
+                level_queue.push(target);
+                levels[target] = level;
+            }
+        }
+    };
+
+    // compute the levels of level graph using BFS
+    while (!level_queue.empty())
+    {
+        relax_node(level_queue.front());
+        level_queue.pop();
+    }
+
+    return levels;
+}
+
+std::size_t DinicMaxFlow::BlockingFlow(FlowEdges &flow,
+                                       LevelGraph &levels,
+                                       const BisectionGraphView &view,
+                                       const SourceSinkNodes &source_nodes,
+                                       const std::vector<NodeID> &border_sink_nodes) const
+{
+    // track the number of augmenting paths (which in sum will equal the number of unique border
+    // edges) (since our graph is undirected)
+    std::size_t flow_increase = 0;
+
+    // augment the flow along a path in the level graph
+    const auto augment_flow = [&flow, &view](const std::vector<NodeID> &path) {
+
+        // add/remove flow edges from the current residual graph
+        const auto augment_one = [&flow, &view](const NodeID from, const NodeID to) {
+            // check if there is flow in the opposite direction
+            auto existing_edge = flow[to].find(from);
+            if (existing_edge != flow[to].end())
+                flow[to].erase(existing_edge); // remove flow from reverse edges first
+            else
+                flow[from].insert(to); // only add flow if no opposite flow exists
+
+            // do augmentation on all pairs, never stop early:
+            return false;
+        };
+
+        // augment all adjacent edges
+        std::adjacent_find(path.begin(), path.end(), augment_one);
+    };
+
+    const auto augment_all_paths = [&](const NodeID sink_node_id) {
+        // only augment sinks
+        if (levels[sink_node_id] == INVALID_LEVEL)
+            return;
+
+        while (true)
+        {
+            // as long as there are augmenting paths from the sink, add them
+            const auto path = GetAugmentingPath(levels, sink_node_id, view, flow, source_nodes);
+            if (path.empty())
+                break;
+            else
+            {
+                augment_flow(path);
+                ++flow_increase;
+            }
+        }
+    };
+
+    std::for_each(border_sink_nodes.begin(), border_sink_nodes.end(), augment_all_paths);
+    BOOST_ASSERT(flow_increase > 0);
+    return flow_increase;
+}
+
+// performs a dfs in the level graph, by adjusting levels that don't offer any further paths to
+// INVALID_LEVEL and by following the level graph, this looks at every edge at most `c` times (O(E))
+std::vector<NodeID> DinicMaxFlow::GetAugmentingPath(LevelGraph &levels,
+                                                    const NodeID node_id,
+                                                    const BisectionGraphView &view,
+                                                    const FlowEdges &flow,
+                                                    const SourceSinkNodes &source_nodes) const
+{
+    std::vector<NodeID> path;
+    BOOST_ASSERT(source_nodes.find(node_id) == source_nodes.end());
+
+    // Keeps the local state of the DFS in forms of the iterators
+    struct DFSState
+    {
+        BisectionGraph::ConstEdgeIterator edge_iterator;
+        const BisectionGraph::ConstEdgeIterator end_iterator;
+    };
+
+    std::stack<DFSState> dfs_stack;
+    DFSState initial_state = {view.BeginEdges(node_id), view.EndEdges(node_id)};
+    dfs_stack.push(std::move(initial_state));
+    path.push_back(node_id);
+
+    while (!dfs_stack.empty())
+    {
+        // the dfs_stack and the path have to be kept in sync
+        BOOST_ASSERT(dfs_stack.size() == path.size());
+
+        while (dfs_stack.top().edge_iterator != dfs_stack.top().end_iterator)
+        {
+            const auto target = dfs_stack.top().edge_iterator->target;
+
+            // look at every edge only once, so advance the state of the current node (last in
+            // path)
+            dfs_stack.top().edge_iterator++;
+
+            // check if the edge is valid
+            const auto has_capacity = flow[target].count(path.back()) == 0;
+            const auto descends_level_graph = levels[target] + 1 == levels[path.back()];
+
+            if (has_capacity && descends_level_graph)
+            {
+                // recurse
+                path.push_back(target);
+
+                // termination
+                if (source_nodes.find(target) != source_nodes.end())
+                {
+                    std::reverse(path.begin(), path.end());
+                    return path;
+                }
+
+                // start next iteration
+                dfs_stack.push({view.BeginEdges(target), view.EndEdges(target)});
+            }
+        }
+
+        // backtrack - mark that there is no way to the target
+        levels[path.back()] = -1;
+        path.pop_back();
+        dfs_stack.pop();
+    }
+    BOOST_ASSERT(path.empty());
+    return path;
+}
+
+bool DinicMaxFlow::Validate(const BisectionGraphView &view,
+                            const SourceSinkNodes &source_nodes,
+                            const SourceSinkNodes &sink_nodes) const
+{
+    // sink and source cannot share a common node
+    const auto separated =
+        std::find_if(source_nodes.begin(), source_nodes.end(), [&sink_nodes](const auto node) {
+            return sink_nodes.count(node);
+        }) == source_nodes.end();
+
+    const auto invalid_id = [&view](const NodeID nid) { return nid >= view.NumberOfNodes(); };
+    const auto in_range_source =
+        std::find_if(source_nodes.begin(), source_nodes.end(), invalid_id) == source_nodes.end();
+    const auto in_range_sink =
+        std::find_if(sink_nodes.begin(), sink_nodes.end(), invalid_id) == sink_nodes.end();
+
+    return separated && in_range_source && in_range_sink;
+}
+
+} // namespace partitioner
+} // namespace osrm

src/partitioner/inertial_flow.cpp

@@ -0,0 +1,163 @@
+#include "partitioner/inertial_flow.hpp"
+#include "partitioner/bisection_graph.hpp"
+#include "partitioner/bisection_graph_view.hpp"
+#include "partitioner/reorder_first_last.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <iterator>
+#include <mutex>
+#include <set>
+#include <tuple>
+#include <unordered_set>
+#include <utility>
+#include <vector>
+
+#include <tbb/blocked_range.h>
+#include <tbb/parallel_for.h>
+
+namespace osrm
+{
+namespace partitioner
+{
+namespace
+{
+// 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(const BisectionGraphView &view, const double ratio, const double slope)
+{
+    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());
+
+    // adress of the very first node
+    const auto node_zero = &(*view.Begin());
+
+    std::transform(view.Begin(), view.End(), std::back_inserter(embedding), [&](const auto &node) {
+        const auto node_id = static_cast<NodeID>(&node - node_zero);
+        return NodeWithCoordinate{node_id, node.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);
+
+    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;
+}
+
+// Makes n cuts with different spatial orders and returns the best.
+DinicMaxFlow::MinCut bestMinCut(const BisectionGraphView &view,
+                                const std::size_t n,
+                                const double ratio,
+                                const double balance)
+{
+    DinicMaxFlow::MinCut best;
+    best.num_edges = -1;
+
+    const auto get_balance = [&view, balance](const auto num_nodes_source) {
+        const auto perfect_balance = view.NumberOfNodes() / 2;
+        const auto allowed_balance = balance * perfect_balance;
+        const auto bigger_side =
+            std::max(num_nodes_source, view.NumberOfNodes() - num_nodes_source);
+
+        if (bigger_side > allowed_balance)
+            return bigger_side / static_cast<double>(allowed_balance);
+        else
+            return 1.0;
+    };
+
+    auto best_balance = 1;
+
+    std::mutex lock;
+
+    tbb::blocked_range<std::size_t> range{0, n, 1};
+
+    const auto balance_delta = [&view](const auto num_nodes_source) {
+        const std::int64_t difference =
+            static_cast<std::int64_t>(view.NumberOfNodes()) / 2 - num_nodes_source;
+        return std::abs(difference);
+    };
+
+    tbb::parallel_for(range, [&](const auto &chunk) {
+        for (auto round = chunk.begin(), end = chunk.end(); round != end; ++round)
+        {
+            const auto slope = -1. + round * (2. / n);
+
+            auto order = makeSpatialOrder(view, ratio, slope);
+            auto cut = DinicMaxFlow()(view, order.sources, order.sinks);
+            auto cut_balance = get_balance(cut.num_nodes_source);
+
+            {
+                std::lock_guard<std::mutex> guard{lock};
+
+                // Swap to keep the destruction of the old object outside of critical section.
+                if (cut.num_edges * cut_balance < best.num_edges * best_balance ||
+                    (cut.num_edges == best.num_edges &&
+                     balance_delta(cut.num_nodes_source) < balance_delta(best.num_nodes_source)))
+                {
+                    best_balance = cut_balance;
+                    std::swap(best, cut);
+                }
+            }
+            // cut gets destroyed here
+        }
+    });
+
+    return best;
+}
+}
+
+DinicMaxFlow::MinCut computeInertialFlowCut(const BisectionGraphView &view,
+                                            const std::size_t num_slopes,
+                                            const double balance,
+                                            const double source_sink_rate)
+{
+    return bestMinCut(view, num_slopes, source_sink_rate, balance);
+}
+
+} // namespace partitioner
+} // namespace osrm

src/partitioner/partitioner.cpp

@@ -0,0 +1,234 @@
+#include "partitioner/partitioner.hpp"
+#include "partitioner/bisection_graph.hpp"
+#include "partitioner/bisection_to_partition.hpp"
+#include "partitioner/cell_storage.hpp"
+#include "partitioner/compressed_node_based_graph_reader.hpp"
+#include "partitioner/edge_based_graph_reader.hpp"
+#include "partitioner/files.hpp"
+#include "partitioner/multi_level_partition.hpp"
+#include "partitioner/recursive_bisection.hpp"
+#include "partitioner/remove_unconnected.hpp"
+#include "partitioner/renumber.hpp"
+
+#include "extractor/files.hpp"
+
+#include "util/coordinate.hpp"
+#include "util/geojson_debug_logger.hpp"
+#include "util/geojson_debug_policies.hpp"
+#include "util/integer_range.hpp"
+#include "util/json_container.hpp"
+#include "util/log.hpp"
+#include "util/mmap_file.hpp"
+
+#include <algorithm>
+#include <iterator>
+#include <vector>
+
+#include <boost/assert.hpp>
+#include <boost/filesystem/operations.hpp>
+
+#include "util/geojson_debug_logger.hpp"
+#include "util/geojson_debug_policies.hpp"
+#include "util/json_container.hpp"
+#include "util/timing_util.hpp"
+
+namespace osrm
+{
+namespace partitioner
+{
+
+void LogGeojson(const std::string &filename, const std::vector<std::uint32_t> &bisection_ids)
+{
+    // reload graph, since we destroyed the old one
+    auto compressed_node_based_graph = LoadCompressedNodeBasedGraph(filename);
+
+    util::Log() << "Loaded compressed node based graph: "
+                << compressed_node_based_graph.edges.size() << " edges, "
+                << compressed_node_based_graph.coordinates.size() << " nodes";
+
+    groupEdgesBySource(begin(compressed_node_based_graph.edges),
+                       end(compressed_node_based_graph.edges));
+
+    auto graph =
+        makeBisectionGraph(compressed_node_based_graph.coordinates,
+                           adaptToBisectionEdge(std::move(compressed_node_based_graph.edges)));
+
+    const auto get_level = [](const std::uint32_t lhs, const std::uint32_t rhs) {
+        auto xored = lhs ^ rhs;
+        std::uint32_t level = log(xored) / log(2.0);
+        return level;
+    };
+
+    std::vector<std::vector<util::Coordinate>> border_vertices(33);
+
+    for (NodeID nid = 0; nid < graph.NumberOfNodes(); ++nid)
+    {
+        const auto source_id = bisection_ids[nid];
+        for (const auto &edge : graph.Edges(nid))
+        {
+            const auto target_id = bisection_ids[edge.target];
+            if (source_id != target_id)
+            {
+                auto level = get_level(source_id, target_id);
+                border_vertices[level].push_back(graph.Node(nid).coordinate);
+                border_vertices[level].push_back(graph.Node(edge.target).coordinate);
+            }
+        }
+    }
+
+    util::ScopedGeojsonLoggerGuard<util::CoordinateVectorToMultiPoint> guard(
+        "border_vertices.geojson");
+    std::size_t level = 0;
+    for (auto &bv : border_vertices)
+    {
+        if (!bv.empty())
+        {
+            std::sort(bv.begin(), bv.end(), [](const auto lhs, const auto rhs) {
+                return std::tie(lhs.lon, lhs.lat) < std::tie(rhs.lon, rhs.lat);
+            });
+            bv.erase(std::unique(bv.begin(), bv.end()), bv.end());
+
+            util::json::Object jslevel;
+            jslevel.values["level"] = util::json::Number(level++);
+            guard.Write(bv, jslevel);
+        }
+    }
+}
+
+auto getGraphBisection(const PartitionerConfig &config)
+{
+    auto compressed_node_based_graph =
+        LoadCompressedNodeBasedGraph(config.GetPath(".osrm.cnbg").string());
+
+    util::Log() << "Loaded compressed node based graph: "
+                << compressed_node_based_graph.edges.size() << " edges, "
+                << compressed_node_based_graph.coordinates.size() << " nodes";
+
+    groupEdgesBySource(begin(compressed_node_based_graph.edges),
+                       end(compressed_node_based_graph.edges));
+
+    auto graph =
+        makeBisectionGraph(compressed_node_based_graph.coordinates,
+                           adaptToBisectionEdge(std::move(compressed_node_based_graph.edges)));
+
+    util::Log() << " running partition: " << config.max_cell_sizes.front() << " " << config.balance
+                << " " << config.boundary_factor << " " << config.num_optimizing_cuts << " "
+                << config.small_component_size
+                << " # max_cell_size balance boundary cuts small_component_size";
+    RecursiveBisection recursive_bisection(graph,
+                                           config.max_cell_sizes.front(),
+                                           config.balance,
+                                           config.boundary_factor,
+                                           config.num_optimizing_cuts,
+                                           config.small_component_size);
+
+    // Return bisection ids, keyed by node based graph nodes
+    return recursive_bisection.BisectionIDs();
+}
+
+int Partitioner::Run(const PartitionerConfig &config)
+{
+    const std::vector<BisectionID> &node_based_partition_ids = getGraphBisection(config);
+
+    // Up until now we worked on the compressed node based graph.
+    // But what we actually need is a partition for the edge based graph to work on.
+    // The following loads a mapping from node based graph to edge based graph.
+    // Then loads the edge based graph tanslates the partition and modifies it.
+    // For details see #3205
+
+    std::vector<extractor::NBGToEBG> mapping;
+    extractor::files::readNBGMapping(config.GetPath(".osrm.cnbg_to_ebg").string(), mapping);
+    util::Log() << "Loaded node based graph to edge based graph mapping";
+
+    auto edge_based_graph = LoadEdgeBasedGraph(config.GetPath(".osrm.ebg").string());
+    util::Log() << "Loaded edge based graph for mapping partition ids: "
+                << edge_based_graph.GetNumberOfEdges() << " edges, "
+                << edge_based_graph.GetNumberOfNodes() << " nodes";
+
+    // Partition ids, keyed by edge based graph nodes
+    std::vector<NodeID> edge_based_partition_ids(edge_based_graph.GetNumberOfNodes(),
+                                                 SPECIAL_NODEID);
+
+    // Only resolve all easy cases in the first pass
+    for (const auto &entry : mapping)
+    {
+        const auto u = entry.u;
+        const auto v = entry.v;
+        const auto forward_node = entry.forward_ebg_node;
+        const auto backward_node = entry.backward_ebg_node;
+
+        // This heuristic strategy seems to work best, even beating chosing the minimum
+        // border edge bisection ID
+        edge_based_partition_ids[forward_node] = node_based_partition_ids[u];
+        if (backward_node != SPECIAL_NODEID)
+            edge_based_partition_ids[backward_node] = node_based_partition_ids[v];
+    }
+
+    std::vector<Partition> partitions;
+    std::vector<std::uint32_t> level_to_num_cells;
+    std::tie(partitions, level_to_num_cells) =
+        bisectionToPartition(edge_based_partition_ids, config.max_cell_sizes);
+
+    auto num_unconnected = removeUnconnectedBoundaryNodes(edge_based_graph, partitions);
+    util::Log() << "Fixed " << num_unconnected << " unconnected nodes";
+
+    util::Log() << "Edge-based-graph annotation:";
+    for (std::size_t level = 0; level < level_to_num_cells.size(); ++level)
+    {
+        util::Log() << "  level " << level + 1 << " #cells " << level_to_num_cells[level]
+                    << " bit size " << std::ceil(std::log2(level_to_num_cells[level] + 1));
+    }
+
+    TIMER_START(renumber);
+    auto permutation = makePermutation(edge_based_graph, partitions);
+    renumber(edge_based_graph, permutation);
+    renumber(partitions, permutation);
+    {
+        renumber(mapping, permutation);
+        extractor::files::writeNBGMapping(config.GetPath(".osrm.cnbg_to_ebg").string(), mapping);
+    }
+    {
+        boost::iostreams::mapped_file segment_region;
+        auto segments = util::mmapFile<extractor::EdgeBasedNodeSegment>(
+            config.GetPath(".osrm.fileIndex"), segment_region);
+        renumber(segments, permutation);
+    }
+    {
+        extractor::EdgeBasedNodeDataContainer node_data;
+        extractor::files::readNodeData(config.GetPath(".osrm.ebg_nodes"), node_data);
+        renumber(node_data, permutation);
+        extractor::files::writeNodeData(config.GetPath(".osrm.ebg_nodes"), node_data);
+    }
+    if (boost::filesystem::exists(config.GetPath(".osrm.hsgr")))
+    {
+        util::Log(logWARNING) << "Found existing .osrm.hsgr file, removing. You need to re-run "
+                                 "osrm-contract after osrm-partition.";
+        boost::filesystem::remove(config.GetPath(".osrm.hsgr"));
+    }
+    TIMER_STOP(renumber);
+    util::Log() << "Renumbered data in " << TIMER_SEC(renumber) << " seconds";
+
+    TIMER_START(packed_mlp);
+    MultiLevelPartition mlp{partitions, level_to_num_cells};
+    TIMER_STOP(packed_mlp);
+    util::Log() << "MultiLevelPartition constructed in " << TIMER_SEC(packed_mlp) << " seconds";
+
+    TIMER_START(cell_storage);
+    CellStorage storage(mlp, edge_based_graph);
+    TIMER_STOP(cell_storage);
+    util::Log() << "CellStorage constructed in " << TIMER_SEC(cell_storage) << " seconds";
+
+    TIMER_START(writing_mld_data);
+    files::writePartition(config.GetPath(".osrm.partition"), mlp);
+    files::writeCells(config.GetPath(".osrm.cells"), storage);
+    extractor::files::writeEdgeBasedGraph(config.GetPath(".osrm.ebg"),
+                                          edge_based_graph.GetNumberOfNodes(),
+                                          graphToEdges(edge_based_graph));
+    TIMER_STOP(writing_mld_data);
+    util::Log() << "MLD data writing took " << TIMER_SEC(writing_mld_data) << " seconds";
+
+    return 0;
+}
+
+} // namespace partitioner
+} // namespace osrm

src/partitioner/recursive_bisection.cpp

@@ -0,0 +1,111 @@
+#include "partitioner/recursive_bisection.hpp"
+#include "partitioner/inertial_flow.hpp"
+
+#include "partitioner/bisection_graph_view.hpp"
+#include "partitioner/recursive_bisection_state.hpp"
+
+#include "util/log.hpp"
+#include "util/timing_util.hpp"
+
+#include <tbb/parallel_do.h>
+
+#include <algorithm>
+#include <climits> // for CHAR_BIT
+#include <cstddef>
+#include <iterator>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+namespace osrm
+{
+namespace partitioner
+{
+
+RecursiveBisection::RecursiveBisection(BisectionGraph &bisection_graph_,
+                                       const std::size_t maximum_cell_size,
+                                       const double balance,
+                                       const double boundary_factor,
+                                       const std::size_t num_optimizing_cuts,
+                                       const std::size_t small_component_size)
+    : bisection_graph(bisection_graph_), internal_state(bisection_graph_)
+{
+    auto components = internal_state.PrePartitionWithSCC(small_component_size);
+    BOOST_ASSERT(!components.empty());
+
+    // Parallelize recursive bisection trees. Root cut happens serially (well, this is a lie:
+    // since we handle big components in parallel, too. But we don't know this and
+    // don't have to. TBB's scheduler handles nested parallelism just fine).
+    //
+    //     [   |   ]
+    //      /     \         root cut
+    //  [ | ]     [ | ]
+    //  /   \     /   \     descend, do cuts in parallel
+    //
+    // https://www.threadingbuildingblocks.org/docs/help/index.htm#reference/algorithms/parallel_do_func.html
+
+    struct TreeNode
+    {
+        BisectionGraphView graph;
+        std::uint64_t depth;
+    };
+
+    // Build a recursive bisection tree for all big components independently in parallel.
+    // Last GraphView is all small components: skip for bisection.
+    auto first = begin(components);
+    auto last = end(components) - 1;
+
+    // We construct the trees on the fly: the root node is the entry point.
+    // All tree branches depend on the actual cut and will be generated while descending.
+    std::vector<TreeNode> forest;
+    forest.reserve(last - first);
+
+    std::transform(first, last, std::back_inserter(forest), [this](auto graph) {
+        return TreeNode{std::move(graph), internal_state.SCCDepth()};
+    });
+
+    using Feeder = tbb::parallel_do_feeder<TreeNode>;
+
+    TIMER_START(bisection);
+
+    // Bisect graph into two parts. Get partition point and recurse left and right in parallel.
+    tbb::parallel_do(begin(forest), end(forest), [&](const TreeNode &node, Feeder &feeder) {
+        const auto cut =
+            computeInertialFlowCut(node.graph, num_optimizing_cuts, balance, boundary_factor);
+        const auto center = internal_state.ApplyBisection(
+            node.graph.Begin(), node.graph.End(), node.depth, cut.flags);
+
+        const auto terminal = [&](const auto &node) {
+            const auto maximum_depth = sizeof(BisectionID) * CHAR_BIT;
+            const auto too_small = node.graph.NumberOfNodes() < maximum_cell_size;
+            const auto too_deep = node.depth >= maximum_depth;
+            return too_small || too_deep;
+        };
+
+        BisectionGraphView left_graph{bisection_graph, node.graph.Begin(), center};
+        TreeNode left_node{std::move(left_graph), node.depth + 1};
+
+        if (!terminal(left_node))
+            feeder.add(std::move(left_node));
+
+        BisectionGraphView right_graph{bisection_graph, center, node.graph.End()};
+        TreeNode right_node{std::move(right_graph), node.depth + 1};
+
+        if (!terminal(right_node))
+            feeder.add(std::move(right_node));
+    });
+
+    TIMER_STOP(bisection);
+
+    util::Log() << "Full bisection done in " << TIMER_SEC(bisection) << "s";
+}
+
+const std::vector<BisectionID> &RecursiveBisection::BisectionIDs() const
+{
+    return internal_state.BisectionIDs();
+}
+
+std::uint32_t RecursiveBisection::SCCDepth() const { return internal_state.SCCDepth(); }
+
+} // namespace partitioner
+} // namespace osrm

src/partitioner/recursive_bisection_state.cpp

@@ -0,0 +1,185 @@
+#include "partitioner/recursive_bisection_state.hpp"
+#include "extractor/tarjan_scc.hpp"
+#include "partitioner/tarjan_graph_wrapper.hpp"
+
+#include <algorithm>
+#include <climits> // for CHAR_BIT
+#include <numeric>
+#include <set>
+#include <unordered_map>
+
+namespace osrm
+{
+namespace partitioner
+{
+
+RecursiveBisectionState::RecursiveBisectionState(BisectionGraph &bisection_graph_)
+    : scc_levels(0), bisection_graph(bisection_graph_)
+{
+    bisection_ids.resize(bisection_graph.NumberOfNodes(), BisectionID{0});
+}
+
+RecursiveBisectionState::~RecursiveBisectionState() {}
+
+BisectionID RecursiveBisectionState::GetBisectionID(const NodeID node) const
+{
+    return bisection_ids[node];
+}
+
+RecursiveBisectionState::NodeIterator
+RecursiveBisectionState::ApplyBisection(const NodeIterator const_begin,
+                                        const NodeIterator const_end,
+                                        const std::size_t depth,
+                                        const std::vector<bool> &partition)
+{
+    BOOST_ASSERT(depth >= scc_levels);
+    // ensure that the iterators belong to the graph
+    BOOST_ASSERT(bisection_graph.GetID(*const_begin) < bisection_graph.NumberOfNodes() &&
+                 bisection_graph.GetID(*const_begin) + std::distance(const_begin, const_end) <=
+                     bisection_graph.NumberOfNodes());
+    // augment the partition ids
+    const auto flag = BisectionID{1} << (sizeof(BisectionID) * CHAR_BIT - depth - 1);
+    for (auto itr = const_begin; itr != const_end; ++itr)
+    {
+        const auto nid = std::distance(const_begin, itr);
+        if (partition[nid])
+            bisection_ids[itr->original_id] |= flag;
+    }
+
+    // Keep items with `0` as partition id to the left, move other to the right
+    auto by_flag_bit = [this, flag](const auto &node) {
+        return BisectionID{0} == (bisection_ids[node.original_id] & flag);
+    };
+
+    auto begin = bisection_graph.Begin() + std::distance(bisection_graph.CBegin(), const_begin);
+    const auto end = begin + std::distance(const_begin, const_end);
+
+    // remap the edges
+    std::vector<NodeID> mapping(std::distance(const_begin, const_end), SPECIAL_NODEID);
+    // calculate a mapping of all node ids
+    std::size_t lesser_id = 0, upper_id = 0;
+    std::transform(const_begin,
+                   const_end,
+                   mapping.begin(),
+                   [by_flag_bit, &lesser_id, &upper_id](const auto &node) {
+                       return by_flag_bit(node) ? lesser_id++ : upper_id++;
+                   });
+
+    // erase all edges that point into different partitions
+    std::for_each(begin, end, [&](auto &node) {
+        const auto node_flag = by_flag_bit(node);
+        bisection_graph.RemoveEdges(node, [&](const BisectionGraph::EdgeT &edge) {
+            const auto target_flag = by_flag_bit(*(const_begin + edge.target));
+            return (node_flag != target_flag);
+        });
+    });
+
+    auto center = std::stable_partition(begin, end, by_flag_bit);
+
+    // remap all remaining edges
+    std::for_each(const_begin, const_end, [&](const auto &node) {
+        for (auto &edge : bisection_graph.Edges(node))
+            edge.target = mapping[edge.target];
+    });
+
+    return const_begin + std::distance(begin, center);
+}
+
+std::vector<BisectionGraphView>
+RecursiveBisectionState::PrePartitionWithSCC(const std::size_t small_component_size)
+{
+    // since our graphs are unidirectional, we don't realy need the scc. But tarjan is so nice and
+    // assigns IDs and counts sizes
+    TarjanGraphWrapper wrapped_graph(bisection_graph);
+    extractor::TarjanSCC<TarjanGraphWrapper> scc_algo(wrapped_graph);
+    scc_algo.Run();
+
+    // Map Edges to Sccs
+    const auto in_small = [&scc_algo, small_component_size](const NodeID node_id) {
+        return scc_algo.GetComponentSize(scc_algo.GetComponentID(node_id)) <= small_component_size;
+    };
+
+    const constexpr std::size_t small_component_id = -1;
+    std::unordered_map<std::size_t, std::size_t> component_map;
+    const auto transform_id = [&](const NodeID node_id) -> std::size_t {
+        if (in_small(node_id))
+            return small_component_id;
+        else
+            return scc_algo.GetComponentID(node_id);
+    };
+
+    std::vector<NodeID> mapping(bisection_graph.NumberOfNodes(), SPECIAL_NODEID);
+    for (const auto &node : bisection_graph.Nodes())
+        mapping[node.original_id] = component_map[transform_id(node.original_id)]++;
+
+    // needs to remove edges, if we should ever switch to directed graphs here
+    std::stable_sort(
+        bisection_graph.Begin(), bisection_graph.End(), [&](const auto &lhs, const auto &rhs) {
+            return transform_id(lhs.original_id) < transform_id(rhs.original_id);
+        });
+
+    // remap all remaining edges
+    std::for_each(bisection_graph.Begin(), bisection_graph.End(), [&](const auto &node) {
+        for (auto &edge : bisection_graph.Edges(node))
+            edge.target = mapping[edge.target];
+    });
+
+    std::vector<BisectionGraphView> views;
+    auto last = bisection_graph.CBegin();
+    auto last_id = transform_id(bisection_graph.Begin()->original_id);
+    std::set<std::size_t> ordered_component_ids;
+    for (auto itr = bisection_graph.CBegin(); itr != bisection_graph.CEnd(); ++itr)
+    {
+        auto itr_id = transform_id(itr->original_id);
+        ordered_component_ids.insert(itr_id);
+        if (last_id != itr_id)
+        {
+            views.push_back(BisectionGraphView(bisection_graph, last, itr));
+            last_id = itr_id;
+            last = itr;
+        }
+    }
+    views.push_back(BisectionGraphView(bisection_graph, last, bisection_graph.CEnd()));
+
+    bool has_small_component = [&]() {
+        for (std::size_t i = 0; i < scc_algo.GetNumberOfComponents(); ++i)
+            if (scc_algo.GetComponentSize(i) <= small_component_size)
+                return true;
+        return false;
+    }();
+
+    if (!has_small_component)
+        views.push_back(
+            BisectionGraphView(bisection_graph, bisection_graph.CEnd(), bisection_graph.CEnd()));
+
+    // apply scc as bisections, we need scc_level bits for this with scc_levels =
+    // ceil(log_2(components))
+    scc_levels = ceil(log(views.size()) / log(2.0));
+
+    const auto conscutive_component_id = [&](const NodeID nid) {
+        const auto component_id = transform_id(nid);
+        const auto itr = ordered_component_ids.find(component_id);
+        BOOST_ASSERT(itr != ordered_component_ids.end());
+        BOOST_ASSERT(static_cast<std::size_t>(std::distance(ordered_component_ids.begin(), itr)) <
+                     ordered_component_ids.size());
+        return std::distance(ordered_component_ids.begin(), itr);
+    };
+
+    const auto shift = sizeof(BisectionID) * CHAR_BIT - scc_levels;
+
+    // store the component ids as first part of the bisection id
+    for (const auto &node : bisection_graph.Nodes())
+        bisection_ids[node.original_id] = conscutive_component_id(node.original_id) << shift;
+
+    return views;
+}
+
+const std::vector<BisectionID> &RecursiveBisectionState::BisectionIDs() const
+{
+    return bisection_ids;
+}
+
+std::uint32_t RecursiveBisectionState::SCCDepth() const { return scc_levels; }
+
+} // namespace partitioner
+} // namespace osrm

src/partitioner/renumber.cpp

@@ -0,0 +1,70 @@
+#include "partitioner/renumber.hpp"
+
+#include "util/permutation.hpp"
+
+namespace osrm
+{
+namespace partitioner
+{
+namespace
+{
+// Returns a vector that is indexed by node ID marking the level at which it is a border node
+std::vector<LevelID> getHighestBorderLevel(const DynamicEdgeBasedGraph &graph,
+                                           const std::vector<Partition> &partitions)
+{
+    std::vector<LevelID> border_level(graph.GetNumberOfNodes(), 0);
+
+    for (const auto level : util::irange<LevelID>(1, partitions.size() + 1))
+    {
+        const auto &partition = partitions[level - 1];
+        for (auto node : util::irange<NodeID>(0, graph.GetNumberOfNodes()))
+        {
+            for (auto edge : graph.GetAdjacentEdgeRange(node))
+            {
+                auto target = graph.GetTarget(edge);
+                if (partition[node] != partition[target])
+                {
+                    // level is monotone increasing so we wil
+                    // always overwrite here with a value equal
+                    // or greater then the current border_level
+                    border_level[node] = level;
+                    border_level[target] = level;
+                }
+            }
+        }
+    }
+
+    return border_level;
+}
+}
+
+std::vector<std::uint32_t> makePermutation(const DynamicEdgeBasedGraph &graph,
+                                           const std::vector<Partition> &partitions)
+{
+    std::vector<std::uint32_t> ordering(graph.GetNumberOfNodes());
+    std::iota(ordering.begin(), ordering.end(), 0);
+
+    // Sort the nodes by cell ID recursively:
+    // Nodes in the same cell will be sorted by cell ID on the level below
+    for (const auto &partition : partitions)
+    {
+        std::stable_sort(
+            ordering.begin(), ordering.end(), [&partition](const auto lhs, const auto rhs) {
+                return partition[lhs] < partition[rhs];
+            });
+    }
+
+    // Now sort the nodes by the level at which they are a border node, descening.
+    // That means nodes that are border nodes on the highest level will have a very low ID,
+    // whereas nodes that are nerver border nodes are sorted to the end of the array.
+    // Note: Since we use a stable sort that preserves the cell sorting within each level
+    auto border_level = getHighestBorderLevel(graph, partitions);
+    std::stable_sort(
+        ordering.begin(), ordering.end(), [&border_level](const auto lhs, const auto rhs) {
+            return border_level[lhs] > border_level[rhs];
+        });
+
+    return util::orderingToPermutation(ordering);
+}
+}
+}

src/partitioner/tarjan_graph_wrapper.cpp

@@ -0,0 +1,26 @@
+#include "partitioner/tarjan_graph_wrapper.hpp"
+
+namespace osrm
+{
+namespace partitioner
+{
+
+TarjanGraphWrapper::TarjanGraphWrapper(const BisectionGraph &bisection_graph_)
+    : bisection_graph(bisection_graph_)
+{
+}
+
+std::size_t TarjanGraphWrapper::GetNumberOfNodes() const { return bisection_graph.NumberOfNodes(); }
+
+util::range<EdgeID> TarjanGraphWrapper::GetAdjacentEdgeRange(const NodeID nid) const
+{
+    return util::irange<EdgeID>(bisection_graph.BeginEdgeID(nid), bisection_graph.EndEdgeID(nid));
+}
+
+NodeID TarjanGraphWrapper::GetTarget(const EdgeID eid) const
+{
+    return bisection_graph.Edge(eid).target;
+}
+
+} // namespace partitioner
+} // namespace osrm