#783735 (C++17) No.5007 Steiner Space Travel

提出ソース
結果

問題	No.5007 Steiner Space Travel
ユーザー	qLethon
提出日時	2022-08-02 02:56:59
言語	C++17 (gcc 12.3.0 + boost 1.83.0)
結果	AC
実行時間	988 ms / 1,000 ms
コード長	15,942 bytes
コンパイル時間	3,451 ms
実行使用メモリ	6,952 KB
スコア	8,458,939
最終ジャッジ日時	2022-08-02 02:57:36
合計ジャッジ時間	35,303 ms
ジャッジサーバーID （参考情報）	judge13 / judge12
	純コード判定しない問題か言語
このコードへのチャレンジ
（要ログイン）
テストケース

テストケース表示
入力	結果	実行時間実行使用メモリ
testcase_00	AC	987 ms 4,904 KB
testcase_01	AC	985 ms 4,904 KB
testcase_02	AC	985 ms 4,900 KB
testcase_03	AC	988 ms 4,904 KB
testcase_04	AC	987 ms 4,904 KB
testcase_05	AC	986 ms 4,904 KB
testcase_06	AC	985 ms 4,900 KB
testcase_07	AC	987 ms 6,948 KB
testcase_08	AC	985 ms 4,900 KB
testcase_09	AC	986 ms 6,952 KB
testcase_10	AC	987 ms 4,900 KB
testcase_11	AC	988 ms 4,904 KB
testcase_12	AC	987 ms 4,900 KB
testcase_13	AC	988 ms 4,900 KB
testcase_14	AC	985 ms 4,900 KB
testcase_15	AC	986 ms 4,904 KB
testcase_16	AC	987 ms 4,904 KB
testcase_17	AC	987 ms 4,900 KB
testcase_18	AC	987 ms 4,900 KB
testcase_19	AC	986 ms 4,900 KB
testcase_20	AC	985 ms 6,948 KB
testcase_21	AC	987 ms 4,900 KB
testcase_22	AC	986 ms 4,904 KB
testcase_23	AC	986 ms 4,908 KB
testcase_24	AC	987 ms 6,948 KB
testcase_25	AC	986 ms 4,904 KB
testcase_26	AC	986 ms 4,904 KB
testcase_27	AC	986 ms 4,900 KB
testcase_28	AC	988 ms 4,900 KB
testcase_29	AC	986 ms 4,904 KB
権限があれば一括ダウンロードができます
ソースコード

raw source code
#include <bits/stdc++.h>

using namespace std;

u_int64_t seed = 91;
u_int64_t xor_shift(){
    seed ^= seed << 1;
    seed ^= seed >> 3;
    seed ^= seed << 10;
    return seed;
}

class Time{
    public:
    chrono::system_clock::time_point start_time;
    Time(){
        start_time = now();
    }

    chrono::system_clock::time_point now(){
        return chrono::system_clock::now();
    }

    int elapsed(){
        auto cur = chrono::system_clock::now();
        return chrono::duration_cast<chrono::milliseconds>(cur - start_time).count();
    }
};

template <class T>
class TSP{
    vector<vector<T>> D;
    vector<int> ans, A, v_to_cycle;
    vector<vector<int>> nearest;
    int n;
    T INF;
    public:
    TSP (vector<vector<T>> D, T INF): D(D), n(D.size()), INF(INF){
        ans.resize(n);
        iota(ans.begin(), ans.end(), 0);
        v_to_cycle.resize(n);
        iota(v_to_cycle.begin(), v_to_cycle.end(), 0);
        A.resize(n);
        iota(A.begin(), A.end(), 0);

        nearest.assign(n, vector<int>(n));
        for (int i = 0; i < n; i++){
            vector<pair<T, int>> B(n);
            for (int j = 0; j < n; j++)
                B[j] = {D[i][j], j};
            sort(B.begin(), B.end());
            for (int j = 0; j < n; j++)
                nearest[i][j] = B[j].second;
        }
    };

    vector<int> solve_perm(){
        vector<int> P(n);
        iota(P.begin(), P.end(), 0);

        T best_score = INF;
        vector<int> best = P;
        do {
            T s = calc_score(P);
            if (s < best_score){
                best = P;
                best_score = s;
            }
        } while (next_permutation(P.begin(), P.end()));
    
        return best;
    }
    
    vector<int> solve(int time_limit){
        if (n < 10)
            return solve_perm();

        Time time;

        // // initial 2-opt
        // while (true){
        //     bool updated = false;
        //     for (int t1 = 0; t1 < n; t1++){
        //         int t2 = (t1 + 1) % n;
        //         for (const auto &v: nearest[ans[t2]]){
        //             int t3 = v_to_cycle[v];
        //             if (dist(t1, t2) <= dist(t2, t3))
        //                 break;
        //             int t4 = (t3 - 1 + n) % n;
        //             int i, j;
        //             if (t2 <= t4)
        //                 tie(i, j) = {t2, t4};
        //             else
        //                 tie(i, j) = {t3, t1};
                    
        //             const auto d = calc_two_opt_score(i, j);
        //             if (d < 0){
        //                 two_opt_swap(i, j);
        //                 updated = true;
        //             }
        //         }
        //     }
        //     if (!updated)
        //         break;
        // }

        // kick (double bridge + 2-opt)
        auto best = ans;
        auto best_v_to_cycle = v_to_cycle;
        T score = calc_score(ans);
        int loop = 0;

        while (time.elapsed() < time_limit){
            loop++;
            ans = best;
            v_to_cycle = best_v_to_cycle;

            // kick
            T diff = 0;
            {
            auto chosen = choose_k(4);
            sort(chosen.begin(), chosen.end());
            diff = calc_double_bridge(chosen[0], chosen[1], chosen[2], chosen[3]);
            double_birdge(chosen[0], chosen[1], chosen[2], chosen[3]);
            }

            // // random 2-opt
            // for (int i = 0; i < 6; i++){
            //     const int a = xor_shift() % n;
            //     const int b = xor_shift() % n;
            //     diff += calc_two_opt_score(a, b);
            //     two_opt_swap(a, b);
            // }
            
            vector<int> Q(n);
            iota(Q.begin(), Q.end(), 0);

            // local optimization
            while (true){
                bool updated = false;
                vector<int> NQ;
                for (auto u: Q){
                    int t1 = v_to_cycle[u];
                    int t2 = (t1 + 1) % n;
                    for (const auto &v: nearest[ans[t2]]){
                        int t3 = v_to_cycle[v];
                        if (dist(t1, t2) <= dist(t2, t3))
                            break;
                        int t4 = (t3 - 1 + n) % n;
                        int i, j;
                        if (t2 <= t4)
                            tie(i, j) = {t2, t4};
                        else
                            tie(i, j) = {t3, t1};
                        
                        const auto d = calc_two_opt_score(i, j);
                        if (d < 0){
                            NQ.push_back(ans[(t1 - 1 + n) % n]);
                            NQ.push_back(ans[t1]);
                            NQ.push_back(ans[t2]);
                            NQ.push_back(ans[(t2 + 1) % n]);
                            NQ.push_back(ans[(t4 - 1 + n) % n]);
                            NQ.push_back(ans[t3]);
                            NQ.push_back(ans[t4]);
                            NQ.push_back(ans[(t3 + 1) % n]);
                            two_opt_swap(i, j);
                            diff += d;
                            updated = true;
                            break;
                        }
                    }
                }
                if (!updated)
                    break;
                Q = NQ;
                sort(Q.begin(), Q.end());
                Q.erase(unique(Q.begin(), Q.end()), Q.end());
            }

            if (diff < 0){
                best = ans;
                best_v_to_cycle = v_to_cycle;
                score += diff;
                // cerr << score << endl;
            }
        }

        // cerr << "loop: " << loop << endl;
        // cerr << score << endl;
        return best;
    }

    T dist(int i, int j){
        // こっちに%n押し付けたいだろ
        if (i < 0)
            i += n;
        if (i >= n)
            i -= n;
        if (j < 0)
            j += n;
        if (j >= n)
            j -= n;
        return D[ans[i]][ans[j]];
    }

    T dist(int i, int j, const vector<int> &ans){
        return D[ans[i]][ans[j]];
    }

    T calc_score(const vector<int> &ans){
        T res = 0;
        assert(ans.size() > 0);
        for (int i = 0; i < ans.size() - 1; i++)
            res += dist(i, i + 1, ans);
        res += dist(ans.size() - 1, 0, ans);
        return res;
    }

    T calc_two_opt_score(int i, int j){
        T res = 0;
        if (i > j)
            swap(i, j);
        if (i == 0 and j == n - 1)
            return 0;

        res += dist((i - 1 + n) % n, j) - dist((i - 1 + n) % n, i);
        res += dist(i, (j + 1) % n) - dist(j, (j + 1) % n);
        
        return res;
    }

    void two_opt_swap(int i, int j){
        if (i > j)
            swap(i, j);
        reverse(ans.begin() + i, ans.begin() + j + 1);
        for (int a = i; a <= j; a++)
            v_to_cycle[ans[a]] = a;
    }

    T three_opt_swap(int i, int j, int k){
        assert(i < j and j < k);
        array<T, 5> d;

        i--; j--; k--;  // バカ

        d[0] = dist(i, i + 1) + dist(j, j + 1) + dist(k, k + 1);
        d[1] = dist(i, j) + dist(i + 1, j + 1) + dist(k, k + 1);
        d[2] = dist(i, i + 1) + dist(j, k) + dist(j + 1, k + 1);
        d[3] = dist(i, j + 1) + dist(k, i + 1) + dist(j, k + 1);
        d[4] = dist(k + 1, i + 1) + dist(j, j + 1) + dist(k, i);

        i++; j++; k++;

        if (d[0] > d[1]){
            reverse(ans.begin() + i, ans.begin() + j);
            return -d[0] + d[1];
        } else if (d[0] > d[2]){
            reverse(ans.begin() + j, ans.begin() + k);
            return -d[0] + d[2];
        } else if (d[0] > d[4]){
            reverse(ans.begin() + i, ans.begin() + k);
            return -d[0] + d[4];
        } else if (d[0] > d[3]){
            vector<int> tmp(k - i);
            copy(ans.begin() + j, ans.begin() + k, tmp.begin());
            copy(ans.begin() + i, ans.begin() + j, tmp.begin() + k - j);
            copy(tmp.begin(), tmp.end(), ans.begin() + i);
            return -d[0] + d[3];
        }
        return 0;
    }

    T calc_double_bridge(int i, int j, int k, int l){
        assert(i < j and j < k and k < l);
        T res = 0;
        res -= dist(i, i + 1) + dist(j, j + 1) + dist(k, k + 1) + dist(l, (l + 1) % n);
        res += dist(i, k + 1) + dist(j, (l + 1) % n) + dist(k, i + 1) + dist(l, j + 1);

        return res;
    }

    void double_birdge(int i, int j, int k, int l){
        assert(i < j and j < k and k < l);
        assert(l < n);
        vector<int> res(n);
        copy(ans.begin(), ans.begin() + i + 1, res.begin());
        copy(ans.begin() + k + 1, ans.begin() + l + 1, res.begin() + i + 1);
        copy(ans.begin() + j + 1, ans.begin() + k + 1, res.begin() + i + l - k + 1);
        copy(ans.begin() + i + 1, ans.begin() + j + 1, res.begin() + i + l - j + 1);
        copy(ans.begin() + l + 1, ans.end(), res.begin() + l + 1);
        ans = res;

        for (int a = 0; a < n; a++)
            v_to_cycle[ans[a]] = a;
    }

    vector<int> choose_k(int k){
        const int m = A.size();
        assert(k <= m);
        vector<int> res(k);
        for (int i = 0; i < k; i++){
            int r = i + xor_shift() % (m - i);
            swap(A[i], A[r]);
            res[i] = A[i];
        }

        return res;
    }
};

template <class T>
pair<vector<T>, vector<T>> k_means(vector<T> X, vector<T> Y, int k){
    const int n = X.size();
    vector<int> C(n), cluster_size(k);

    for (auto &c: C){
        c = xor_shift() % k;
        cluster_size[c]++;
    }
    
    vector<int> prev_C;
    do {
        vector<T> vx(k), vy(k);
        for (int i = 0; i < n; i++){
            vx[C[i]] += X[i];
            vy[C[i]] += Y[i];
        }
        for (int i = 0; i < k; i++){
            if (cluster_size[i] == 0){
                vx[i] = xor_shift() % 1000;
                vy[i] = xor_shift() % 1000;
                continue;
            }
            vx[i] /= cluster_size[i];
            vy[i] /= cluster_size[i];
        }

        prev_C = C;
        for (int i = 0; i < n; i++){
            T mi = (vx[0] - X[i]) * (vx[0] - X[i]) + (vy[0] - Y[i]) * (vy[0] - Y[i]);
            int mixi = 0;
            for (int j = 1; j < k; j++){
                T d = (vx[j] - X[i]) * (vx[j] - X[i]) + (vy[j] - Y[i]) * (vy[j] - Y[i]);
                if (d < mi){
                    mi = d;
                    mixi = j;
                }
            }
            C[i] = mixi;
        }
        fill(cluster_size.begin(), cluster_size.end(), 0);
        for (const auto &c: C)
            cluster_size[c]++;
    } while (C != prev_C);

    vector<T> vx(k), vy(k);
    for (int i = 0; i < n; i++){
        vx[C[i]] += X[i];
        vy[C[i]] += Y[i];
    }
    for (int i = 0; i < k; i++){
        if (cluster_size[i] == 0)
            continue;
        vx[i] /= cluster_size[i];
        vy[i] /= cluster_size[i];
    }

    return {vx, vy};
}

int n, m;
vector<int64_t> A, B;
constexpr int64_t alpha = 5;

vector<vector<int64_t>> calc_dist(vector<int64_t> vx, vector<int64_t> vy){
    vector<vector<int64_t>> D(n + m, vector<int64_t>(n + m));
    for (int i = 0; i < n; i++){
        for (int j = 0; j < n; j++)
            D[j][i] = D[i][j] = ((A[i] - A[j]) * (A[i] - A[j]) + (B[i] - B[j]) * (B[i] - B[j])) * alpha * alpha;
        for (int j = 0; j < m; j++)
            D[n + j][i] = D[i][n + j] = ((A[i] - vx[j]) * (A[i] - vx[j]) + (B[i] - vy[j]) * (B[i] - vy[j])) * alpha;
    }

    for (int i = 0; i < m; i++){
        for (int j = 0; j < m; j++)
            D[n + j][n + i] = D[n + i][n + j] = (vx[i] - vx[j]) * (vx[i] - vx[j]) + (vy[i] - vy[j]) * (vy[i] - vy[j]);
    }

    for (int k = 0; k < n + m; k++){
        for (int i = 0; i < n + m; i++){
            for (int j = 0; j < n + m; j++)
                D[i][j] = min(D[i][j], D[i][k] + D[k][j]);
        }
    }

    vector<vector<int64_t>> res(n, vector<int64_t>(n));
    for (int i = 0; i < n; i++){
        for (int j = 0; j < n; j++)
            res[i][j] = D[i][j];
    }

    return res;
}

pair<int, vector<int>> find_best_route(int start, int goal, const vector<vector<int64_t>> &D){
    using state = array<int64_t, 3>;  // d, v
    priority_queue<state, vector<state>, greater<state>> Q;
    vector<int> prev(n + m, -1);
    Q.push({0, start, start});
    int dist;
    while (!Q.empty()){
        auto [d, v, from] = Q.top();
        Q.pop();
        if (prev[v] != -1)
            continue;
        prev[v] = from;
        if (v == goal){
            dist = d;
            break;
        }
        for (int to = 0; to < n + m; to++)
            Q.push({d + D[v][to], to, v});
    }

    vector<int> route;
    int cur = goal;
    while (cur != start){
        route.push_back(cur);
        cur = prev[cur];
    }
    route.push_back(start);
    reverse(route.begin(), route.end());
    
    return {dist, route};
}

vector<pair<int, int>> make_ans(vector<int> tour, vector<int64_t> vx, vector<int64_t> vy){
    vector<vector<int64_t>> D(n + m, vector<int64_t>(n + m));
    for (int i = 0; i < n; i++){
        for (int j = 0; j < n; j++)
            D[j][i] = D[i][j] = ((A[i] - A[j]) * (A[i] - A[j]) + (B[i] - B[j]) * (B[i] - B[j])) * alpha * alpha;
        for (int j = 0; j < m; j++)
            D[n + j][i] = D[i][n + j] = ((A[i] - vx[j]) * (A[i] - vx[j]) + (B[i] - vy[j]) * (B[i] - vy[j])) * alpha;
    }

    for (int i = 0; i < m; i++){
        for (int j = 0; j < m; j++)
            D[n + j][n + i] = D[n + i][n + j] = (vx[i] - vx[j]) * (vx[i] - vx[j]) + (vy[i] - vy[j]) * (vy[i] - vy[j]);
    }

    vector<pair<int, int>> res;
    int64_t d = 0;
    for (int i = 0; i < n; i++){
        auto [dist, route] = find_best_route(tour[i], tour[(i + 1) % n], D);
        for (const auto &r: route){
            const int s = r >= n ? 2 : 1;
            const int to = r >= n ? r - n : r;
            if (res.size() > 0 and res.back() == make_pair(s, to))
                continue;
            res.push_back({s, to});
        }
        d += dist;
    }

    // cerr << "d: " << d << endl;
    return res;
}

int main(){
    cin >> n >> m;
    A.resize(n);
    B.resize(n);
    for (int i = 0; i < n; i++)
        cin >> A[i] >> B[i];

    Time time;
    int64_t mi = 1e18;
    auto [mi_vx, mi_vy] = k_means(A, B, m);
    vector<int> mi_tour;
    constexpr int kinbo = 25;
    std::random_device seed_gen;
    std::default_random_engine engine(seed_gen());
    // 平均0.0、標準偏差1.0で分布させる
    constexpr double sigma = 50;
    std::normal_distribution<> dist(0.0, sigma);
    
    while (time.elapsed() < 960){
        auto vx = mi_vx;
        auto vy = mi_vy;
        for (int i = 0; i < m; i++){
            vx[i] += dist(engine) - sigma / 2;
            vy[i] += dist(engine) - sigma / 2;
            vx[i] = min<int64_t>(1000, max<int64_t>(0, vx[i]));
            vy[i] = min<int64_t>(1000, max<int64_t>(0, vy[i]));
        }

        auto D = calc_dist(vx, vy);

        TSP tsp(D, int64_t(1e15));
        auto tour = tsp.solve(1);
        const auto d = tsp.calc_score(tour);
        if (d < mi){
            mi = d;
            mi_vx = vx;
            mi_vy = vy;
            mi_tour = tour;
        }
    }

    auto vx = mi_vx;
    auto vy = mi_vy;
    auto D = calc_dist(vx, vy);
    TSP tsp(D, int64_t(1e15));
    auto tour = tsp.solve(20);
    const auto d = tsp.calc_score(tour);
    cerr << round(1e9/(1000 + sqrt(mi))) << endl;
    cerr << round(1e9/(1000 + sqrt(d))) << endl;
    mi = d;


    for (int i = 0; i < n; i++){
        if (tour[i] == 0)
            rotate(tour.begin(), tour.begin() + i, tour.end());
    }
    assert(tour[0] == 0);

    for (int i = 0; i < m; i++)
        cout << vx[i] << " " << vy[i] << endl;

    auto ans = make_ans(tour, vx, vy);
    cout << ans.size() << endl;
    for (const auto &[a, b]: ans)
        cout << a << " " << b + 1 << endl;

}
yukicoder

結果

テストケース

ソースコード
