Manhattan distance MST
(Mylib/Graph/MinimumSpanningTree/manhattan_minimum_spanning_tree.cpp)

• View this file on GitHub
• Last update: 2021-04-23 23:44:44+09:00
• Link: View error logs on GitHub Actions

Operations

Requirements

Notes

Problems

References

• https://knuu.github.io/manhattan_mst.html

Depends on

• Min monoid (Mylib/AlgebraicStructure/Monoid/min.cpp)
• Segment tree (Mylib/DataStructure/SegmentTree/segment_tree.cpp)
• Basic graph (Mylib/Graph/Template/graph.cpp)

Verified with

• test/yosupo-judge/manhattanmst/main.test.cpp

Code

#pragma once
#include <algorithm>
#include <cmath>
#include <numeric>
#include <utility>
#include <vector>
#include "Mylib/AlgebraicStructure/Monoid/min.cpp"
#include "Mylib/DataStructure/SegmentTree/segment_tree.cpp"
#include "Mylib/Graph/Template/graph.cpp"

namespace haar_lib {
  template <typename T, typename MST>
  std::vector<edge<T>> manhattan_minimum_spanning_tree(std::vector<T> x, std::vector<T> y, MST mst) {
    const int N = x.size();
    graph<T> g(N);
    segment_tree<min_monoid<std::pair<T, int>>> seg(N);

    auto f =
        [&]() {
          std::vector<T> Y(y);
          std::sort(Y.begin(), Y.end());
          Y.erase(std::unique(Y.begin(), Y.end()), Y.end());

          seg.init(std::nullopt);

          std::vector<int> ord(N);
          std::iota(ord.begin(), ord.end(), 0);
          std::sort(
              ord.begin(), ord.end(),
              [&](int i, int j) {
                if (y[i] - x[i] == y[j] - x[j]) return x[i] > x[j];
                return y[i] - x[i] < y[j] - x[j];
              });

          for (int i : ord) {
            int lb = std::lower_bound(Y.begin(), Y.end(), y[i]) - Y.begin();

            if (auto res = seg.fold(lb, N); res) {
              auto j = res->second;
              T c    = std::abs(x[i] - x[j]) + std::abs(y[i] - y[j]);
              g.add_edge(i, j, c);
            }

            if (auto res = seg[lb]; not res or x[i] + y[i] < res->first) {
              seg.set(lb, {{x[i] + y[i], i}});
            }
          }
        };

    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        for (int k = 0; k < 2; ++k) {
          f();
          for (int l = 0; l < N; ++l) std::swap(x[l], y[l]);
        }
        for (int l = 0; l < N; ++l) x[l] = -x[l];
      }
      for (int l = 0; l < N; ++l) y[l] = -y[l];
    }

    return mst(g);
  }
}  // namespace haar_lib

#line 2 "Mylib/Graph/MinimumSpanningTree/manhattan_minimum_spanning_tree.cpp"
#include <algorithm>
#include <cmath>
#include <numeric>
#include <utility>
#include <vector>
#line 3 "Mylib/AlgebraicStructure/Monoid/min.cpp"
#include <optional>

namespace haar_lib {
  template <typename T>
  struct min_monoid {
    using value_type = std::optional<T>;

    value_type operator()() const { return {}; }
    value_type operator()(const value_type &a, const value_type &b) const {
      if (not a) return b;
      if (not b) return a;
      return {std::min(*a, *b)};
    }
  };
}  // namespace haar_lib
#line 3 "Mylib/DataStructure/SegmentTree/segment_tree.cpp"
#include <cassert>
#include <functional>
#line 6 "Mylib/DataStructure/SegmentTree/segment_tree.cpp"

namespace haar_lib {
  template <typename Monoid>
  class segment_tree {
  public:
    using value_type = typename Monoid::value_type;

  private:
    Monoid M_;
    int depth_, size_, hsize_;
    std::vector<value_type> data_;

  public:
    segment_tree() {}
    segment_tree(int n) : depth_(n > 1 ? 32 - __builtin_clz(n - 1) + 1 : 1),
                          size_(1 << depth_),
                          hsize_(size_ / 2),
                          data_(size_, M_()) {}

    auto operator[](int i) const {
      assert(0 <= i and i < hsize_);
      return data_[hsize_ + i];
    }

    auto fold(int l, int r) const {
      assert(0 <= l and l <= r and r <= hsize_);
      value_type ret_left  = M_();
      value_type ret_right = M_();

      int L = l + hsize_, R = r + hsize_;
      while (L < R) {
        if (R & 1) ret_right = M_(data_[--R], ret_right);
        if (L & 1) ret_left = M_(ret_left, data_[L++]);
        L >>= 1, R >>= 1;
      }

      return M_(ret_left, ret_right);
    }

    auto fold_all() const {
      return data_[1];
    }

    void set(int i, const value_type &x) {
      assert(0 <= i and i < hsize_);
      i += hsize_;
      data_[i] = x;
      while (i > 1) i >>= 1, data_[i] = M_(data_[i << 1 | 0], data_[i << 1 | 1]);
    }

    void update(int i, const value_type &x) {
      assert(0 <= i and i < hsize_);
      i += hsize_;
      data_[i] = M_(data_[i], x);
      while (i > 1) i >>= 1, data_[i] = M_(data_[i << 1 | 0], data_[i << 1 | 1]);
    }

    template <typename T>
    void init_with_vector(const std::vector<T> &val) {
      data_.assign(size_, M_());
      for (int i = 0; i < (int) val.size(); ++i) data_[hsize_ + i] = val[i];
      for (int i = hsize_; --i >= 1;) data_[i] = M_(data_[i << 1 | 0], data_[i << 1 | 1]);
    }

    template <typename T>
    void init(const T &val) {
      init_with_vector(std::vector<value_type>(hsize_, val));
    }

  private:
    template <bool Lower, typename F>
    int bound(const int l, const int r, value_type x, F f) const {
      std::vector<int> pl, pr;
      int L = l + hsize_;
      int R = r + hsize_;
      while (L < R) {
        if (R & 1) pr.push_back(--R);
        if (L & 1) pl.push_back(L++);
        L >>= 1, R >>= 1;
      }

      std::reverse(pr.begin(), pr.end());
      pl.insert(pl.end(), pr.begin(), pr.end());

      value_type a = M_();

      for (int i : pl) {
        auto b = M_(a, data_[i]);

        if ((Lower and not f(b, x)) or (not Lower and f(x, b))) {
          while (i < hsize_) {
            const auto c = M_(a, data_[i << 1 | 0]);
            if ((Lower and not f(c, x)) or (not Lower and f(x, c))) {
              i = i << 1 | 0;
            } else {
              a = c;
              i = i << 1 | 1;
            }
          }

          return i - hsize_;
        }

        a = b;
      }

      return r;
    }

  public:
    template <typename F = std::less<value_type>>
    int lower_bound(int l, int r, value_type x, F f = F()) const {
      return bound<true>(l, r, x, f);
    }

    template <typename F = std::less<value_type>>
    int upper_bound(int l, int r, value_type x, F f = F()) const {
      return bound<false>(l, r, x, f);
    }
  };
}  // namespace haar_lib
#line 2 "Mylib/Graph/Template/graph.cpp"
#include <iostream>
#line 4 "Mylib/Graph/Template/graph.cpp"

namespace haar_lib {
  template <typename T>
  struct edge {
    int from, to;
    T cost;
    int index = -1;
    edge() {}
    edge(int from, int to, T cost) : from(from), to(to), cost(cost) {}
    edge(int from, int to, T cost, int index) : from(from), to(to), cost(cost), index(index) {}
  };

  template <typename T>
  struct graph {
    using weight_type = T;
    using edge_type   = edge<T>;

    std::vector<std::vector<edge<T>>> data;

    auto& operator[](size_t i) { return data[i]; }
    const auto& operator[](size_t i) const { return data[i]; }

    auto begin() const { return data.begin(); }
    auto end() const { return data.end(); }

    graph() {}
    graph(int N) : data(N) {}

    bool empty() const { return data.empty(); }
    int size() const { return data.size(); }

    void add_edge(int i, int j, T w, int index = -1) {
      data[i].emplace_back(i, j, w, index);
    }

    void add_undirected(int i, int j, T w, int index = -1) {
      add_edge(i, j, w, index);
      add_edge(j, i, w, index);
    }

    template <size_t I, bool DIRECTED = true, bool WEIGHTED = true>
    void read(int M) {
      for (int i = 0; i < M; ++i) {
        int u, v;
        std::cin >> u >> v;
        u -= I;
        v -= I;
        T w = 1;
        if (WEIGHTED) std::cin >> w;
        if (DIRECTED)
          add_edge(u, v, w, i);
        else
          add_undirected(u, v, w, i);
      }
    }
  };

  template <typename T>
  using tree = graph<T>;
}  // namespace haar_lib
#line 10 "Mylib/Graph/MinimumSpanningTree/manhattan_minimum_spanning_tree.cpp"

namespace haar_lib {
  template <typename T, typename MST>
  std::vector<edge<T>> manhattan_minimum_spanning_tree(std::vector<T> x, std::vector<T> y, MST mst) {
    const int N = x.size();
    graph<T> g(N);
    segment_tree<min_monoid<std::pair<T, int>>> seg(N);

    auto f =
        [&]() {
          std::vector<T> Y(y);
          std::sort(Y.begin(), Y.end());
          Y.erase(std::unique(Y.begin(), Y.end()), Y.end());

          seg.init(std::nullopt);

          std::vector<int> ord(N);
          std::iota(ord.begin(), ord.end(), 0);
          std::sort(
              ord.begin(), ord.end(),
              [&](int i, int j) {
                if (y[i] - x[i] == y[j] - x[j]) return x[i] > x[j];
                return y[i] - x[i] < y[j] - x[j];
              });

          for (int i : ord) {
            int lb = std::lower_bound(Y.begin(), Y.end(), y[i]) - Y.begin();

            if (auto res = seg.fold(lb, N); res) {
              auto j = res->second;
              T c    = std::abs(x[i] - x[j]) + std::abs(y[i] - y[j]);
              g.add_edge(i, j, c);
            }

            if (auto res = seg[lb]; not res or x[i] + y[i] < res->first) {
              seg.set(lb, {{x[i] + y[i], i}});
            }
          }
        };

    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        for (int k = 0; k < 2; ++k) {
          f();
          for (int l = 0; l < N; ++l) std::swap(x[l], y[l]);
        }
        for (int l = 0; l < N; ++l) x[l] = -x[l];
      }
      for (int l = 0; l < N; ++l) y[l] = -y[l];
    }

    return mst(g);
  }
}  // namespace haar_lib

Back to top page