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

• 問題: No.5007 Steiner Space Travel
• ユーザー: eijirou
• 提出日時: 2022-07-30 17:56:51
• 言語: C++17(gcc12) (gcc 12.3.0 + boost 1.87.0)
• 結果: AC
• 実行時間: 999 ms / 1,000 ms
• コード長: 11,049 bytes
• コンパイル時間: 4,046 ms
• 実行使用メモリ: 4,136 KB
• スコア: 8,344,452
• 最終ジャッジ日時: 2022-07-30 17:57:28
• 合計ジャッジ時間: 36,690 ms
• ジャッジサーバーID （参考情報）: judge12 / judge13

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

ソースコード

#pragma GCC target("avx2")
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")

#include <bits/stdc++.h>
//#include <atcoder/all>
using namespace std;
//using namespace atcoder;

template <class T>
T sum(vector<T>& array) {
    return accumulate(array.begin(), array.end(), (T)0);
}

template <class T>
int index(vector<T>& data, T element) {
    return distance(data.begin(), find(data.begin(),data.end(),element));
}

template <class T>
int contain(vector<T>& data, T element) {
    return !(index(data, element) == (int)data.size());
}

struct Randomizer {
    mt19937_64 engine;
    uniform_int_distribution<> dist = uniform_int_distribution<>(0, INT_MAX);

    Randomizer() {
        random_device seed_gen;
        engine = mt19937_64(seed_gen());
    }

    Randomizer(int seed) {
        engine = mt19937_64(seed);
    }

    double uniform(double l, double r) {
        assert(l <= r);
        return l + (r-l)*dist(engine)/INT_MAX;
    }

    int randrange(int l, int r) {
        assert(l < r);
        return l + dist(engine)%(r-l);
    }

    template <class T>
    T choice(vector<T>& data) {
        return data.at(randrange(0, data.size()));
    }

    vector<int> sample(int l, int r, int num) {
        assert(0 < num && num <= r - l);
        vector<int> ret;
        while ((int)ret.size() < num) {
            int x = randrange(l, r);
            if (!contain(ret, x)) {
                ret.push_back(x);
            }
        }
        sort(ret.begin(), ret.end());
        return ret;
    }
};

struct Timer {
    timespec start;
    Randomizer randomizer;

    Timer() {}

    Timer(int seed) {
        randomizer = Randomizer(seed);
    }

    void begin() {
        clock_gettime(CLOCK_REALTIME, &start);
    }

    double stopwatch() {
        timespec end;
        clock_gettime(CLOCK_REALTIME, &end);
        double sec = end.tv_sec - start.tv_sec;
        double nsec = end.tv_nsec - start.tv_nsec;
        return sec + nsec/1000000000;
    }

    bool scheduler(double delta, double time_limit, double t0, double t1) {
        assert(0.0 < time_limit && 0.0 <= t1 && t1 <= t0);
        if (0.0 <= delta) {
            return true;
        } else {
            double ratio = stopwatch() / time_limit;
            double t = pow(t0, 1.0-ratio) * pow(t1, ratio);
            return randomizer.uniform(0.0, 1.0) < pow(2.0, delta/t);
        }
    }
};

vector<int> tsp_solver(vector<vector<int>> cost_matrix) {
    // n を始点および終点とするTSPの厳密解を求める.
    int n = cost_matrix.size() - 1;
    int s = n;
    int inf = 1 << 30;
 
    // dp
    vector<vector<int>> dp(n, vector<int>(1<<n,inf));
    for (int i = 0; i < n; i++) {
        dp[i][1<<i] = cost_matrix[s][i];
    }
    for (int mask = 1; mask < 1<<n; mask++) {
        for (int i = 0; i < n; i++) {
            if ((mask>>i & 1) == 0) {
                continue;
            }
            for (int j = 0; j < n; j++) {
                if (mask>>j & 1) {
                    continue;
                }
                dp[j][mask^(1<<j)] = min(dp[j][mask^(1<<j)], dp[i][mask]+cost_matrix[i][j]);
            }
        }
    }
    // 経路復元
    int last = 0;
    int cost = inf;
    for (int i = 0; i < n; i++) {
        if (dp[i].back() + cost_matrix[i][s] < cost) {
            last = i;
            cost = dp[i].back() + cost_matrix[i][s];
        }
    }
    vector<int> ret = {last};
    int mask = (1<<n) - 1;
    mask ^= 1 << last;
    while (mask) {
        int j = ret.back();
        for (int i = 0; i < n; i++) {
            if (mask>>i & 1 && dp[i][mask] + cost_matrix[i][j] == dp[j][mask^1<<j]) {
                ret.push_back(i);
                mask ^= 1 << i;
                break;
            }
        }
        assert(ret.back() != j);
    }
    ret.push_back(s);
    reverse(ret.begin(), ret.end());
    return ret;
}

const int alpha = 5;
const double time_limit = 0.997;
const int inf = 1 << 30;

struct Answer {
    int score;
    vector<int> c, d, t, r;

    void print() {
        for (int i = 0; i < (int)c.size(); i++) {
            cout << c[i] << ' ' << d[i] << '\n';
        }
        cout << t.size() << '\n';
        for (int i = 0; i < (int)t.size(); i++) {
            cout << t[i] << ' ' << r[i] + 1 << '\n';
        }
    }
};

struct Solver {
    Timer timer;
    Randomizer randomizer;
    int n, m;
    vector<int> a, b;
    vector<vector<int>> groups;
    vector<int> c, d;

    Solver() {
        timer.begin();
    }

    void input() {
        cin >> n >> m;
        a = vector<int>(n);
        b = vector<int>(n);
        for (int i = 0; i < n; i++) {
            cin >> a[i] >> b[i];
        }
    }

    pair<int,int> get_center(vector<int>& group) {
        if (group.empty()) {
            int j = randomizer.randrange(0, n);
            return {a[j], b[j]};
        }
        int x = 0;
        int y = 0;
        for (int i : group) {
            x += a[i];
            y += b[i];
        }
        return {x/(int)group.size(), y/(int)group.size()};
    }

    void k_means() {
        groups = vector<vector<int>>(m);
        for (int i = 0; i < m; i++) {
            groups[i].push_back(i);
        }
        for (int i = m; i < n; i++) {
            groups[randomizer.randrange(0,m)].push_back(i);
        }
        while (true) {
            vector<pair<int,int>> centers(m);
            for (int i = 0; i < m; i++) {
                centers[i] = get_center(groups[i]);
            }
            vector<vector<int>> new_groups(m);
            for (int i = 0; i < n; i++) {
                int best = 0;
                int min_cost = inf;
                for (int j = 0; j < m; j++) {
                    auto [x, y] = centers[j];
                    int cost = (a[i]-x)*(a[i]-x) + (b[i]-y)*(b[i]-y);
                    if (cost < min_cost) {
                        best = j;
                        min_cost = cost;
                    }
                }
                new_groups[best].push_back(i);
            }
            if (groups == new_groups) {
                break;
            }
            groups = new_groups;
        }
    }

    Answer make() {
        k_means();
        c = vector<int>(m);
        d = vector<int>(m);
        for (int i = 0; i < m; i++) {
            pair<int,int> p = get_center(groups[i]);
            c[i] = p.first;
            d[i] = p.second;
        }
        vector<vector<int>> cost_matrix(m, vector<int>(m));
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < m; j++) {
                int dist = (c[i]-c[j])*(c[i]-c[j]) + (d[i]-d[j])*(d[i]-d[j]);
                cost_matrix[i][j] = dist;
                cost_matrix[j][i] = dist;
            }
        }
        vector<int> order = tsp_solver(cost_matrix);
        int start = -1;
        for (int i = 0; i < m; i++) {
            if (contain(groups[order[i]], 0)) {
                start = i;
            }
        }
        vector<int> new_order(m);
        for (int i = 0; i < m; i++) {
            new_order[i] = order[(start+i) % m];
        }
        order = new_order;

        for (int i = 0; i < m; i++) {
            vector<pair<int,pair<int,int>>> data;
            for (int j : groups[i]) {
                data.push_back({j, {a[j]-c[i], b[j]-d[i]}});
            }
            sort(data.begin(), data.end(), [](const auto &a, const auto &b) {
                return atan2(a.second.second, a.second.first) < atan2(b.second.second, b.second.first);
            });
            int start = 0;
            if (contain(groups[i], 0)) {
                for (int j = 0; j < (int)data.size(); j++) {
                    if (data[j].first == 0) {
                        start = j;
                        break;
                    }
                }
            }
            for (int j = 0; j < (int)data.size(); j++) {
                groups[i][j] = data[(start+j) % (int)data.size()].first;
            }
        }
        vector<int> t = {1};
        vector<int> r = {0};
        for (int i = 0; i < m; i++) {
            int station = order[i];
            for (int star : groups[station]) {
                if (star == 0) {
                    continue;
                }
                if (t.back() == 1) {
                    int e1 = alpha * alpha * ((a[star]-a[r.back()])*(a[star]-a[r.back()]) + (b[star]-b[r.back()])*(b[star]-b[r.back()]));
                    int e2 = alpha * ((a[star]-c[station])*(a[star]-c[station]) + (b[star]-d[station])*(b[star]-d[station]));
                    e2 += alpha * ((a[r.back()]-c[station])*(a[r.back()]-c[station]) + (b[r.back()]-d[station])*(b[r.back()]-d[station]));
                    if (e2 < e1) {
                        t.push_back(2);
                        r.push_back(station);
                    }
                }
                t.push_back(1);
                r.push_back(star);
            }
            if (t.back() == 1) {
                t.push_back(2);
                r.push_back(station);
            }
            int next_station = order[(i+1) % m];
            t.push_back(2);
            r.push_back(next_station);
        }
        t.push_back(1);
        r.push_back(0);

        vector<int> sum_x(m, 0);
        vector<int> sum_y(m, 0);
        vector<int> counts(m, 0);
        for (int i = 0; i < (int)t.size() - 1; i++) {
            int station = -1;
            int star = -1;
            if (t[i] == 1 && t[i+1] == 2) {
                station = r[i+1];
                star = r[i];
            } else if (t[i] == 2 && t[i+1] == 1) {
                station = r[i];
                star = r[i+1];
            } else {
                continue;
            }
            sum_x[station] += a[star];
            sum_y[station] += b[star];
            counts[station] += 1;
        }
        for (int i = 0; i < m; i++) {
            c[i] = sum_x[i] / counts[i];
            d[i] = sum_y[i] / counts[i];
        }
        int score = 0;
        for (int i = 0; i < (int)t.size() - 1; i++) {
            int sx, sy, tx, ty;
            if (t[i] == 1) {
                sx = a[r[i]];
                sy = b[r[i]];
            } else {
                sx = c[r[i]];
                sy = d[r[i]];
            }
            if (t[i+1] == 1) {
                tx = a[r[i+1]];
                ty = b[r[i+1]];
            } else {
                tx = c[r[i+1]];
                ty = d[r[i+1]];
            }
            int dist = (sx-tx)*(sx-tx) + (sy-ty)*(sy-ty);
            if (t[i] + t[i+1] == 2) {
                dist *= alpha * alpha;
            } else if (t[i] + t[i+1] == 3) {
                dist *= alpha;
            }
            score += dist;
        }
        return {score, c, d, t, r};
    }

    void solve() {
        Answer best = make();
        while (timer.stopwatch() < time_limit) {
            Answer ans = make();
            if (ans.score < best.score) {
                best = ans;
            }
        }
        best.print();
    }
};

int main() {
    Solver solver;
    solver.input();
    solver.solve();
    return 0;
}