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

• 問題: No.5007 Steiner Space Travel
• ユーザー: milanis48663220
• 提出日時: 2022-07-30 16:23:52
• 言語: C++17(gcc12) (gcc 12.3.0 + boost 1.87.0)
• 結果: AC
• 実行時間: 936 ms / 1,000 ms
• コード長: 10,206 bytes
• コンパイル時間: 1,808 ms
• 実行使用メモリ: 6,952 KB
• スコア: 8,481,660
• 最終ジャッジ日時: 2022-07-30 16:24:24
• 合計ジャッジ時間: 32,202 ms
• ジャッジサーバーID （参考情報）: judge11 / judge13

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

ソースコード

#include <iostream>
#include <algorithm>
#include <iomanip>
#include <vector>
#include <queue>
#include <deque>
#include <set>
#include <map>
#include <tuple>
#include <cmath>
#include <numeric>
#include <functional>
#include <cassert>
#include <random>

std::random_device rnd;
std::mt19937 mt(rnd());

int randint(const int l, const int r){
    return mt()%(r - l) + l;
}

#define debug_value(x) cerr << "line" << __LINE__ << ":<" << __func__ << ">:" << #x << "=" << x << endl;
#define debug(x) cerr << "line" << __LINE__ << ":<" << __func__ << ">:" << x << endl;

template<class T> inline bool chmax(T& a, T b) { if (a < b) { a = b; return 1; } return 0; }
template<class T> inline bool chmin(T& a, T b) { if (a > b) { a = b; return 1; } return 0; }

using namespace std;
typedef long long ll;

template<typename T>
vector<vector<T>> vec2d(int n, int m, T v){
    return vector<vector<T>>(n, vector<T>(m, v));
}

template<typename T>
vector<vector<vector<T>>> vec3d(int n, int m, int k, T v){
    return vector<vector<vector<T>>>(n, vector<vector<T>>(m, vector<T>(k, v)));
}

template<typename T>
void print_vector(vector<T> v, char delimiter=' '){
    if(v.empty()) {
        cout << endl;
        return;
    }
    for(int i = 0; i+1 < v.size(); i++) cout << v[i] << delimiter;
    cout << v.back() << endl;
}

const ll infll = 1e+18;

vector<vector<ll>> warshallfloyd(vector<vector<ll>> g){
    int n = g.size();
    vector<vector<ll>> dist(n, vector<ll>(n));
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            dist[i][j] = g[i][j];
        }
        dist[i][i] = 0;
    }
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            for(int k = 0; k < n; k++){
                ll new_len = dist[j][i] + dist[i][k];
                dist[j][k] = min(new_len, dist[j][k]);
            }
        }
    }
    return dist;
}

const int n = 100;
const int m = 8;
int a[n], b[n];

using P = pair<int, int>;
using Pdd = pair<double, double>;

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

double dist(double x0, double y0, double x1, double y1){
    double dx = x0-x1;
    double dy = y0-y1;
    return sqrt(dx*dx+dy*dy);
}

vector<P> k_means(vector<P> points, int k=m, int iter=100){
    vector<ll> x_sum(k), y_sum(k), cnt(k);
    vector<int> perm(points.size());
    iota(perm.begin(), perm.end(), 0);
    shuffle(perm.begin(), perm.end(), mt);
    for(int i = 0; i < points.size(); i++){
        x_sum[i%k] += points[i].first;
        y_sum[i%k] += points[i].second;
        cnt[i%k]++;
    }
    vector<Pdd> cur_ans;
    for(int i = 0; i < k; i++) cur_ans.push_back(P(x_sum[i]/cnt[i], y_sum[i]/cnt[i]));
    for(int i = 0; i < iter; i++){
        vector<vector<P>> clusters(k);
        for(auto [x, y]: points){
            int idx = -1;
            double min_dist = 1e9;
            for(int j = 0; j < k; j++){
                auto [xx, yy] = cur_ans[j];
                if(chmin(min_dist, dist(xx, yy, x, y))) idx = j;
            }
            clusters[idx].push_back(P(x, y));
        }
        for(int i = 0; i < k; i++){
            int sz = clusters[i].size();
            if(sz == 0) continue;
            int x_sum = 0, y_sum = 0;
            for(auto [x, y]: clusters[i]){
                x_sum += x;
                y_sum += y;
            }
            cur_ans[i] = P(double(x_sum)/double(sz), double(y_sum)/double(sz));
        }
    }
    vector<P> ans;
    for(int i = 0; i < k; i++){
        auto [xx, yy] = cur_ans[i];
        int x = max(0, (int)round(xx));
        int y = max(0, (int)round(yy));
        ans.push_back(P(x, y));
    }
    return ans;
}

bool seen_shortest_path[n+m];
vector<int> shortest_path(const vector<vector<ll>> &dist, const vector<vector<ll>> &g, const int s, const int t){
    int sz = dist.size();
    int cur = t;
    vector<int> ans;
    while(cur != s){
        ans.push_back(cur);
        seen_shortest_path[cur] = true;
        for(int i = 0; i < sz; i++){
            if(i == cur) continue;
            if(seen_shortest_path[i]) continue;
            if(dist[s][cur] == dist[s][i]+g[i][cur]){
                cur = i;
                break;
            }
            assert(i != sz-1);
        }
    }
    for(int v: ans) seen_shortest_path[v] = false;
    ans.push_back(s);
    reverse(ans.begin(), ans.end());
    return ans;
}

vector<int> greedy_tour(const vector<vector<ll>> &dist, int s){
    vector<bool> visited(n);
    int cur = s;
    vector<int> ans;
    for(int i = 0; i < n; i++){
        ans.push_back(cur);
        visited[cur] = true;
        ll min_dist = infll;
        int nx = -1;
        for(int j = 0; j < n; j++){
            if(visited[j]) continue;
            if(chmin(min_dist, dist[cur][j])){
                nx = j;
            }
        }
        cur = nx;
    }
    return ans;
}

class TspSolution{
    public:
    int sz;
    ll len;
    vector<int> tour;
    vector<vector<ll>> dist;
    ll tour_len(){
        ll len = 0;
        for(int i = 0; i < sz; i++){
            int j = (i+1)%sz;
            len += dist[tour[i]][tour[j]];
        }
        return len;
    }
    TspSolution(vector<int> tour, vector<vector<ll>> dist): sz(tour.size()), tour(tour), dist(dist){
        len = tour_len();
    }
     pair<int, int> suggest_two_opt(){
        int i = (mt()%(sz-1)) + 1;
        int j = (mt()%(sz-1)) + 1;
        while(i == j) j = (mt()%(m-1)) + 1;
        if(i > j) swap(i, j);
        return make_pair(i, j);
    }
    ll swapped_diff(int i, int j){
        int increase = dist[tour[i]][tour[j]] + dist[tour[(i+1)%sz]][tour[(j+1)%sz]];
        int decrease = dist[tour[i]][tour[(i+1)%sz]] + dist[tour[j]][tour[(j+1)%sz]];
        return increase-decrease;
    };
    void apply_swap(int i, int j){
        vector<int> new_tour;
        for(int k = 0; k <= i; k++) new_tour.push_back(tour[k]);
        for(int k = j; k >= i+1; k--) new_tour.push_back(tour[k]);
        for(int k = j+1; k < sz; k++) new_tour.push_back(tour[k]);
        tour = new_tour;
        len = tour_len();
    };
};

vector<int> solve_tsp(const vector<vector<ll>> &dist, int s, double time_limit, double start_temp = 50.0, double end_temp = 0.0){
    clock_t start = clock();
    auto sol = TspSolution(greedy_tour(dist, s), dist);

    // sa settings
    uniform_real_distribution<float> rand_dist(0.0, 1.0);
    auto accept_proba = [&](int diff, double temp){
        return min(1.0, exp((double)-diff/temp));
    };

    auto try_two_opt = [&](double temperature){
        auto [i, j] = sol.suggest_two_opt();
        int diff = sol.swapped_diff(i, j);
        if(rand_dist(mt) <= accept_proba(diff, temperature)) sol.apply_swap(i, j);
        // if(diff < 0) sol.apply_swap(i, j);
    };

    double remaining_time = time_limit - (double)(clock()-start)/CLOCKS_PER_SEC;
    double temperature = start_temp;
    for(int iter = 0;; iter++){
        if(iter%100 == 0){
            clock_t cur_time = clock();
            remaining_time = time_limit - (double)(cur_time-start)/CLOCKS_PER_SEC;
            temperature = end_temp + (start_temp-end_temp) * (remaining_time/time_limit);
        }
        if(remaining_time < 0.0) {
            break;
        }
        try_two_opt(temperature);
    }
    auto ans = sol.tour;
    ans.push_back(s);
    return ans;
}

class Solution{
    public:
    vector<P> stations;
    vector<int> tp;
    vector<int> idx;
    Solution(vector<P> stations, vector<int> tp, vector<int> idx): stations(stations), tp(tp), idx(idx) {}
    void output(){
        for(auto [x, y]: stations) cout << x << ' ' << y << endl;
        int v = tp.size();
        assert(idx.size() == v);
        cout << v << endl;
        for(int i = 0; i < v; i++) {
            cout << tp[i] << ' ' << idx[i]+1 << endl;
        }
    }
    ll eval(){
        int v = tp.size();
        auto get_x = [&](int tp, int idx){
            if(tp == 1) return a[idx];
            return stations[idx].first;
        };
        auto get_y = [&](int tp, int idx){
            if(tp == 1) return b[idx];
            return stations[idx].second;
        };
        ll s = 0;
        for(int i = 0; i < v-1; i++) {
            ll dx = get_x(tp[i], idx[i])-get_x(tp[i+1], idx[i+1]);
            ll dy = get_y(tp[i], idx[i])-get_y(tp[i+1], idx[i+1]);
            ll coef = (tp[i] == 1 ? 5 : 1)*(tp[i+1] == 1 ? 5 : 1);
            s += (dx*dx+dy*dy)*coef;
        }
        ll ans = round((double)1000000000/(1000+sqrt(s)));
        return ans;
    }
};

Solution solve(){
    auto g = vec2d(n+m, n+m, infll);
    auto coef = [](int i, int j){
        if(i < n && j < n) return 25;
        if(i < n || j < n) return 5;
        return 1;
    };
    vector<P> planets(n);
    for(int i = 0; i < n; i++) planets[i] = P(a[i], b[i]);
    auto stations = k_means(planets);
    auto get_x = [&](int i){
        if(i < n) return a[i];
        return stations[i-n].first;
    };
    auto get_y = [&](int i){
        if(i < n) return b[i];
        return stations[i-n].second;
    };
    for(int i = 0; i < n+m; i++){
        for(int j = i+1; j < n+m; j++){
            ll dx = get_x(i)-get_x(j), dy = get_y(i)-get_y(j);
            ll dist = (dx*dx+dy*dy)*coef(i, j);
            g[i][j] = dist;
            g[j][i] = dist;
        }
    }
    auto dist = warshallfloyd(g);
    auto ord = solve_tsp(dist, 0, 0.9);
    assert(ord.size() == n+1);
    auto sol = Solution(stations, vector<int>(), vector<int>());
    for(int i = 0; i < n; i++){
        int s = ord[i], t = ord[i+1];
        auto path = shortest_path(dist, g, s, t);
        auto add = [&](int idx){
            if(idx < n) {
                sol.tp.push_back(1);
                sol.idx.push_back(idx);
            }else{
                sol.tp.push_back(2);
                sol.idx.push_back(idx-n);
            }
        };
        if(i == 0){
            add(path[0]);
        }
        for(int j = 1; j < path.size(); j++){
            add(path[j]);
        }
    }
    return sol;
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout << setprecision(10) << fixed;
    get_input();
    auto sol = solve();
    sol.output();
    cerr << sol.eval() << endl;
}