tsp.hpp

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <limits>
#include <py2cpp/range.hpp>
#include <py2cpp/set.hpp>
#include <set>
#include <tuple>
#include <vector>
#include <xnetwork/classes/graph.hpp>
 
// ---------------------------------------------------------------------------
// Helper: calculate_total_distance
// ---------------------------------------------------------------------------
 
template <typename Node, typename WeightFunc>
double calculate_total_distance(const std::vector<Node>& path, WeightFunc&& weight) {
    double dist = 0.0;
    for (size_t i = 0; i + 1 < path.size(); ++i) {
        dist += std::forward<WeightFunc>(weight)(path[i], path[i + 1]);
    }
    return dist;
}
 
// ---------------------------------------------------------------------------
// Internal helpers for the Christofides algorithm
// ---------------------------------------------------------------------------
 
namespace detail {
 
    template <typename Node, typename WeightFunc> auto prim_mst(size_t n, WeightFunc&& weight)
        -> std::vector<std::pair<Node, Node>> {
        constexpr double INF = std::numeric_limits<double>::max();
 
        std::vector<bool> in_mst(n, false);
        std::vector<double> key(n, INF);
        std::vector<Node> parent(n, Node{0});
 
        key[0] = 0.0;
 
        for (size_t count = 0; count < n; ++count) {
            // Find the vertex with smallest key not yet in MST
            double best = INF;
            Node u = Node{0};
            for (Node v = Node{0}; v < static_cast<Node>(n); ++v) {
                if (!in_mst[static_cast<size_t>(v)] && key[static_cast<size_t>(v)] < best) {
                    best = key[static_cast<size_t>(v)];
                    u = v;
                }
            }
            if (best == INF) break;  // isolated vertex (shouldn't happen for complete graphs)
 
            in_mst[static_cast<size_t>(u)] = true;
 
            // Relaxation
            for (Node v = Node{0}; v < static_cast<Node>(n); ++v) {
                if (!in_mst[static_cast<size_t>(v)]) {
                    const double w = std::forward<WeightFunc>(weight)(u, v);
                    auto& kv = key[static_cast<size_t>(v)];
                    if (w < kv) {
                        kv = w;
                        parent[static_cast<size_t>(v)] = u;
                    }
                }
            }
        }
 
        std::vector<std::pair<Node, Node>> edges;
        edges.reserve(n - 1);
        for (Node v = Node{1}; v < static_cast<Node>(n); ++v) {
            const auto p = parent[static_cast<size_t>(v)];
            if (p != v)  // safety
                edges.emplace_back(p, v);
        }
        return edges;
    }
 
    template <typename Node>
    auto find_odd_degree_nodes(const std::vector<std::pair<Node, Node>>& mst_edges, size_t n)
        -> std::vector<Node> {
        std::vector<int> deg(n, 0);
        for (const auto& [u, v] : mst_edges) {
            ++deg[static_cast<size_t>(u)];
            ++deg[static_cast<size_t>(v)];
        }
        std::vector<Node> odd;
        for (Node v = Node{0}; v < static_cast<Node>(n); ++v) {
            if (deg[static_cast<size_t>(v)] % 2 != 0) odd.push_back(v);
        }
        return odd;
    }
 
    template <typename Node, typename WeightFunc>
    auto greedy_min_weight_matching(const std::vector<Node>& odd_nodes, WeightFunc&& weight)
        -> std::vector<std::pair<Node, Node>> {
        const size_t k = odd_nodes.size();
        if (k < 2) return {};
 
        // Build edge list sorted by weight
        std::vector<std::tuple<double, Node, Node>> edges;
        edges.reserve(k * (k - 1) / 2);
        for (size_t i = 0; i < k; ++i) {
            for (size_t j = i + 1; j < k; ++j) {
                edges.emplace_back(std::forward<WeightFunc>(weight)(odd_nodes[i], odd_nodes[j]),
                                   odd_nodes[i], odd_nodes[j]);
            }
        }
        std::sort(edges.begin(), edges.end(),
                  [](const auto& a, const auto& b) { return std::get<0>(a) < std::get<0>(b); });
 
        std::vector<bool> used(k, false);
        std::vector<std::pair<Node, Node>> matching;
        matching.reserve(k / 2);
 
        for (const auto& [w, u, v] : edges) {
            // Map node -> index in odd_nodes (linear scan, k is small)
            size_t iu = 0;
            size_t iv = 0;
            for (size_t t = 0; t < k; ++t) {
                if (odd_nodes[t] == u) iu = t;
                if (odd_nodes[t] == v) iv = t;
            }
            if (!used[iu] && !used[iv]) {
                matching.emplace_back(u, v);
                used[iu] = true;
                used[iv] = true;
            }
        }
        return matching;
    }
 
    template <typename Node>
    auto build_multigraph(size_t n, const std::vector<std::pair<Node, Node>>& mst_edges,
                          const std::vector<std::pair<Node, Node>>& matching_edges)
        -> std::vector<std::multiset<Node>> {
        std::vector<std::multiset<Node>> adj(n);
        const auto add = [&](Node a, Node b) {
            adj[static_cast<size_t>(a)].insert(b);
            adj[static_cast<size_t>(b)].insert(a);
        };
        for (const auto& [u, v] : mst_edges) add(u, v);
        for (const auto& [u, v] : matching_edges) add(u, v);
        return adj;
    }
 
    template <typename Node> auto hierholzer(std::vector<std::multiset<Node>> adj, Node start)
        -> std::vector<std::pair<Node, Node>>;
 
    template <typename Node>
    auto shortcut_eulerian(const std::vector<std::pair<Node, Node>>& eulerian_circuit)
        -> std::vector<Node> {
        std::vector<Node> path;
        py::set<Node> visited;
        for (const auto& [u, v] : eulerian_circuit) {
            if (!visited.contains(u)) {
                path.push_back(u);
                visited.insert(u);
            }
        }
        if (!path.empty()) path.push_back(path.front());  // close the loop
        return path;
    }
 
}  // namespace detail
 
// ---------------------------------------------------------------------------
// Christofides TSP - 3/2-approximation
// ---------------------------------------------------------------------------
 
template <typename Graph, typename WeightFunc>
auto christofides_tsp(const Graph& G, WeightFunc&& weight) -> std::vector<typename Graph::node_t> {
    using Node = typename Graph::node_t;
    const size_t n = G.number_of_nodes();
 
    if (n == 0) return {};
    if (n == 1) return {Node{0}, Node{0}};
 
    // 1. Minimum Spanning Tree
    const auto mst_edges = detail::prim_mst<Node>(n, weight);
 
    // 2. Odd-degree vertices in the MST
    const auto odd_nodes = detail::find_odd_degree_nodes(mst_edges, n);
 
    // 3. Minimum-weight perfect matching on odd vertices
    const auto matching
        = detail::greedy_min_weight_matching(odd_nodes, std::forward<WeightFunc>(weight));
 
    // 4. Build the Eulerian multigraph (MST + matching)
    auto adj = detail::build_multigraph<Node>(n, mst_edges, matching);
 
    // 5. Eulerian circuit via Hierholzer
    const auto eulerian_circuit = detail::hierholzer<Node>(std::move(adj), Node{0});
 
    // 6. Shortcut -> Hamiltonian cycle
    return detail::shortcut_eulerian(eulerian_circuit);
}
 
// ---------------------------------------------------------------------------
// 2-Opt local search
// ---------------------------------------------------------------------------
 
template <typename Graph, typename WeightFunc>
auto two_opt(std::vector<typename Graph::node_t> path, const Graph& G, WeightFunc&& weight)
    -> std::vector<typename Graph::node_t> {
    (void)G;  // only used for type deduction
 
    bool improved = true;
    while (improved) {
        improved = false;
        for (size_t i = 1; i + 2 < path.size(); ++i) {
            for (size_t j = i + 1; j + 1 < path.size(); ++j) {
                if (j - i == 1) continue;
 
                const double delta = -weight(path[i - 1], path[i]) - weight(path[j], path[j + 1])
                                     + weight(path[i - 1], path[j])
                                     + std::forward<WeightFunc>(weight)(path[i], path[j + 1]);
 
                if (delta < -1.0e-12) {
                    std::reverse(path.begin() + static_cast<ptrdiff_t>(i),
                                 path.begin() + static_cast<ptrdiff_t>(j + 1));
                    improved = true;
                }
            }
        }
    }
    return path;
}
 
// ---------------------------------------------------------------------------
// Combined solver: Christofides + 2-Opt
// ---------------------------------------------------------------------------
 
template <typename Graph, typename WeightFunc>
auto solve_christofides_2opt_tsp(const Graph& G, WeightFunc&& weight)
    -> std::vector<typename Graph::node_t> {
    auto path = christofides_tsp(G, std::forward<WeightFunc>(weight));
    if (path.size() <= 3) return path;
    return two_opt(std::move(path), G, weight);
}