#include // clang-format off using namespace std; using ll=long long; using ull=unsigned long long; using pll=pair; const ll INF=4e18; void print0(){}; template void print0(H h,T... t){cout<void print(H h,T... t){print0(h);if(sizeof...(T)>0)print0(" ");print(t...);} void perr0(){}; template void perr0(H h,T... t){cerr<void perr(H h,T... t){perr0(h);if(sizeof...(T)>0)perr0(" ");perr(t...);} void ioinit() { cout<; 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 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 ops, vector 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> kmeans() { // k-means // 初期値 vector perm(N); iota(perm.begin(), perm.end(), 0); shuffle(perm.begin(), perm.end(), engine); vector stations; for (int i = 0; i < M; i++) { stations.push_back(planets[perm[i]]); } vector pre_assign(N, -1); for (int iter = 0; iter < 10000; iter++) { vector 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 sums(M, {0, 0}); vector 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 cluster_order(M); // クラスタをたどる順を全探索で探す(bitdpの方が速い) { int sid_init = pre_assign[0]; vector 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 ops; { ops.push_back({PLANET, 0}); for (auto sid : cluster_order) { ops.push_back({STATION, sid}); // 座圧? vector i2p; map 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 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& line, vector>& 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> get_mindist(point pa, point pb, int paid, vector& stations) { // 2つのplanet間の最短距離 (経由する最大station数を2とする) int mindist = ALPHA * ALPHA * distance2(pa, pb); // 直接行く tuple 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>& dists, vector& 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>& dists, vector& 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> updates(vector init_ops, vector 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 line; for (auto op : init_ops) { if (op.first == PLANET) { line.push_back(op.second); } } vector> dists(N, vector(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 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 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 best_km_ops; vector best_km_stations; int best_km_cost = 1e9; while (now() < END_TIME * 0.2) { vector ops; vector 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 best_ops; vector 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; }