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

• 問題: No.5007 Steiner Space Travel
• ユーザー: wanui
• 提出日時: 2022-07-30 19:30:37
• 言語: C++17(gcc12) (gcc 12.3.0 + boost 1.87.0)
• 結果: AC
• 実行時間: 928 ms / 1,000 ms
• コード長: 11,884 bytes
• コンパイル時間: 2,314 ms
• 実行使用メモリ: 6,952 KB
• スコア: 8,793,690
• 最終ジャッジ日時: 2022-07-30 19:31:10
• 合計ジャッジ時間: 32,333 ms
• ジャッジサーバーID （参考情報）: judge16 / judge10

ファイルパターン	結果
other	AC * 30

ソースコード

#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 のように移動
    vector<pii> ops;
    {
        int sid = pre_assign[0];
        int sid_init = sid;
        for (auto sid : cluster_order) {
            for (int pid = 0; pid < N; pid++) {
                if (pre_assign[pid] == sid) {
                    ops.push_back({PLANET, pid});
                    ops.push_back({STATION, sid});
                }
            }
            ops.push_back({STATION, sid});
        }
        ops.push_back({STATION, sid_init});
        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};
}
bool valid(point p) {
    return (0 <= p.i && p.i <= 1000 && 0 <= p.j && p.j <= 1000);
}
pair<vector<pii>, vector<point>> updates(vector<pii> init_ops, vector<point> init_stations) {
    // double start_time = now();
    auto stations = init_stations;
    auto ops = init_ops;

    const int OIT = 80;
    double begin_time = now();
    for (int oiter = 0; true; oiter++) {
        // stationの位置を微調整
        {
            ll old_score = -calc_cost(ops, stations);
            for (int iter = 0; iter < 1000; iter++) {
                int m = engine() % M;
                int dx = randint(-10, 10);
                int dy = randint(-10, 10);
                stations[m].i += dx;
                stations[m].j += dy;
                ll new_score = valid(stations[m]) ? (-calc_cost(ops, stations)) : (-1e18);
                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;
                }
            }
        }
        if (oiter == OIT) {
            break;
        }
        // stationの位置が完全に定まっているとする
        // 2つのplanet間の距離を定めることができる (経由する最大station数を2とする)
        // これはもうただのTSP
        vector<vector<int>> dists(N, vector<int>(N));
        {
            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;
                }
            }
        }

        vector<int> line;
        {
            // 初期解は元々のops
            for (auto op : ops) {
                if (op.first == PLANET) {
                    line.push_back(op.second);
                }
            }

            double et = END_TIME * (oiter + 1) / OIT;
            double end_time = END_TIME;
            double begin_temp = 10.0;
            double end_temp = 0.5;
            ll old_score = calc_tsp_score(line, dists);
            for (int iter = 0; now() < et; iter++) {
                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);

                // 2-opt
                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);
                }
            }
        }

        // 復元
        {
            ops.clear();
            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});
        }
    }

    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;
}