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

• 問題: No.5007 Steiner Space Travel
• ユーザー: wanui
• 提出日時: 2022-07-30 16:26:20
• 言語: C++17 (gcc 12.3.0 + boost 1.83.0)
• 結果: AC
• 実行時間: 934 ms / 1,000 ms
• コード長: 9,821 bytes
• コンパイル時間: 2,731 ms
• 実行使用メモリ: 6,952 KB
• スコア: 8,191,062
• 最終ジャッジ日時: 2022-07-30 16:26:53
• 合計ジャッジ時間: 33,045 ms
• ジャッジサーバーID （参考情報）: judge11 / judge15

入力	結果	実行時間; 実行使用メモリ
testcase_00	AC	907 ms; 6,948 KB
testcase_01	AC	933 ms; 6,952 KB
testcase_02	AC	906 ms; 4,900 KB
testcase_03	AC	914 ms; 4,900 KB
testcase_04	AC	928 ms; 4,908 KB
testcase_05	AC	905 ms; 4,900 KB
testcase_06	AC	934 ms; 6,952 KB
testcase_07	AC	908 ms; 4,904 KB
testcase_08	AC	926 ms; 4,900 KB
testcase_09	AC	919 ms; 4,904 KB
testcase_10	AC	906 ms; 6,952 KB
testcase_11	AC	934 ms; 4,904 KB
testcase_12	AC	906 ms; 4,900 KB
testcase_13	AC	934 ms; 4,908 KB
testcase_14	AC	907 ms; 4,900 KB
testcase_15	AC	907 ms; 4,904 KB
testcase_16	AC	932 ms; 4,900 KB
testcase_17	AC	906 ms; 4,904 KB
testcase_18	AC	931 ms; 4,904 KB
testcase_19	AC	906 ms; 6,952 KB
testcase_20	AC	907 ms; 4,904 KB
testcase_21	AC	932 ms; 4,900 KB
testcase_22	AC	907 ms; 6,952 KB
testcase_23	AC	932 ms; 6,948 KB
testcase_24	AC	906 ms; 4,900 KB
testcase_25	AC	921 ms; 6,952 KB
testcase_26	AC	916 ms; 4,908 KB
testcase_27	AC	913 ms; 6,952 KB
testcase_28	AC	933 ms; 4,900 KB
testcase_29	AC	906 ms; 4,904 KB

ソースコード

#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};
}
pair<vector<pii>, vector<point>> anneal(vector<pii> init_ops, vector<point> init_stations) {
    // double start_time = now();
    auto stations = init_stations;
    auto ops = init_ops;
    ll old_score = -calc_cost(ops, stations);
    int iter = 0;
    for (iter = 0; now() < END_TIME; 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 = -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;
        }
    }
    perr("anneal_iter=", iter);
    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) = anneal(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;
}