tsp.cpp

// tsp.cpp
// -----------------------------------------------------------------
// TSP（旅行推銷員問題）：
//   1. Held–Karp 位元 DP：精確解 O(2^n · n^2)，含路徑還原
//   2. 最近鄰啟發式（Nearest Neighbor）：O(n^2) 快速建構
//   3. 2-opt 局部搜尋：反覆拆兩邊重接，改善 NN 的解
//
// Compile: g++ -std=c++17 -O2 -Wall -Wextra -o tsp tsp.cpp
// -----------------------------------------------------------------
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <random>
#include <chrono>
#include <iomanip>

using namespace std;

using Matrix = vector<vector<double>>;

// ---------- 1. Held–Karp ----------
// dp[S][v] = 從城市 0 出發、恰好經過集合 S、目前在 v 的最短距離
pair<double, vector<int>> heldKarp(const Matrix& d) {
    int n = (int)d.size();
    const double INF = 1e18;
    vector<vector<double>> dp(1 << n, vector<double>(n, INF));
    vector<vector<int>> par(1 << n, vector<int>(n, -1));
    dp[1][0] = 0;                                       // 從 0 出發
    for (int S = 1; S < (1 << n); ++S) {
        if (!(S & 1)) continue;                         // 必含起點 0
        for (int v = 0; v < n; ++v) {
            if (dp[S][v] >= INF) continue;
            for (int u = 0; u < n; ++u) {
                if (S & (1 << u)) continue;
                int S2 = S | (1 << u);
                double nd = dp[S][v] + d[v][u];
                if (nd < dp[S2][u]) {
                    dp[S2][u] = nd;
                    par[S2][u] = v;
                }
            }
        }
    }
    int full = (1 << n) - 1;
    double best = INF;
    int last = -1;
    for (int v = 1; v < n; ++v)                          // 回到 0 收尾
        if (dp[full][v] + d[v][0] < best) {
            best = dp[full][v] + d[v][0];
            last = v;
        }
    vector<int> tour;                                    // 還原
    int S = full, v = last;
    while (v != -1) {
        tour.push_back(v);
        int p = par[S][v];
        S &= ~(1 << v);
        v = p;
    }
    reverse(tour.begin(), tour.end());                   // 0 ... last
    tour.push_back(0);                                   // 閉合
    return {best, tour};
}

// ---------- 2. 最近鄰啟發式 ----------
pair<double, vector<int>> nearestNeighbor(const Matrix& d, int start = 0) {
    int n = (int)d.size();
    vector<bool> used(n, false);
    vector<int> tour = {start};
    used[start] = true;
    double len = 0;
    for (int step = 1; step < n; ++step) {
        int cur = tour.back(), next = -1;
        for (int v = 0; v < n; ++v)
            if (!used[v] && (next == -1 || d[cur][v] < d[cur][next])) next = v;
        used[next] = true;
        len += d[cur][next];
        tour.push_back(next);
    }
    len += d[tour.back()][start];
    tour.push_back(start);
    return {len, tour};
}

// ---------- 3. 2-opt 局部搜尋 ----------
// 拆掉 (t[i],t[i+1]) 和 (t[j],t[j+1])，把中段反轉重接；有改善就做
pair<double, vector<int>> twoOpt(const Matrix& d, vector<int> tour) {
    int n = (int)tour.size() - 1;                        // 尾巴是回到起點
    bool improved = true;
    while (improved) {
        improved = false;
        for (int i = 0; i < n - 1; ++i)
            for (int j = i + 1; j < n; ++j) {
                int a = tour[i], b = tour[i + 1];
                int c = tour[j], e = tour[j + 1];
                double delta = d[a][c] + d[b][e] - d[a][b] - d[c][e];
                if (delta < -1e-10) {
                    reverse(tour.begin() + i + 1, tour.begin() + j + 1);
                    improved = true;
                }
            }
    }
    double len = 0;
    for (size_t i = 0; i + 1 < tour.size(); ++i)
        len += d[tour[i]][tour[i + 1]];
    return {len, tour};
}

static Matrix randomEuclidean(int n, mt19937& rng,
                              vector<pair<double,double>>* ptsOut = nullptr) {
    uniform_real_distribution<double> coord(0, 100);
    vector<pair<double,double>> pts;
    for (int i = 0; i < n; ++i) pts.push_back({coord(rng), coord(rng)});
    Matrix d(n, vector<double>(n, 0));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            d[i][j] = hypot(pts[i].first - pts[j].first,
                            pts[i].second - pts[j].second);
    if (ptsOut) *ptsOut = pts;
    return d;
}

static void printTour(const vector<int>& t) {
    for (size_t i = 0; i < t.size(); ++i)
        cout << t[i] << (i + 1 < t.size() ? "-" : "\n");
}

int main() {
    cout << fixed << setprecision(2);

    // ---- 範例 1：手算小例（4 城市） ----
    cout << "=== Example 1: 4-city worked example ===\n";
    Matrix m4 = {{0, 10, 15, 20},
                 {10, 0, 35, 25},
                 {15, 35, 0, 30},
                 {20, 25, 30, 0}};
    auto [opt4, tour4] = heldKarp(m4);
    cout << "Held-Karp optimal = " << opt4 << " (expect 80)  tour: ";
    printTour(tour4);
    auto [nn4, nnTour4] = nearestNeighbor(m4);
    cout << "Nearest neighbor  = " << nn4 << "  tour: ";
    printTour(nnTour4);
    cout << "\n";

    // ---- 範例 2：NN 可以多糟 + 2-opt 補救 ----
    cout << "=== Example 2: NN vs NN+2opt vs optimal (n=12, Euclidean) ===\n";
    mt19937 rng(20260724);
    Matrix d12 = randomEuclidean(12, rng);
    auto t0 = chrono::steady_clock::now();
    auto [opt, optTour] = heldKarp(d12);
    auto t1 = chrono::steady_clock::now();
    auto [nnLen, nnTour] = nearestNeighbor(d12);
    auto [toLen, toTour] = twoOpt(d12, nnTour);
    auto t2 = chrono::steady_clock::now();
    cout << "optimal (Held-Karp): " << opt << "  ("
         << chrono::duration<double, milli>(t1 - t0).count() << " ms)\n";
    cout << "nearest neighbor   : " << nnLen
         << "  (+" << 100 * (nnLen / opt - 1) << "% above optimal)\n";
    cout << "NN + 2-opt         : " << toLen
         << "  (+" << 100 * (toLen / opt - 1) << "% above optimal, "
         << chrono::duration<double, milli>(t2 - t1).count() << " ms)\n";
    cout << "optimal tour: ";
    printTour(optTour);
    cout << "2-opt tour  : ";
    printTour(toTour);
    (void)optTour;
    cout << "\n";

    // ---- 範例 3：多起點 NN（廉價的改善技巧） ----
    cout << "=== Example 3: multi-start nearest neighbor ===\n";
    double bestNN = 1e18;
    int bestStart = 0;
    for (int s = 0; s < 12; ++s) {
        auto [len, t] = nearestNeighbor(d12, s);
        if (len < bestNN) { bestNN = len; bestStart = s; }
        (void)t;
    }
    cout << "best of 12 starts: " << bestNN << " (start=" << bestStart
         << ", single-start was " << nnLen << ")\n\n";

    // ---- 範例 4：規模牆——Held-Karp 的極限 ----
    cout << "=== Example 4: the scaling wall ===\n";
    for (int n : {13, 16, 19}) {
        Matrix d = randomEuclidean(n, rng);
        auto s0 = chrono::steady_clock::now();
        auto [o, ot] = heldKarp(d);
        auto s1 = chrono::steady_clock::now();
        auto [nl, nt] = nearestNeighbor(d);
        auto [tl, tt] = twoOpt(d, nt);
        auto s2 = chrono::steady_clock::now();
        (void)ot; (void)tt;
        cout << "n=" << n
             << " | HK " << chrono::duration<double, milli>(s1 - s0).count()
             << " ms (opt=" << o << ")"
             << " | NN+2opt " << chrono::duration<double, milli>(s2 - s1).count()
             << " ms (len=" << tl << ", gap +"
             << 100 * (tl / o - 1) << "%)\n";
        (void)nl;
    }
    cout << "Held-Karp memory: 2^n * n doubles -> n=25 needs ~6.7 GB. Game over.\n";
    return 0;
}