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

提出ソース
結果

問題	No.5007 Steiner Space Travel
ユーザー	wanui
提出日時	2022-07-30 17:57:47
言語	C++17 (gcc 12.3.0 + boost 1.83.0)
結果	WA
実行時間	-
コード長	15,273 bytes
コンパイル時間	2,598 ms
実行使用メモリ	6,952 KB
スコア	7,689,220
最終ジャッジ日時	2022-07-30 17:58:52
合計ジャッジ時間	33,628 ms
ジャッジサーバーID （参考情報）	judge16 / judge9
このコードへのチャレンジ
（要ログイン）
テストケース

テストケース表示
入力	結果	実行時間実行使用メモリ
testcase_00	AC	940 ms 6,948 KB
testcase_01	AC	920 ms 4,908 KB
testcase_02	WA	-
testcase_03	AC	968 ms 6,952 KB
testcase_04	AC	920 ms 4,900 KB
testcase_05	AC	921 ms 6,948 KB
testcase_06	AC	939 ms 6,948 KB
testcase_07	AC	916 ms 6,948 KB
testcase_08	AC	920 ms 6,952 KB
testcase_09	AC	949 ms 4,900 KB
testcase_10	AC	938 ms 4,904 KB
testcase_11	AC	920 ms 4,908 KB
testcase_12	AC	939 ms 6,952 KB
testcase_13	AC	917 ms 4,900 KB
testcase_14	AC	934 ms 4,900 KB
testcase_15	AC	919 ms 6,952 KB
testcase_16	AC	917 ms 4,904 KB
testcase_17	AC	940 ms 4,904 KB
testcase_18	AC	920 ms 4,900 KB
testcase_19	AC	918 ms 4,904 KB
testcase_20	AC	919 ms 6,952 KB
testcase_21	AC	931 ms 4,900 KB
testcase_22	AC	920 ms 4,900 KB
testcase_23	AC	941 ms 4,904 KB
testcase_24	AC	918 ms 4,904 KB
testcase_25	AC	937 ms 6,952 KB
testcase_26	AC	921 ms 4,904 KB
testcase_27	AC	925 ms 6,948 KB
testcase_28	AC	916 ms 4,904 KB
testcase_29	AC	927 ms 4,900 KB
権限があれば一括ダウンロードができます
ソースコード

raw source code
#include <bits/stdc++.h>
// clang-format off
using namespace std; using ll=long long; using ull=unsigned long long; using pll=pair<ll,ll>; const ll INF=4e18;
void print0(){}; template<typename H,typename... T> void print0(H h,T... t){cout<<h;print0(t...);}
void print(){print0("\n");}; template<typename H,typename... T>void print(H h,T... t){print0(h);if(sizeof...(T)>0)print0(" ");print(t...);}
void perr0(){}; template<typename H,typename... T> void perr0(H h,T... t){cerr<<h;perr0(t...);}
void perr(){perr0("\n");}; template<typename H,typename... T>void perr(H h,T... t){perr0(h);if(sizeof...(T)>0)perr0(" ");perr(t...);}
void ioinit() { cout<<fixed<<setprecision(15); cerr<<fixed<<setprecision(6); ios_base::sync_with_stdio(0); cin.tie(0); }
// clang-format on

// K-meansやりたい！
using pii = pair<int, int>;

double END_TIME = 900;
mt19937 engine(0);
clock_t start_time;
double now() {
    return 1000.0 * (clock() - start_time) / CLOCKS_PER_SEC;
}
void marathon_init() {
    start_time = clock();
    random_device seed_gen;
    engine.seed(seed_gen());
}
int randint(int mn, int mx) {
    int rng = mx - mn + 1;
    return mn + (engine() % rng);
}
double uniform(double x, double y) {
    const int RND = 1e8;
    double mean = (x + y) / 2.0;
    double dif = y - mean;
    double p = double(engine() % RND) / RND;
    return mean + dif * (1.0 - 2.0 * p);
}
bool anneal_accept(double new_score, double old_score, double cur_time, double begin_time, double end_time, double begin_temp, double end_temp) {
    const int ANNEAL_RND = 1e8;
    const double ANNEAL_EPS = 1e-6;
    double temp = cur_time * (end_temp - begin_temp) / (end_time - begin_time) + (end_time * begin_temp - end_temp * begin_time) / (end_time - begin_time);
    return (exp((new_score - old_score) / temp) > double(engine() % ANNEAL_RND) / ANNEAL_RND + ANNEAL_EPS);
}

struct point {
    int i;
    int j;
};

const int PLANET = 1;
const int STATION = 2;
const int N = 100;
const int M = 8;
const int ALPHA = 5;

vector<point> planets(N);

int distance2(point a, point b) {
    return (a.i - b.i) * (a.i - b.i) + (a.j - b.j) * (a.j - b.j);
}
ll calc_cost(vector<pii> ops, vector<point> stations) {
    // round(10^9/(1000+sqrt(S)))
    ll S = 0;
    int m = ops.size();
    for (int i = 1; i < m; i++) {
        pii pre = ops[i - 1];
        pii cur = ops[i];
        int ratio = 1;
        point pp;
        point cp;
        if (pre.first == STATION) {
            pp = stations[pre.second];
        } else {
            pp = planets[pre.second];
            ratio *= ALPHA;
        }
        if (cur.first == STATION) {
            cp = stations[cur.second];
        } else {
            cp = planets[cur.second];
            ratio *= ALPHA;
        }
        S += ratio * distance2(cp, pp);
    }
    return S;
}
pair<vector<pii>, vector<point>> kmeans() {
    // k-means

    // 初期値
    vector<int> perm(N);
    iota(perm.begin(), perm.end(), 0);
    shuffle(perm.begin(), perm.end(), engine);
    vector<point> stations;
    for (int i = 0; i < M; i++) {
        stations.push_back(planets[perm[i]]);
    }

    vector<int> pre_assign(N, -1);
    for (int iter = 0; iter < 10000; iter++) {
        vector<int> assign(N);
        for (int pid = 0; pid < N; pid++) {
            auto pl = planets[pid];
            int mindist = 1e9;
            int minstation = -1;
            for (int sid = 0; sid < M; sid++) {
                auto st = stations[sid];
                int d = distance2(pl, st);
                if (mindist > d) {
                    mindist = d;
                    minstation = sid;
                }
            }
            assign[pid] = minstation;
        }
        if (pre_assign == assign) {
            break;
        }
        vector<point> sums(M, {0, 0});
        vector<int> counts(M);
        for (int pid = 0; pid < N; pid++) {
            auto pl = planets[pid];
            auto sid = assign[pid];
            sums[sid].i += pl.i;
            sums[sid].j += pl.j;
            counts[sid]++;
        }
        for (int sid = 0; sid < M; sid++) {
            int cnt = counts[sid];
            auto su = sums[sid];
            if (cnt == 0) {
                // cnt==0 はないと思うが... stationの位置に変更なし
                continue;
            }

            stations[sid] = {su.i / cnt, su.j / cnt};
        }
        pre_assign = assign;
    }
    vector<int> cluster_order(M);
    // クラスタをたどる順を全探索で探す(bitdpの方が速い)
    {
        int sid_init = pre_assign[0];
        vector<int> perm(M);
        iota(perm.begin(), perm.end(), 0);

        int mindist = 1e9;
        do {
            if (perm[0] != sid_init) continue;
            int s = 0;
            for (ll i = 0; i < M; i++) {
                s += distance2(stations[perm[i]], stations[perm[(i + 1) % M]]);
            }
            if (mindist > s) {
                mindist = s;
                cluster_order = perm;
            }
        } while (next_permutation(perm.begin(), perm.end()));
    }
    // station1->planet1.1->station1->planet1.2->station1->station2 のように移動
    // だが、planet1.1->planet1.2 のようにやるほうが短いならそれを選んでもよい
    vector<pii> ops;
    {
        ops.push_back({PLANET, 0});
        for (auto sid : cluster_order) {
            ops.push_back({STATION, sid});

            // 座圧？
            vector<int> i2p;
            map<int, int> p2i;
            for (int pid = 1; pid < N; pid++) {  // pid==0 は特別扱い
                if (pre_assign[pid] == sid) {
                    p2i[pid] = i2p.size();
                    i2p.push_back(pid);
                }
            }
            if (i2p.size() > 0) {
                // 距離に偏りのあるTSP. まずは貪欲 / TODO 焼きなましとかしたほうがよさそう
                int m = i2p.size();
                point st = stations[sid];
                vector<bool> done(m);
                int donenum = 0;
                while (true) {
                    pii cur = ops.back();
                    if (cur.first == STATION) {
                        // 適当に選ぶ
                        for (int i = 0; i < m; i++) {
                            if (!done[i]) {
                                ops.push_back({PLANET, i2p[i]});
                                break;
                            }
                        }
                        continue;
                    }
                    int pid = cur.second;
                    int i = p2i[pid];
                    done[i] = true;
                    donenum++;
                    if (donenum == m) {
                        ops.push_back({STATION, sid});
                        break;
                    }
                    int d_station = ALPHA * distance2(planets[pid], st);
                    int mind = d_station;
                    int minj = -1;
                    for (int j = 0; j < m; j++) {
                        if (done[j]) continue;
                        int dj = ALPHA * ALPHA * distance2(planets[pid], planets[i2p[j]]);
                        if (mind > dj) {
                            mind = dj;
                            minj = j;
                        }
                    }
                    if (minj == -1) {
                        ops.push_back({STATION, sid});
                        continue;
                    } else {
                        ops.push_back({PLANET, i2p[minj]});
                        continue;
                    }
                }
            }
        }
        ops.push_back({STATION, cluster_order[0]});
        ops.push_back({PLANET, 0});
    }
    return {ops, stations};
}
ll calc_tsp_score(vector<int>& line, vector<vector<int>>& dists) {
    int m = line.size();
    ll d = 0;
    for (int i = 0; i < m - 1; i++) {
        d += dists[line[i]][line[i + 1]];
    }
    return -d;
}
pair<int, tuple<pii, pii, pii>> get_mindist(point pa, point pb, int paid, vector<point>& stations) {
    // 2つのplanet間の最短距離 (経由する最大station数を2とする)

    int mindist = ALPHA * ALPHA * distance2(pa, pb);  // 直接行く
    tuple<pii, pii, pii> curops = {{PLANET, paid}, {-1, -1}, {-1, -1}};
    for (int sid = 0; sid < M; sid++) {
        auto ss = stations[sid];
        int d = ALPHA * (distance2(pa, ss) + distance2(ss, pb));  // station s 経由
        if (mindist > d) {
            mindist = d;
            curops = {{PLANET, paid}, {STATION, sid}, {-1, -1}};
        }
        for (int tid = 0; tid < M; tid++) {
            auto st = stations[tid];
            int d = ALPHA * distance2(pa, ss) + distance2(ss, st) + ALPHA * distance2(st, pb);  // station s-t 経由
            if (mindist > d) {
                mindist = d;
                curops = {{PLANET, paid}, {STATION, sid}, {STATION, tid}};
            }
        }
    }
    return {mindist, curops};
}
void setdists(vector<vector<int>>& dists, vector<point>& stations) {
    for (int a = 0; a < N; a++) {
        dists[a][a] = 0;
        for (int b = a + 1; b < N; b++) {
            int mindist = get_mindist(planets[a], planets[b], a, stations).first;
            dists[a][b] = mindist;
            dists[b][a] = mindist;
        }
    }
}
void update_dists(vector<vector<int>>& dists, vector<point>& stations, int tgt_sid) {
    for (int a = 0; a < N; a++) {
        dists[a][a] = 0;
        for (int b = a + 1; b < N; b++) {
            int mindist = dists[a][b];
            auto pa = planets[a];
            auto pb = planets[b];
            for (int sid = 0; sid < M; sid++) {
                auto ss = stations[sid];
                if (sid == tgt_sid) {
                    for (int tid = 0; tid < M; tid++) {
                        auto st = stations[tid];
                        int d = ALPHA * distance2(pa, ss) + distance2(ss, st) + ALPHA * distance2(st, pb);  // station s-t 経由
                        mindist = min(mindist, d);
                    }
                } else {
                    int tid = tgt_sid;
                    auto st = stations[tid];
                    int d = ALPHA * distance2(pa, ss) + distance2(ss, st) + ALPHA * distance2(st, pb);  // station s-t 経由
                    mindist = min(mindist, d);
                }
            }
            dists[a][b] = mindist;
            dists[b][a] = mindist;
        }
    }
}
pair<vector<pii>, vector<point>> updates(vector<pii> init_ops, vector<point> init_stations) {
    // double start_time = now();
    auto stations = init_stations;

    double begin_time = now();
    double end_time = END_TIME;
    double begin_temp = 10.0;
    double end_temp = 0.5;

    // stationの位置が完全に定まっているとする
    // 2つのplanet間の距離を定めることができる (経由する最大station数を2とする)
    // これはもうただのTSP

    // 初期解はk-meansで求めたops,stations
    vector<int> line;
    for (auto op : init_ops) {
        if (op.first == PLANET) {
            line.push_back(op.second);
        }
    }

    vector<vector<int>> dists(N, vector<int>(N));
    setdists(dists, stations);

    ll old_score = calc_tsp_score(line, dists);
    int iter;
    for (iter = 1; now() < end_time; iter++) {
        if (iter % 200 != 0) {
            // 2-opt
            int x = 1 + engine() % (int(line.size()) - 2);
            int y = 1 + engine() % (int(line.size()) - 2);
            if (x == y) continue;
            if (x > y) swap(x, y);

            vector<int> newline = line;
            for (int i = x; i <= y; i++) {
                int j = -i + y + x;
                if (i >= j) break;
                swap(newline[i], newline[j]);
            }
            ll new_score = calc_tsp_score(newline, dists);
            if (anneal_accept(new_score, old_score, now(), begin_time, end_time, begin_temp, end_temp)) {
                old_score = new_score;
                swap(line, newline);
            }
        } else {
            // stationの位置を微調整
            int m = engine() % M;
            int dx = randint(-10, 10);
            int dy = randint(-10, 10);

            auto newstations = stations;
            newstations[m].i += dx;
            newstations[m].j += dy;

            update_dists(dists, newstations, m);  // TODO 差分更新
            ll new_score = calc_tsp_score(line, dists);
            if (anneal_accept(new_score, old_score, now(), begin_time, end_time, begin_temp, end_temp)) {
                // if (new_score >= old_score) {
                old_score = new_score;
                swap(stations, newstations);
            } else {
                update_dists(dists, stations, m);
            }
        }
    }
    perr("anneal_iter=", iter);

    // 復元
    vector<pii> ops;
    {
        int m = line.size();
        for (int i = 0; i < m - 1; i++) {
            int paid = line[i];
            auto pa = planets[line[i]];
            auto pb = planets[line[i + 1]];

            pii op1, op2, op3;
            tie(op1, op2, op3) = get_mindist(pa, pb, paid, stations).second;

            if (op1.first >= 0) ops.push_back(op1);
            if (op2.first >= 0) ops.push_back(op2);
            if (op3.first >= 0) ops.push_back(op3);
        }
        ops.push_back({PLANET, 0});
    }

    // stationを微調整
    {
        ll old_score = -calc_cost(ops, stations);
        for (int iter = 0; iter < 20000; iter++) {
            int m = engine() % M;
            int dx = randint(-5, 5);
            int dy = randint(-5, 5);
            stations[m].i += dx;
            stations[m].j += dy;
            ll new_score = -calc_cost(ops, stations);
            if (new_score > old_score) {
                old_score = new_score;
                // perr("update!", iter, -old_score, -new_score);
            } else {
                stations[m].i -= dx;
                stations[m].j -= dy;
            }
        }
    }
    return {ops, stations};
}

void solve() {
    // kmeansを何度か回して、よさげな解を見つける
    vector<pii> best_km_ops;
    vector<point> best_km_stations;
    int best_km_cost = 1e9;
    while (now() < END_TIME * 0.2) {
        vector<pii> ops;
        vector<point> stations;
        tie(ops, stations) = kmeans();
        ll cost = calc_cost(ops, stations);
        if (cost < best_km_cost) {
            best_km_ops = ops;
            best_km_stations = stations;
            best_km_cost = cost;
        }
    }

    // k-meansを初期解にして微調整
    vector<pii> best_ops;
    vector<point> best_stations;
    tie(best_ops, best_stations) = updates(best_km_ops, best_km_stations);

    // 出力
    for (auto st : best_stations) {
        print(st.i, st.j);
    }
    print(best_ops.size());
    for (auto op : best_ops) {
        print(op.first, op.second + 1);
    }

    // デバッグ用
    int sc = int(0.5 + 1e9 / (1000.0 + sqrt(calc_cost(best_ops, best_stations))));
    perr("score=", sc);
    perr("time=", int(now()));
}
int main() {
    marathon_init();
    ioinit();
    int n, m;
    cin >> n >> m;
    for (int i = 0; i < n; i++) {
        int a, b;
        cin >> a >> b;
        planets[i] = {i : a, j : b};
    }
    solve();
    return 0;
}
yukicoder

結果

テストケース

ソースコード
