test/Heuristic/approx_tsp_solver/my_test.cpp

View this file on GitHub
Last update: 2025-07-28 00:16:46+09:00

Depends on

Code

// competitive-verifier: STANDALONE

#include "Heuristic/approx_tsp_solver.hpp"
#include <bits/stdc++.h>

using namespace std;

long double solve(const int N, const vector<int> &X, const vector<int> &Y)
{
    approx_tsp::ApproxTspSolver<approx_tsp::Euclidean, approx_tsp::PathMode::Closed> tsp_solver;
    vector<array<int, 2>> points(X.size());
    for (int i = 0; i < X.size(); i++)
    {
        points[i] = {X[i], Y[i]};
    }

    auto ans = tsp_solver.solve(points);
    return ans.first;
}

void solve1()
{
    const int N = 4;
    const vector<int> X = {0, 0, 1, 1}, Y = {0, 1, 0, 1};
    const long double my_ans = solve(N, X, Y);
    const long double true_ans = 4.000000000000;

    assert(abs(my_ans - true_ans) < 0.001);
}

void solve2()
{
    const int N = 7;
    const vector<int> X = {2, 2, 4, 1, 7, 5, 6}, Y = {5, 3, 1, 1, 2, 3, 5};
    const long double my_ans = solve(N, X, Y);
    const long double true_ans = 18.870481592668;

    assert(abs(my_ans - true_ans) < 0.001);
}

void solve3()
{
    const int N = 11;
    const vector<int> X = {642, 16, 898, 264, 751, 92, 963, 152, 333, 632, 27};
    const vector<int> Y = {364, 717, 53, 824, 558, 496, 277, 618, 743, 559, 40};
    const long double my_ans = solve(N, X, Y);
    const long double true_ans = 3270.811560094203;

    assert(abs(my_ans - true_ans) < 0.001);
}

void solve4()
{
    const int N = 15;
    const vector<int> X = {776, 218, 839, 975, 192, 90, 493, 285, 441, 445, 560, 784, 570, 982, 452};
    const vector<int> Y = {273, 998, 577, 670, 465, 329, 586, 494, 175, 612, 777, 266, 778, 130, 599};
    const long double my_ans = solve(N, X, Y);
    const long double true_ans = 3441.451976970622;

    assert(abs(my_ans - true_ans) < 0.001);
}

int main()
{

    solve1();
    solve2();
    solve3();
    solve4();

    return 0;
}

#line 1 "test/Heuristic/approx_tsp_solver/my_test.cpp"
// competitive-verifier: STANDALONE

#line 1 "Heuristic/approx_tsp_solver.hpp"
#include <bits/stdc++.h>
using namespace std;

#line 2 "Heuristic/fast_random_engine.hpp"
using namespace std;

class FastRandomEngine
{
private:
    // 0以上UINT_MAX(0xffffffff)以下の整数をとる乱数 xorshift https://ja.wikipedia.org/wiki/Xorshift
    static uint32_t randXor()
    {
        static uint32_t x = 123456789;
        static uint32_t y = 362436069;
        static uint32_t z = 521288629;
        static uint32_t w = 88675123;
        uint32_t t;

        t = x ^ (x << 11);
        x = y;
        y = z;
        z = w;
        return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
    }

    // 0以上1未満の小数をとる乱数
    static double rand01()
    {
        return (randXor() + 0.5) * (1.0 / UINT_MAX);
    }

public:
    using result_type = uint32_t;

    FastRandomEngine(uint32_t seed = 0) : state(seed) {}

    uint32_t operator()()
    {
        return randXor();
    }

    static constexpr uint32_t min()
    {
        return std::numeric_limits<uint32_t>::min();
    }

    static constexpr uint32_t max()
    {
        return std::numeric_limits<uint32_t>::max();
    }

    static long long random_int(pair<long long, long long> min_max)
    {
        return random_int(min_max.first, min_max.second);
    }

    static long long random_int(long long min, long long max)
    {
        long long width = max - min + 1;

        return randXor() % width + min;
    }

    static long double random_real(pair<long double, long double> min_max)
    {
        return random_real(min_max.first, min_max.second);
    }

    static long double random_real(long double min, long double max)
    {
        long double width = max - min;

        return rand01() * width + min;
    }

private:
    uint32_t state;
};
FastRandomEngine rand_engine(0);

template <class T>
static void shuffle_vector(vector<T> &v)
{
    shuffle(v.begin(), v.end(), rand_engine);
}
#line 5 "Heuristic/approx_tsp_solver.hpp"

namespace approx_tsp
{
    // ---------------------
    // 距離関数オブジェクト
    // ---------------------
    struct Euclidean
    {
        using CostType = double;
        template <typename T, size_t D>
        double operator()(const array<T, D> &a, const array<T, D> &b) const
        {
            CostType sum = 0;
            for (size_t i = 0; i < D; i++)
            {
                CostType d = a[i] - b[i];
                sum += d * d;
            }
            return sqrt(sum);
        }
    };

    struct Manhattan
    {
        using CostType = int;
        template <typename T, size_t D>
        CostType operator()(const array<T, D> &a, const array<T, D> &b) const
        {
            CostType sum = 0;
            for (size_t i = 0; i < D; i++)
            {
                sum += abs(a[i] - b[i]);
            }
            return sum;
        }
    };

    // ---------------------
    // パスモード列挙型
    // ---------------------
    enum class PathMode
    {
        Closed, // 1周する必要がある 1 → 2　→ ...　→ N　→ 1
        Open    // 1周する必要がない 1 → 2　→ ...　→ N
    };

    // 参考 https://future-architect.github.io/articles/20211201a/
    template <
        typename Distance = Euclidean,
        PathMode Mode = PathMode::Closed>
    class ApproxTspSolver
    {
        using CostType = Distance::CostType;
        Distance dist;

        vector<vector<pair<CostType, int>>> cost_table;

        CostType double_bridge(vector<int> &nexts, vector<int> &prevs, const vector<vector<CostType>> &costs)
        {
            int N = nexts.size();
            if (nexts.size() < 8)
            {
                return 0;
            }

            static vector<int> random_ids;
            if (random_ids.size() != nexts.size())
            {
                random_ids = nexts;
            }

            shuffle_vector(random_ids);

            unordered_set<int> used;
            array<array<int, 2>, 4> bridge_edges;

            int bridge_cnt = 0;
            // ダブルブリッジするidを8つ選ぶ。 これは隣接する2項を4つ選ぶ。
            for (auto curr : random_ids)
            {
                int next = nexts[curr];

                if (used.contains(curr) || used.contains(next))
                {
                    continue;
                }

                used.insert(curr), used.insert(next);

                bridge_edges[bridge_cnt++] = {curr, next};
                if (bridge_cnt == 4)
                {
                    break;
                }
            }

            if (bridge_cnt != 4)
            {
                return 0;
            }

            vector<int> order(N);
            int now = 0;
            for (int i = 0; i < N; i++)
            {
                order[now] = i;
                now = nexts[now];
            }

            sort(bridge_edges.begin(), bridge_edges.end(), [&](const array<int, 2> &a, const array<int, 2> &b) -> bool
                 { return order[a[0]] < order[b[0]]; });

            CostType diff_cost = 0;
            for (int i = 0; i < 4; i++)
            {
                diff_cost -= costs[bridge_edges[i][0]][bridge_edges[i][1]];
            }

            for (int i = 0; i < 4; i++)
            {
                nexts[bridge_edges[i][0]] = bridge_edges[(i + 2) % 4][1];
                prevs[bridge_edges[(i + 2) % 4][1]] = bridge_edges[i][0];
            }

            for (int i = 0; i < 4; i++)
            {
                diff_cost += costs[bridge_edges[i][0]][bridge_edges[(i + 2) % 4][1]];
            }

            return diff_cost;
        }

        // l ~ r までの向きを逆にする(lとrは含まない)
        void reverse_node(vector<int> &nexts, vector<int> &prevs, int l, int r)
        {
            nexts[r] = prevs[r];
            int now = nexts[r];
            while (now != l)
            {
                swap(nexts[now], prevs[now]);
                now = nexts[now];
            }
            prevs[l] = nexts[l];
        }

        CostType two_opt(vector<int> &nexts, vector<int> &prevs, const vector<vector<CostType>> &costs)
        {
            int N = nexts.size();
            if (nexts.size() < 4)
            {
                return 0;
            }

            CostType diff_cost = 0;

            auto try_swap = [&](int u1, int u2) -> bool
            {
                for (auto [cost, u3] : cost_table[u1])
                {
                    int u4 = nexts[u3];
                    if (u2 == u3)
                    {
                        break;
                    }
                    if (costs[u1][u2] + costs[u3][u4] <= costs[u1][u3] + costs[u2][u4])
                    {
                        continue;
                    }

                    diff_cost -= costs[u1][u2], diff_cost -= costs[u3][u4];
                    diff_cost -= costs[u2][nexts[u2]], diff_cost -= costs[u3][prevs[u3]];

                    reverse_node(nexts, prevs, u2, u3);
                    nexts[u1] = u3, prevs[u3] = u1;
                    nexts[u2] = u4, prevs[u4] = u2;

                    diff_cost += costs[u1][u3], diff_cost += costs[u2][u4];
                    diff_cost += costs[u3][nexts[u3]], diff_cost += costs[u2][prevs[u2]];

                    return true;
                }

                for (auto [cost, u4] : cost_table[u2])
                {
                    int u3 = prevs[u4];
                    if (u1 == u4)
                    {
                        break;
                    }
                    if (costs[u1][u2] + costs[u3][u4] <= costs[u1][u3] + costs[u2][u4])
                    {
                        continue;
                    }

                    diff_cost -= costs[u1][u2], diff_cost -= costs[u3][u4];
                    diff_cost -= costs[u2][nexts[u2]], diff_cost -= costs[u3][prevs[u3]];

                    reverse_node(nexts, prevs, u2, u3);
                    nexts[u1] = u3, prevs[u3] = u1;
                    nexts[u2] = u4, prevs[u4] = u2;

                    diff_cost += costs[u1][u3], diff_cost += costs[u2][u4];
                    diff_cost += costs[u3][nexts[u3]], diff_cost += costs[u2][prevs[u2]];

                    return true;
                }

                return false;
            };

            while (true)
            {
                bool update = false;
                for (int u1 = 0; u1 < N; u1++)
                {
                    int u2 = nexts[u1];

                    if (try_swap(u1, u2))
                    {
                        update = true;
                        break;
                    }
                }

                if (!update)
                {
                    break;
                }
            }

            return diff_cost;
        }

    public:
        template <typename PointType, size_t D>
        pair<CostType, vector<int>> solve(const vector<array<PointType, D>> &points)
        {
            const int N = points.size();
            const int M = (Mode == PathMode::Closed ? N : N + 1);

            vector costs(M, vector<CostType>(M));

            for (int i = 0; i < N; i++)
            {
                for (int j = 0; j < N; j++)
                {
                    if (i == j)
                    {
                        continue;
                    }

                    costs[i][j] = dist(points[i], points[j]);
                }
            }

            return solve(costs);
        }

        pair<CostType, vector<int>> solve(const vector<vector<CostType>> &costs)
        {
            int N = costs.size();

            // i番目の前後の頂点が何番かを管理する
            vector<int> nexts(N), prevs(N);
            for (int i = 0; i < N; i++)
            {
                nexts[i] = (i + 1) % N;
                prevs[i] = (i - 1 + N) % N;
            }

            cost_table = vector(N, vector<pair<CostType, int>>());
            for (int i = 0; i < N; i++)
            {
                for (int j = 0; j < N; j++)
                {
                    if (i == j)
                    {
                        continue;
                    }
                    cost_table[i].push_back(pair(costs[i][j], j));
                }
                sort(cost_table[i].begin(), cost_table[i].end());
            }

            vector<int> global_nexts = nexts, global_prevs = prevs;
            for (int i = 0; i < 100; i++)
            {
                nexts = global_nexts, prevs = global_prevs;

                CostType diff_score = 0;
                diff_score += double_bridge(nexts, prevs, costs);
                diff_score += two_opt(nexts, prevs, costs);

                if (diff_score < 0)
                {
                    global_nexts = nexts, global_prevs = prevs;
                }
            }

            vector<int> result_path = {global_nexts.back()};
            while (result_path.back() != N - 1)
            {
                result_path.push_back(global_nexts[result_path.back()]);
            }

            if constexpr (Mode == PathMode::Open)
            {
                result_path.pop_back();
            }

            CostType result_cost = 0;
            for (int i = 0; i < result_path.size() - 1; i++)
            {
                result_cost += costs[result_path[i]][result_path[i + 1]];
            }

            if constexpr (Mode == PathMode::Closed)
            {
                result_cost += costs[result_path.back()][result_path.front()];
            }

            return pair(result_cost, result_path);
        }
    };
}
#line 5 "test/Heuristic/approx_tsp_solver/my_test.cpp"

using namespace std;

long double solve(const int N, const vector<int> &X, const vector<int> &Y)
{
    approx_tsp::ApproxTspSolver<approx_tsp::Euclidean, approx_tsp::PathMode::Closed> tsp_solver;
    vector<array<int, 2>> points(X.size());
    for (int i = 0; i < X.size(); i++)
    {
        points[i] = {X[i], Y[i]};
    }

    auto ans = tsp_solver.solve(points);
    return ans.first;
}

void solve1()
{
    const int N = 4;
    const vector<int> X = {0, 0, 1, 1}, Y = {0, 1, 0, 1};
    const long double my_ans = solve(N, X, Y);
    const long double true_ans = 4.000000000000;

    assert(abs(my_ans - true_ans) < 0.001);
}

void solve2()
{
    const int N = 7;
    const vector<int> X = {2, 2, 4, 1, 7, 5, 6}, Y = {5, 3, 1, 1, 2, 3, 5};
    const long double my_ans = solve(N, X, Y);
    const long double true_ans = 18.870481592668;

    assert(abs(my_ans - true_ans) < 0.001);
}

void solve3()
{
    const int N = 11;
    const vector<int> X = {642, 16, 898, 264, 751, 92, 963, 152, 333, 632, 27};
    const vector<int> Y = {364, 717, 53, 824, 558, 496, 277, 618, 743, 559, 40};
    const long double my_ans = solve(N, X, Y);
    const long double true_ans = 3270.811560094203;

    assert(abs(my_ans - true_ans) < 0.001);
}

void solve4()
{
    const int N = 15;
    const vector<int> X = {776, 218, 839, 975, 192, 90, 493, 285, 441, 445, 560, 784, 570, 982, 452};
    const vector<int> Y = {273, 998, 577, 670, 465, 329, 586, 494, 175, 612, 777, 266, 778, 130, 599};
    const long double my_ans = solve(N, X, Y);
    const long double true_ans = 3441.451976970622;

    assert(abs(my_ans - true_ans) < 0.001);
}

int main()
{

    solve1();
    solve2();
    solve3();
    solve4();

    return 0;
}