結果
| 問題 | No.5007 Steiner Space Travel |
| ユーザー |
 wanui
|
| 提出日時 | 2022-07-30 18:03:09 |
| 言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
| 結果 |
AC
|
| 実行時間 | 941 ms / 1,000 ms |
| コード長 | 15,528 bytes |
| コンパイル時間 | 3,594 ms |
| 実行使用メモリ | 5,156 KB |
| スコア | 8,031,162 |
| 最終ジャッジ日時 | 2022-07-30 18:03:43 |
| 合計ジャッジ時間 | 34,085 ms |
|
ジャッジサーバーID （参考情報） |
judge13 / judge14 |
| 純コード判定しない問題か言語 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 941 ms
3,864 KB |
| testcase_01 | AC | 917 ms
3,968 KB |
| testcase_02 | AC | 937 ms
3,908 KB |
| testcase_03 | AC | 918 ms
4,008 KB |
| testcase_04 | AC | 915 ms
3,824 KB |
| testcase_05 | AC | 940 ms
3,920 KB |
| testcase_06 | AC | 919 ms
3,940 KB |
| testcase_07 | AC | 916 ms
3,964 KB |
| testcase_08 | AC | 916 ms
3,864 KB |
| testcase_09 | AC | 916 ms
3,968 KB |
| testcase_10 | AC | 941 ms
4,100 KB |
| testcase_11 | AC | 918 ms
4,056 KB |
| testcase_12 | AC | 915 ms
4,156 KB |
| testcase_13 | AC | 916 ms
3,820 KB |
| testcase_14 | AC | 915 ms
3,828 KB |
| testcase_15 | AC | 940 ms
4,004 KB |
| testcase_16 | AC | 917 ms
3,820 KB |
| testcase_17 | AC | 920 ms
4,108 KB |
| testcase_18 | AC | 917 ms
3,832 KB |
| testcase_19 | AC | 917 ms
4,012 KB |
| testcase_20 | AC | 941 ms
3,972 KB |
| testcase_21 | AC | 918 ms
3,864 KB |
| testcase_22 | AC | 926 ms
4,104 KB |
| testcase_23 | AC | 918 ms
4,900 KB |
| testcase_24 | AC | 916 ms
4,900 KB |
| testcase_25 | AC | 940 ms
4,908 KB |
| testcase_26 | AC | 917 ms
5,156 KB |
| testcase_27 | AC | 929 ms
4,904 KB |
| testcase_28 | AC | 916 ms
4,900 KB |
| testcase_29 | AC | 915 ms
4,900 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};
}
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};
}
void setdists(vector<vector<int>>& dists, vector<point>& stations) {
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;
}
}
}
void update_dists(vector<vector<int>>& dists, vector<point>& stations, int tgt_sid) {
for (int a = 0; a < N; a++) {
dists[a][a] = 0;
for (int b = a + 1; b < N; b++) {
int mindist = dists[a][b];
auto pa = planets[a];
auto pb = planets[b];
for (int sid = 0; sid < M; sid++) {
auto ss = stations[sid];
if (sid == tgt_sid) {
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 経由
mindist = min(mindist, d);
}
} else {
int tid = tgt_sid;
auto st = stations[tid];
int d = ALPHA * distance2(pa, ss) + distance2(ss, st) + ALPHA * distance2(st, pb); // station s-t 経由
mindist = min(mindist, d);
}
}
dists[a][b] = mindist;
dists[b][a] = mindist;
}
}
}
pair<vector<pii>, vector<point>> updates(vector<pii> init_ops, vector<point> init_stations) {
// double start_time = now();
auto stations = init_stations;
double begin_time = now();
double end_time = END_TIME;
double begin_temp = 10.0;
double end_temp = 0.5;
// stationの位置が完全に定まっているとする
// 2つのplanet間の距離を定めることができる (経由する最大station数を2とする)
// これはもうただのTSP
// 初期解はk-meansで求めたops,stations
vector<int> line;
for (auto op : init_ops) {
if (op.first == PLANET) {
line.push_back(op.second);
}
}
vector<vector<int>> dists(N, vector<int>(N));
setdists(dists, stations);
ll old_score = calc_tsp_score(line, dists);
int iter;
for (iter = 1; now() < end_time; iter++) {
if (iter % 200 != 0) {
// 2-opt
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);
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);
}
} else {
// stationの位置を微調整
int m = engine() % M;
int dx = randint(-10, 10);
int dy = randint(-10, 10);
auto newstations = stations;
newstations[m].i += dx;
newstations[m].j += dy;
update_dists(dists, newstations, m); // TODO 差分更新
ll new_score = calc_tsp_score(line, dists);
if (newstations[m].i < 0 || newstations[m].i > 1000 || newstations[m].j < 0 || newstations[m].j > 1000) new_score = -1e9;
if (anneal_accept(new_score, old_score, now(), begin_time, end_time, begin_temp, end_temp)) {
// if (new_score >= old_score) {
old_score = new_score;
swap(stations, newstations);
} else {
update_dists(dists, stations, m);
}
}
}
perr("anneal_iter=", iter);
// 復元
vector<pii> ops;
{
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});
}
// stationを微調整
{
ll old_score = -calc_cost(ops, stations);
for (int iter = 0; iter < 20000; iter++) {
int m = engine() % M;
int dx = randint(-5, 5);
int dy = randint(-5, 5);
stations[m].i += dx;
stations[m].j += dy;
ll new_score = -calc_cost(ops, stations);
if (stations[m].i < 0 || stations[m].i > 1000 || stations[m].j < 0 || stations[m].j > 1000) new_score = -1e9;
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;
}
}
}
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;
}