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

• 問題: No.5007 Steiner Space Travel
• ユーザー: highjump
• 提出日時: 2023-04-27 22:23:35
• 言語: C++17(gcc12) (gcc 12.3.0 + boost 1.87.0)
• 結果: AC
• 実行時間: 954 ms / 1,000 ms
• コード長: 11,473 bytes
• コンパイル時間: 3,102 ms
• コンパイル使用メモリ: 154,680 KB
• 実行使用メモリ: 4,376 KB
• スコア: 8,953,008
• 最終ジャッジ日時: 2023-04-27 22:24:12
• 合計ジャッジ時間: 34,149 ms
• ジャッジサーバーID （参考情報）: judge13 / judge14

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

ソースコード

#pragma GCC optimize("Ofast")
#include <iostream>
#include <algorithm>
#include <string>
#include <cstdlib>
#include <cmath>
#include <iomanip>
#include <cctype>
#include <sstream>
#include <vector>
#include <stack>
#include <deque>
#include <queue>
#include <list>
#include <set>
#include <map>
#include <unordered_map>  
#include <chrono>
#include <random>
#include <bitset>

using namespace std;

using ll = long long;
using P = pair<int, int>;
template <class T>using V = vector<T>;
template <class T>using VV = V<V<T>>;
template <class T, class U>bool chmin(T& t, const U& u) { if (t > u) { t = u; return 1; } else return 0; }
template <class T, class U>bool chmax(T& t, const U& u) { if (t < u) { t = u; return 1; } else return 0; }

#define REP(i,n) for(int i=0;i<int(n);i++)
#define FOR(i,a,b) for(int i=int(a);i<=int(b);i++)
#define FORD(i,a,b) for(int i=int(a);i>=int(b);i--)
#define MP make_pair
#define SZ(x) int(x.size())
#define ALL(x) x.begin(),x.end()
#define INF 10001000
#define endl "\n"

chrono::system_clock::time_point t_start;
const double TIME_LIMIT1 = 100, TIME_LIMIT = 950;
double last_update = 0, t_diff;
double start_temp, end_temp;

mt19937 rng;
using uni_int = uniform_int_distribution<>;
using uni_real = uniform_real_distribution<>;
using u64 = uint64_t;
using u32 = uint32_t;

class XorInt {
public:
    uint32_t x = 2463534242;

    XorInt() {};

    XorInt(uint32_t seed) {
        x = 2463534242 + seed;
    };

    uint32_t next() {
        x ^= x << 13; x ^= x >> 17; x ^= x << 5;
        return x;
    }

    uint32_t operator()(uint32_t l, uint32_t r) {
        return (this->next() % (r - l)) + l;
    }

    uint32_t operator()(uint32_t d) {
        return this->next() % d;
    }
};

XorInt xs;

class Xor64 {
public:
    uint64_t x = 88172645463325252ULL;

    Xor64() {};

    Xor64(uint64_t seed) {
        x = 88172645463325252ULL + seed;
    }

    uint64_t next() {
        x ^= x << 13; x ^= x >> 7; x ^= x << 17;
        return x;
    }

    uint64_t operator()(uint64_t l, uint64_t r) {
        return (this->next() % (r - l)) + l;
    }

    uint64_t operator()(uint64_t d) {
        return this->next() % d;
    }
};

Xor64 xx;

void get_time() {
    auto t_now = chrono::system_clock::now();
    t_diff = chrono::duration_cast<chrono::milliseconds>(t_now - t_start).count();
}

int N, M;
const int A = 5;
int score = 0;
int bestScore = 0;
int cnt = 0,update_cnt=0,cnt2=0;
int tempFuel = 0;


VV<int> candidates;
VV<int> distances;
VV<int> destinations;
V<int> route;

class Pos {
public:
    int x = 0;
    int y = 0;
    int kind = 0;
    int order = 0;

    Pos() {};
    Pos(int sx, int sy) :x(sx), y(sy) {};
    Pos(int sx, int sy, int o) :x(sx), y(sy), order(o) {};


};

class Planet :public Pos {
public:

    Planet() { kind = 1; };
    Planet(int sx, int sy, int o) :Pos(sx, sy, o) { kind = 1; };

};

class Station :public Pos {
public:

    Station() { kind = 2; };
    Station(int sx, int sy, int o) :Pos(sx, sy, o) { kind = 2; };

};

int Diff(int lx, int ly, int rx, int ry) {
    return (lx - rx) * (lx - rx) + (ly - ry) * (ly - ry);
}

int getFuel(Pos& l, Pos& r) {
    int dis = (l.x - r.x) * (l.x - r.x) + (l.y - r.y) * (l.y - r.y);
    if (l.kind != r.kind)
        return dis * A;
    if (l.kind == 1)
        return dis * A * A;
    return dis;
}

V<Planet> planets;
V<Station> stations;

Pos& setPos(int k) {
    if (k < N)return planets[k];
    return stations[k % N];
}

int getScore(V<int>& route) {
    int sumFuel = 0;
    FOR(i, 1, route.size() - 1) {
        Pos& current = setPos(route[i - 1]);
        Pos& target = setPos(route[i]);
        sumFuel += getFuel(current, target);
    }

    return int(round(1000000000.0 / (1000.0 + sqrt(sumFuel))));
}

void calcCenter(V<int>& labels, V<P>& centers) {
    V<int> counter(M);
    REP(i, M)centers[i] = MP(0, 0);
    REP(i, N) {
        Planet& planet = planets[i];
        counter[labels[i]]++;
        centers[labels[i]].first += planet.x;
        centers[labels[i]].second += planet.y;
    }
    REP(i, M) {
        if (counter[i] == 0) {
            continue;
        }
        centers[i].first /= counter[i];
        centers[i].second /= counter[i];
    }
}

int reallocate(V<int>& labels, V<P>& centers) {
    int res = 0;
    REP(i, N) {
        Planet& planet = planets[i];
        int newLabel = labels[i];
        int minDiff = INF;
        REP(label, M) {
            if (chmin(minDiff, Diff(planet.x, planet.y, centers[label].first, centers[label].second)))
                newLabel = label;
        }
        if (newLabel != labels[i])res++;
        labels[i] = newLabel;
    }
    return res;
}

void kmeans(V<int>& xpos, V<int>& ypos) {
    V<int> labels(N);
    REP(i, N)labels[i] = i%M;
    V<P> centers(M);
    calcCenter(labels, centers);
    REP(i, 300) {
        int diff=reallocate(labels, centers);
        if (diff == 0)break;
        calcCenter(labels, centers);
    }
    REP(i, M) {
        if (centers[i].first == 0)continue;
        xpos[i] = centers[i].first;
        ypos[i] = centers[i].second;
    }
}

void initwfroid() {
    REP(i, N)REP(j, N) {
        destinations[i][j] = j;
        if (i == j) {
            distances[i][j] = 0;
            continue;
        }
        distances[i][j] = getFuel(setPos(i), setPos(j));
    }

    REP(k, N) {
        REP(i, N) {
            REP(j, N) {
                if (chmin(distances[i][j], distances[i][k] + distances[k][j]))
                    destinations[i][j] = destinations[i][k];
            }
        }
    }
}

int calc() {
    int res = 0;
    REP(i, N) {
        res += distances[route[i]][route[(i + 1) % N]];
    }
    return res;
}

void wfroid(int range) {
    REP(i, range)REP(j, range) {
        destinations[i][j] = j;
        if (i == j) {
            distances[i][j] = 0;
            continue;
        }
        distances[i][j] = getFuel(setPos(i), setPos(j));
    }

    REP(k, range) {
        REP(i, range) {
            REP(j, range) {
                if (chmin(distances[i][j], distances[i][k] + distances[k][j]))
                    destinations[i][j] = destinations[i][k];
            }
        }
    }
}

void setup() {
    candidates = VV<int>(N);
    planets = V<Planet>(N);
    stations = V<Station>(M);
    distances = VV<int>(N + M, V<int>(N + M,INF));
    destinations = VV<int>(N + M, V<int>(N + M));
    route = V<int>(N);

    REP(i, N) {
        int x, y; cin >> x >> y;
        planets[i] = Planet(x, y, i + 1);
    }

    start_temp = 10000;
    end_temp = 10;
    V<int> xpos = { 250,500,750,250,750,250,500,750 };
    V<int> ypos = { 250,250,250,500,500,750,750,750 };

    //kmeans(xpos, ypos);

    REP(i, M) {
        int x = xpos[i];
        int y = ypos[i];
        stations[i] = Station(x, y, i + 1);
    }

    REP(i, N) {
        Planet& base = planets[i];
        FOR(j, 0, N-1) {
            if (i == j)continue;
            Planet& target = planets[j];
            if (abs(base.x - target.x)* abs(base.x - target.x) + abs(base.y - target.y)* abs(base.y - target.y) <=800*800) {
                candidates[i].push_back(j);
            }
        }
    }

    initwfroid();

    REP(i, N) {
        route[i] = i;
    }
    tempFuel = calc();
}

void findPosition(int& bestx, int& besty) {
    int minDiff = INF;
    REP(i, 100) {
        REP(j, 100) {
            int d = 0;
            int x = 10 * i;
            int y = 10 * j;
            Station station(x, y, 0);
            REP(r, N) {
                int temp = route[r];
                int next = route[(r + 1) % N];
                int dis = getFuel(planets[temp], station) + getFuel(station, planets[next]);
                if (distances[temp][next] > dis) {
                    d += dis - distances[temp][next];
                }
            }
            if (chmin(minDiff, d)) {
                bestx = x;
                besty = y;
            }
        }
    }
}

void install() {
    REP(i, M) {
        int bestx = 500; int besty = 500;
        findPosition(bestx, besty);
        stations[i] = Station(bestx, besty, i + 1);
        wfroid(N + i + 1);
    }
    tempFuel = calc();
}

void twoOpt(double start,double end) {
    int a = xs(N - 1);
    int b = a + 1;

    int c = xs(N - 3);
    if (c >= a) {
        c += 2;
    }
    int d = c + 1;

    if (a > c) {
        swap(a, c);
        swap(b, d);
    }

    int diff = distances[route[a]][route[c]] + distances[route[b]][route[d]] 
                - distances[route[a]][route[b]] - distances[route[c]][route[d]];

    ////焼きなまし
    double temp = start_temp + (end_temp - start_temp) * (t_diff-start) / (end-start);
    double prob = exp(-diff / temp);//最大化
    if (prob > (xs(INF)) / (double)INF) {
    //if (diff<0) {//山登り
        if (c - b < N - d + a) {
            int l = b, r = c;
            REP(i,(c-b+1)/2) {
                swap(route[l], route[r]);
                l++; r--;
            }
        }
        else {
            int l = d, r = a;
            REP(i,(N-d+a+1)/2){
                swap(route[l], route[r]);
                l=(l+1)%N;
                r=(r+N-1)%N;
            }
        }
        tempFuel += diff;
        update_cnt++;
    }
}

void shift() {
    int target = xs(M);
    Station station = stations[target];
    int dx = xs(60) - 30;
    int dy = xs(60) - 30;
    stations[target] = Station(min(max(station.x + dx, 0), 1000), min(max(station.y + dy, 0), 1000), target + 1);
    wfroid(N + M);
    int newFuel = calc();
    if (tempFuel > newFuel) {
        tempFuel = newFuel;
        REP(i, 10000) {
            twoOpt(0, TIME_LIMIT);
        }
    }
    else {
        stations[target] = station;
    }
}

void justification(double limit) {
    while (t_diff < limit) {
        shift();
        cnt2++;
        get_time();
    }
}

void makeAnswer(V<int>& new_route) {
    int start = 0;
    REP(i, N) {
        if (route[i] == 0) {
            start = i;
            break;
        }
    }
    new_route.push_back(0);
    FOR(i,1,N){
        int temp = route[(start+i - 1)%N];
        int next = route[(start+i) % N];
        int mid = destinations[temp][next];
        while (true) {
            new_route.push_back(mid);
            if (mid == next)break;
            mid = destinations[mid][next];
        }
    }
}

void out(V<int>& route) {
    REP(i, M) {
        cout << stations[i].x << " " << stations[i].y << endl;
    }
    cout << route.size() << endl;
    for (auto p : route) {
        Pos& pos = setPos(p);
        cout << pos.kind << " " << pos.order << endl;
    }
}

signed main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    cout << fixed << setprecision(15);
    cerr << fixed << setprecision(15);
    t_start = chrono::system_clock::now();

    get_time();
    cin >> N >> M;

    setup();
    //route = chokudaiSearch(1);
    while (t_diff < 100) {
        twoOpt(0,100);
        cnt++;
        if (cnt % 1000 == 0)get_time();
    }

    install();

    while (t_diff < 300) {
        twoOpt(100, 300);
        cnt++;
        if (cnt % 1000 == 0)get_time();
    }

    justification(TIME_LIMIT);



    V<int>new_route;
    makeAnswer(new_route);

    score = getScore(new_route);
    out(new_route);
    get_time();

    cerr << "time=" << int(t_diff) << " ms" << endl;
    cerr << "score=" << score << endl;
    //cerr << "last=" << last_update << endl;
    cerr << "cnt=" << cnt << endl;
    cerr << "update=" << update_cnt << endl;
    cerr << "cnt2=" << cnt2 << endl;

    return 0;
}