TSP (2-opt, 3-opt, double-bridge, Held-Karp lower bound) (#274)

* TSP

* TSP: add doc

Author: hitonanode

tsp/calc_lkh_alpha.hpp

@@ -0,0 +1,61 @@
+#pragma once
+
+#include <vector>
+
+// http://webhotel4.ruc.dk/~keld/research/LKH/LKH-2.0/DOC/LKH_REPORT.pdf, p.20
+template <class DistanceMatrix> auto calc_lkh_alpha(const DistanceMatrix &dist) {
+    using T = decltype((*dist.adjacents(0).begin()).second);
+
+    std::vector<std::vector<int>> to(dist.n());
+
+    for (auto [s, t] : mst_edges(dist)) {
+        to.at(s).push_back(t);
+        to.at(t).push_back(s);
+    }
+
+    std::vector ret(dist.n(), std::vector<T>(dist.n()));
+
+    for (int s = 0; s < dist.n(); ++s) {
+        auto rec = [&](auto &&self, int now, int prv, T hi) -> void {
+            ret.at(s).at(now) = dist(s, now) - hi;
+            for (int nxt : to.at(now)) {
+                if (nxt == prv) continue;
+                self(self, nxt, now, std::max(hi, dist(now, nxt)));
+            }
+        };
+        rec(rec, s, -1, T());
+    }
+
+    // Determining special node for the 1-tree
+    // Reference: p.26 of http://webhotel4.ruc.dk/~keld/research/LKH/LKH-2.0/DOC/LKH_REPORT.pdf
+    int best_one = -1;
+    T longest_2nd_nearest = T();
+
+    for (int one = 0; one < dist.n(); ++one) {
+        if (to.at(one).size() != 1) continue;
+        const int ng = to.at(one).front();
+        bool found = false;
+        T second_nearest = T();
+
+        for (auto [v, len] : dist.adjacents(one)) {
+            if (v == ng) continue;
+            if (!found) {
+                found = true, second_nearest = len;
+            } else if (len < second_nearest) {
+                second_nearest = len;
+            }
+        }
+
+        if (found and (best_one < 0 or second_nearest > longest_2nd_nearest))
+            best_one = one, longest_2nd_nearest = second_nearest;
+    }
+
+    if (best_one != -1) {
+        for (auto [v, len] : dist.adjacents(best_one)) {
+            if (v == to.at(best_one).front()) continue;
+            ret.at(best_one).at(v) = ret.at(v).at(best_one) = len - longest_2nd_nearest;
+        }
+    }
+
+    return ret;
+}

tsp/csr_distance_matrix.hpp

@@ -0,0 +1,86 @@
+#pragma once
+
+#include <algorithm>
+#include <utility>
+#include <vector>
+
+template <class T> class csr_distance_matrix {
+
+    int _rows = 0;
+    std::vector<int> begins;
+    std::vector<std::pair<int, T>> vals;
+
+public:
+    csr_distance_matrix() : csr_distance_matrix({}, 0) {}
+
+    csr_distance_matrix(const std::vector<std::tuple<int, int, T>> &init, int rows)
+        : _rows(rows), begins(rows + 1) {
+        std::vector<int> degs(rows);
+        for (const auto &p : init) ++degs.at(std::get<0>(p));
+
+        for (int i = 0; i < rows; ++i) begins.at(i + 1) = begins.at(i) + degs.at(i);
+
+        vals.resize(init.size(), std::make_pair(-1, T()));
+        for (auto [i, j, w] : init) vals.at(begins.at(i + 1) - (degs.at(i)--)) = {j, w};
+    }
+
+    void apply_pi(const std::vector<T> &pi) {
+        for (int i = 0; i < n(); ++i) {
+            for (auto &[j, d] : adjacents(i)) d += pi.at(i) + pi.at(j);
+        }
+    }
+
+    int n() const noexcept { return _rows; }
+
+    struct adjacents_sequence {
+        csr_distance_matrix *ptr;
+        int from;
+
+        using iterator = typename std::vector<std::pair<int, T>>::iterator;
+        iterator begin() { return std::next(ptr->vals.begin(), ptr->begins.at(from)); }
+        iterator end() { return std::next(ptr->vals.begin(), ptr->begins.at(from + 1)); }
+    };
+
+    struct const_adjacents_sequence {
+        const csr_distance_matrix *ptr;
+        const int from;
+
+        using const_iterator = typename std::vector<std::pair<int, T>>::const_iterator;
+        const_iterator begin() const {
+            return std::next(ptr->vals.cbegin(), ptr->begins.at(from));
+        }
+        const_iterator end() const {
+            return std::next(ptr->vals.cbegin(), ptr->begins.at(from + 1));
+        }
+    };
+
+    adjacents_sequence adjacents(int from) { return {this, from}; }
+
+    const_adjacents_sequence adjacents(int from) const { return {this, from}; }
+    const_adjacents_sequence operator()(int from) const { return {this, from}; }
+};
+
+template <class DistanceMatrix> auto build_adjacent_info(const DistanceMatrix &dist, int sz) {
+    using T = decltype((*dist.adjacents(0).begin()).second);
+
+    const std::vector<std::vector<T>> alpha = calc_lkh_alpha(dist);
+
+    std::vector<std::tuple<int, int, T>> adjacent_edges;
+
+    std::vector<std::tuple<T, T, int>> candidates;
+    for (int i = 0; i < dist.n(); ++i) {
+        candidates.clear();
+        for (auto [j, d] : dist.adjacents(i)) {
+            if (i != j) candidates.emplace_back(alpha.at(i).at(j), d, j);
+        }
+
+        const int final_sz = std::min<int>(sz, candidates.size());
+        std::nth_element(candidates.begin(), candidates.begin() + final_sz, candidates.end());
+
+        candidates.resize(final_sz);
+        std::sort(candidates.begin(), candidates.end(),
+                  [&](const auto &l, const auto &r) { return std::get<1>(l) < std::get<1>(r); });
+        for (auto [alpha, dij, j] : candidates) adjacent_edges.emplace_back(i, j, dij);
+    }
+    return csr_distance_matrix(adjacent_edges, dist.n());
+}

tsp/dense_distance_matrix.hpp

@@ -0,0 +1,45 @@
+#pragma once
+
+#include <vector>
+
+template <class T> class dense_distance_matrix {
+    int _n;
+    std::vector<T> _d;
+
+public:
+    dense_distance_matrix(const std::vector<std::vector<T>> &distance_vecvec)
+        : _n(distance_vecvec.size()) {
+        _d.reserve(n() * n());
+        for (const auto &vec : distance_vecvec) _d.insert(end(_d), begin(vec), end(vec));
+    }
+
+    template <class U> void apply_pi(const std::vector<U> &pi) {
+        for (int i = 0; i < n(); ++i) {
+            for (int j = 0; j < n(); ++j) _d.at(i * n() + j) += pi.at(i) + pi.at(j);
+        }
+    }
+
+    int n() const noexcept { return _n; }
+
+    T dist(int s, int t) const { return _d.at(s * n() + t); }
+    T operator()(int s, int t) const { return dist(s, t); }
+
+    struct adjacents_sequence {
+        const dense_distance_matrix *ptr;
+        int from;
+        struct iterator {
+            const dense_distance_matrix *ptr;
+            int from;
+            int to;
+            iterator operator++() { return {ptr, from, to++}; }
+            std::pair<int, T> operator*() const { return {to, ptr->dist(from, to)}; }
+            bool operator!=(const iterator &rhs) const {
+                return to != rhs.to or ptr != rhs.ptr or from != rhs.from;
+            }
+        };
+        iterator begin() const { return iterator{ptr, from, 0}; }
+        iterator end() const { return iterator{ptr, from, ptr->n()}; }
+    };
+
+    adjacents_sequence adjacents(int from) const { return {this, from}; }
+};

tsp/held_karp.hpp

@@ -0,0 +1,55 @@
+#pragma once
+
+#include "csr_distance_matrix.hpp"
+#include "minimum_one_tree.hpp"
+
+#include <algorithm>
+#include <utility>
+#include <vector>
+
+// Held–Karp lower bound
+// http://webhotel4.ruc.dk/~keld/research/LKH/LKH-2.0/DOC/LKH_REPORT.pdf
+// p.26, p.33
+template <class DistanceMatrix> auto held_karp_lower_bound(const DistanceMatrix &dist) {
+    using T = decltype((*dist.adjacents(0).begin()).second);
+
+    std::vector<T> best_pi(dist.n()), pi(dist.n());
+    std::vector<int> V;
+    T W;
+    std::tie(W, V) = minimum_one_tree(dist, pi);
+    if (std::count(V.cbegin(), V.cend(), 0) == dist.n()) return std::make_pair(W, pi);
+
+    std::vector<int> lastV = V;
+
+    T BestW = W;
+    const int initial_period = (dist.n() + 1) / 2;
+    bool is_initial_phase = true;
+    int period = initial_period;
+
+    const auto sparse_subgraph = build_adjacent_info(dist, 50); // p.47
+
+    for (long long t0 = 1; t0 > 0; period /= 2, t0 /= 2) {
+        for (int p = 1; t0 > 0 and p <= period; ++p) {
+            for (int i = 0; i < dist.n(); ++i) {
+                if (V.at(i) != 0) pi.at(i) += t0 * (7 * V.at(i) + 3 * lastV.at(i)) / 10;
+            }
+            std::swap(lastV, V);
+            std::tie(W, V) = minimum_one_tree(sparse_subgraph, pi);
+            if (std::count(begin(V), begin(V) + dist.n(), 0) == dist.n())
+                return std::make_pair(W, pi);
+            if (W > BestW) {
+                BestW = W;
+                best_pi = pi;
+                if (is_initial_phase) t0 *= 2;
+                if (p == period) period = std::min(period * 2, initial_period);
+
+            } else if (is_initial_phase and p > period / 2) {
+                is_initial_phase = false;
+                p = 0;
+                t0 = 3 * t0 / 4;
+            }
+        }
+    }
+    BestW = minimum_one_tree(dist, best_pi).first;
+    return std::make_pair(BestW, best_pi);
+}

tsp/minimum_one_tree.hpp

@@ -0,0 +1,56 @@
+#pragma once
+
+#include <numeric>
+#include <utility>
+#include <vector>
+
+template <class T, class DistanceMatrix>
+std::pair<T, std::vector<int>>
+minimum_one_tree(const DistanceMatrix &dist, const std::vector<T> &pi) {
+
+    assert(dist.n() > 2);
+    std::vector<T> dp(dist.n(), std::numeric_limits<T>::max());
+    std::vector<int> prv(dist.n(), -1);
+    std::vector<int> used(dist.n());
+
+    auto fix_v = [&](int x) -> void {
+        dp.at(x) = std::numeric_limits<T>::max();
+        used.at(x) = 1;
+        for (auto [y, d] : dist.adjacents(x)) {
+            if (used.at(y)) continue;
+            if (T len = pi.at(x) + pi.at(y) + d; len < dp.at(y)) dp.at(y) = len, prv.at(y) = x;
+        }
+    };
+
+    T W = T();
+    std::vector<int> V(dist.n(), -2);
+
+    fix_v(0);
+    for (int t = 0; t < dist.n() - 1; ++t) {
+        int i = std::min_element(dp.cbegin(), dp.cend()) - dp.cbegin();
+        W += dp.at(i);
+        ++V.at(i);
+        ++V.at(prv.at(i));
+        fix_v(i);
+    }
+
+    // p.26, http://webhotel4.ruc.dk/~keld/research/LKH/LKH-2.0/DOC/LKH_REPORT.pdf
+    T wlo = T();
+    int ilo = -1, jlo = -1;
+    for (int i = 0; i < dist.n(); ++i) {
+        if (V.at(i) != -1) continue;
+        T tmp = T();
+        int jtmp = -1;
+        for (auto [j, d] : dist.adjacents(i)) {
+            if (prv.at(i) == j or prv.at(j) == i or i == j) continue;
+            if (T len = pi.at(i) + pi.at(j) + d; jtmp == -1 or tmp > len) tmp = len, jtmp = j;
+        }
+        if (jtmp != -1 and (ilo == -1 or wlo < tmp)) wlo = tmp, ilo = i, jlo = jtmp;
+    }
+    ++V.at(ilo);
+    ++V.at(jlo);
+
+    W += wlo - std::accumulate(pi.cbegin(), pi.cend(), T()) * 2;
+
+    return {W, V};
+}

tsp/mst_edges.hpp

@@ -0,0 +1,32 @@
+
+#pragma once
+
+#include <algorithm>
+#include <limits>
+#include <utility>
+#include <vector>
+
+template <class DistanceMatrix>
+std::vector<std::pair<int, int>> mst_edges(const DistanceMatrix &dist) {
+    using T = decltype((*dist.adjacents(0).begin()).second);
+
+    if (dist.n() <= 1) return {};
+
+    std::vector<T> dp(dist.n(), std::numeric_limits<T>::max());
+    std::vector<int> prv(dist.n(), -1);
+    std::vector<int> used(dist.n());
+    std::vector<std::pair<int, int>> ret(dist.n() - 1);
+
+    for (int t = 0; t < dist.n(); ++t) {
+        int x = std::min_element(dp.cbegin(), dp.cend()) - dp.cbegin();
+        dp.at(x) = std::numeric_limits<T>::max();
+        used.at(x) = 1;
+        if (t > 0) ret.at(t - 1) = {prv.at(x), x};
+
+        for (auto [y, len] : dist.adjacents(x)) {
+            if (!used.at(y) and len < dp.at(y)) dp.at(y) = len, prv.at(y) = x;
+        }
+    }
+
+    return ret;
+}