結果
| 問題 | No.5007 Steiner Space Travel |
| ユーザー |
 Jiro_tech15
|
| 提出日時 | 2023-04-30 17:53:37 |
| 言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
| 結果 |
TLE
|
| 実行時間 | - |
| コード長 | 29,623 bytes |
| コンパイル時間 | 4,741 ms |
| コンパイル使用メモリ | 210,144 KB |
| 実行使用メモリ | 154,980 KB |
| スコア | 8,037,027 |
| 最終ジャッジ日時 | 2023-04-30 17:54:14 |
| 合計ジャッジ時間 | 36,569 ms |
|
ジャッジサーバーID （参考情報） |
judge13 / judge12 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | TLE | - |
| testcase_01 | AC | 928 ms
151,664 KB |
| testcase_02 | AC | 940 ms
152,020 KB |
| testcase_03 | AC | 972 ms
152,276 KB |
| testcase_04 | AC | 982 ms
152,528 KB |
| testcase_05 | AC | 931 ms
151,800 KB |
| testcase_06 | AC | 987 ms
152,948 KB |
| testcase_07 | AC | 986 ms
152,856 KB |
| testcase_08 | AC | 954 ms
152,012 KB |
| testcase_09 | AC | 992 ms
153,376 KB |
| testcase_10 | AC | 928 ms
151,820 KB |
| testcase_11 | AC | 972 ms
152,596 KB |
| testcase_12 | AC | 982 ms
152,724 KB |
| testcase_13 | AC | 967 ms
152,428 KB |
| testcase_14 | AC | 922 ms
151,668 KB |
| testcase_15 | AC | 921 ms
150,936 KB |
| testcase_16 | TLE | - |
| testcase_17 | AC | 980 ms
152,924 KB |
| testcase_18 | AC | 991 ms
152,720 KB |
| testcase_19 | AC | 960 ms
152,848 KB |
| testcase_20 | AC | 985 ms
152,868 KB |
| testcase_21 | AC | 961 ms
152,312 KB |
| testcase_22 | AC | 983 ms
152,956 KB |
| testcase_23 | TLE | - |
| testcase_24 | AC | 887 ms
150,752 KB |
| testcase_25 | AC | 936 ms
151,920 KB |
| testcase_26 | AC | 938 ms
151,932 KB |
| testcase_27 | AC | 968 ms
152,852 KB |
| testcase_28 | AC | 961 ms
152,796 KB |
| testcase_29 | AC | 970 ms
152,272 KB |
コンパイルメッセージ
main.cpp: 静的メンバ関数 ‘static void Pos::StaticInit(int)’ 内:
main.cpp:426:24: 警告: ‘dir’ may be used uninitialized [-Wmaybe-uninitialized]
426 | move_to[dir][p] = {adj_x, adj_y};
| ^
main.cpp:413:17: 備考: ‘dir’ はここで宣言されています
413 | Dir dir;
| ^~~
ソースコード
namespace atcoder {}
#ifdef LOCAL
#define dbg(x) cerr << __LINE__ << " : " << #x << " = " << (x) << endl;
#else
#define NDEBUG
#define dbg(x) true;
#pragma GCC target("avx2")
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")
#endif
#ifdef GTEST
#include <gtest/gtest.h>
#endif
#include <math.h>
#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <cstdlib>
#include <cstring>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <list>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <vector>
#ifdef PERF
#include <gperftools/profiler.h>
#endif
using namespace std;
using namespace atcoder;
#define fast_io \
ios_base::sync_with_stdio(false); \
cin.tie(0); \
cout.tie(0);
#define ll long long int
#define rep(i, n) for (int i = 0; i < (int)(n); i++)
#define reps(i, n) for (int i = 1; i <= (int)(n); i++)
#define REP(i, n) for (int i = n - 1; i >= 0; i--)
#define REPS(i, n) for (int i = n; i > 0; i--)
#define MOD (long long int)(1e9 + 7)
#define INF (int)(1e9)
#define LINF (long long int)(1e18)
#define chmax(a, b) a = (((a) < (b)) ? (b) : (a))
#define chmin(a, b) a = (((a) > (b)) ? (b) : (a))
#define all(v) v.begin(), v.end()
typedef pair<int, int> Pii;
typedef pair<ll, ll> Pll;
constexpr double PI = acos(-1);
#ifdef NDEBUG
#define CHECK(v1, op, v2)
#else
#define CHECK(v1, op, v2) \
if (!((v1)op(v2))) { \
cerr << "ERROR:" << (v1) << " " << (v2) << endl; \
assert((v1)op(v2)); \
}
#endif
long double nCr(const int n, const int r) {
long double ret = 1;
rep(t, r) {
ret *= (n - t);
ret /= (r - t);
}
return ret;
}
template <typename T>
string to_string(const vector<T>& vec) {
string ret = "";
rep(i, vec.size()) {
ret += vec[i].to_string();
if (i + 1 != vec.size()) {
ret += ",";
}
}
return ret;
}
template <typename T>
ostream& operator<<(ostream& os, const vector<T>& vec) {
os << to_string(vec);
return os;
}
uint32_t xorshift() {
static uint32_t x = 12345789;
static uint32_t y = 362436069;
static uint32_t z = 521288629;
static uint32_t w = 88675123;
uint32_t t;
t = x ^ (x << 11);
x = y;
y = z;
z = w;
w ^= t ^ (t >> 8) ^ (w >> 19);
return w;
}
int rand(const uint32_t l, const uint32_t r) {
return xorshift() % (r - l) + l;
}
uint32_t rand_other_than(const uint32_t l, const uint32_t r,
const uint32_t other) {
const uint32_t num = rand(l, r - 1);
return num + (num >= other);
}
template <typename T>
const T& rand_vec(const vector<T>& vec) {
assert(vec.size() > 0);
return vec[rand(0, vec.size())];
}
template <typename T>
void shuffle(vector<T>& vec) {
rep(l, (int)vec.size() - 1) {
const int idx = rand(l, vec.size());
swap(vec[idx], vec[l]);
}
}
/* start */
vector<double> PARAMS = {};
class Timer {
chrono::system_clock::time_point _start, _end;
ll _sum = 0, _count = 0;
public:
void start() { _start = chrono::system_clock::now(); }
void stop() { _end = chrono::system_clock::now(); }
void add() {
const chrono::system_clock::time_point now = chrono::system_clock::now();
_sum += static_cast<double>(
chrono::duration_cast<chrono::nanoseconds>(now - _start).count());
_count++;
}
ll sum() const { return _sum / 1000; }
int count() const { return _count; }
string average() const {
if (_count == 0) {
return "NaN";
}
return to_string(_sum / 1000 / _count);
}
void reset() {
_start = chrono::system_clock::now();
_sum = 0;
_count = 0;
}
inline int ms() const {
const chrono::system_clock::time_point now = chrono::system_clock::now();
return static_cast<double>(
chrono::duration_cast<chrono::microseconds>(now - _start).count() /
1000);
}
inline int ns() const {
const chrono::system_clock::time_point now = chrono::system_clock::now();
return static_cast<double>(
chrono::duration_cast<chrono::microseconds>(now - _start).count());
}
};
#ifdef LOCAL
struct Timers : unordered_map<string, Timer> {
friend ostream& operator<<(ostream& os, const Timers& timers) {
for (const auto& pa : timers) {
os << pa.first << " time: " << pa.second.sum() / 1000
<< " count: " << pa.second.count() << endl;
}
return os;
}
};
#else
struct Timers {
struct Dummy {
void start() const {}
void add() const {}
};
Dummy dummy;
const Dummy& operator[](const std::string& str) { return dummy; }
friend ostream& operator<<(ostream& os, const Timers& timers) { return os; }
};
#endif
Timers global_timers;
// NBD = neighborhood
template <class Solution, class NBDGenerator>
struct SimulatedAnnealing {
public:
SimulatedAnnealing(double end_time, double start_temp, double end_temp)
: end_time_(end_time * 1000), // ns
inv_time_(1.0 / end_time_),
cur_time_(0.0),
start_temp_(start_temp),
end_temp_(end_temp),
diff_temp_(end_temp - start_temp),
cur_temp_(start_temp) {
assert(start_temp >= end_temp);
timer_.start();
}
Solution run(const Solution& initial_sol, NBDGenerator& nbd_generator) {
Solution sol(initial_sol);
int acc = 0;
int g = 0;
for (;; ++g) {
if ((g & 0xf) == 0) {
#ifdef SA_MAX_G
if (g >= SA_MAX_G) {
break;
}
#endif
UpdateTime();
UpdateTemp();
if (cur_time_ > end_time_) {
break;
}
}
const auto nbd = nbd_generator.Generate(&sol, g);
constexpr int mod = 1'000'000'000;
const int r = rand(0, mod);
if (static_cast<double>(r) / mod < exp(-nbd->Diff() / cur_temp_)) {
acc++;
nbd->Update(&sol);
} else {
nbd->Restore(&sol);
}
}
cerr << "g,acc: " << g << " " << acc << endl;
return sol;
}
private:
double end_time_, inv_time_, cur_time_;
double start_temp_, end_temp_, diff_temp_, cur_temp_;
Timer timer_;
void UpdateTime() { cur_time_ = timer_.ns(); }
void UpdateTemp() {
cur_temp_ = max(0.0, start_temp_ - cur_time_ * diff_temp_ * inv_time_);
}
};
/* start */
template <typename Point, typename Dist_T>
class Graph_T {
public:
struct Edge {
Point to;
Dist_T dist;
Edge(const Point _to, const Dist_T _dist) : to(_to), dist(_dist) {}
};
Graph_T() {}
explicit Graph_T(const int n)
: n_(n), adjs_(n), dists_(n, vector<Dist_T>(n)) {}
int N() const { return n_; }
Dist_T Dist(const Point& from, const Point& to) const {
return dists_[from][to];
}
void AddEdge(const Point a, const Point b, const Dist_T dist) {
adjs_[a].emplace_back(b, dist);
dists_[a][b] = dist;
}
const vector<Edge>& Adj(const Point a) const { return adjs_[a]; }
void RemoveEdge(const Point a, const Point b) {
{
const auto itr = std::find_if(
all(adjs_[a]), [&](const Edge& edge) { return edge.to == b; });
assert(itr != adjs_[a].end());
adjs_[a].erase(itr);
dists_[a][b] = INF;
}
}
private:
int n_;
vector<vector<Edge>> adjs_;
vector<vector<Dist_T>> dists_;
};
enum Dir {
// y-1, y+1, x-1, x+1
kU,
kD,
kL,
kR
};
struct Pos {
int idx_;
Pos() {}
explicit Pos(const int _idx) : idx_(_idx) {}
Pos(int _x, int _y) : idx_(_y * N + _x) { assert(N > 0); }
int X() const { return pos_2_x[*this]; }
int Y() const { return pos_2_y[*this]; }
int Idx() const { return idx_; }
operator int() const { return Idx(); }
operator size_t() const { return Idx(); }
const vector<Pos>& Adj() const { return adj_poses[*this]; }
const vector<Dir>& AdjDirs() const { return adj_dirs[*this]; }
int Manhattan(const Pos& other) const {
return abs(X() - other.X()) + abs(Y() - other.Y());
}
int Euclid2(const Pos& other) const {
const int dx = X() - other.X();
const int dy = Y() - other.Y();
return dx * dx + dy * dy;
}
bool Move(const Dir dir) {
if (move_to[dir][*this].IsDummy()) {
return false;
} else {
*this = move_to[dir][*this];
return true;
}
}
bool operator<(const Pos& other) const { return this->Idx() < other.Idx(); }
bool operator==(const Pos& other) const { return this->Idx() == other.Idx(); }
bool operator!=(const Pos& other) const { return this->Idx() != other.Idx(); }
friend ostream& operator<<(ostream& os, const Pos& pos) {
os << pos.X() << " " << pos.Y();
return os;
}
bool IsDummy() const { return this->idx_ < 0; }
static Pos Dummy() {
Pos p;
p.idx_ = -1;
return p;
}
static Pos TryCreate(const int x, const int y, const int z) {
if (y < 0 || y >= N || x < 0 || x >= N) {
return Pos::Dummy();
} else {
return Pos(x, y);
}
}
static void StaticInit(const int n) {
N = n;
N2 = N * N;
N4 = N2 * N2;
adj_poses = vector<vector<Pos>>(N2);
adj_dirs = vector<vector<Dir>>(N2);
move_to = vector<vector<Pos>>(4, vector<Pos>(N2, Pos::Dummy()));
pos_2_y.resize(N2);
pos_2_x.resize(N2);
rep(y, N) {
rep(x, N) {
const Pos p{x, y};
pos_2_y[p] = y;
pos_2_x[p] = x;
for (int dy = -1; dy <= 1; ++dy) {
for (int dx = -1; dx <= 1; ++dx) {
if (abs(dy) + abs(dx) != 1) continue;
const int adj_y = y + dy;
const int adj_x = x + dx;
if (adj_y < 0 || adj_y >= N) continue;
if (adj_x < 0 || adj_x >= N) continue;
adj_poses[p].emplace_back(adj_x, adj_y);
Dir dir;
if (dy == -1) {
dir = kU;
} else if (dy == 1) {
dir = kD;
} else if (dx == -1) {
dir = kL;
} else if (dx == 1) {
dir = kR;
} else {
assert(false);
}
move_to[dir][p] = {adj_x, adj_y};
adj_dirs[p].emplace_back(dir);
}
}
}
}
}
static int N, N2, N4;
static vector<vector<Pos>> adj_poses;
static vector<vector<Dir>> adj_dirs;
static vector<vector<Pos>> move_to;
static vector<int> pos_2_x, pos_2_y;
};
vector<vector<Pos>> Pos::adj_poses;
vector<vector<Dir>> Pos::adj_dirs;
vector<vector<Pos>> Pos::move_to;
vector<int> Pos::pos_2_x, Pos::pos_2_y;
int Pos::N, Pos::N2, Pos::N4;
/* start */
constexpr int kPlanetType = 1;
constexpr int kStationType = 2;
struct Solution {
using Graph = Graph_T<int, int>;
Solution() {
orders.resize(N);
std::iota(all(orders), 0);
stations.resize(M);
rep(i, M) { stations[i] = Pos{rand(0, Pos::N + 1), rand(0, Pos::N + 1)}; }
}
static void StaticInit(const int _N, const int _M,
const vector<Pos>& _planets, const Graph& _graph) {
N = _N;
M = _M;
planets = _planets;
graph = _graph;
}
// [l, r)のreverse
int DiffSwap(const int l, const int r) const {
// l-1, l, r-1, r
// l-1, r-1, l, r
const int i0 = (l - 1 + N) % N;
const int i1 = l % N;
const int i2 = (r - 1 + N) % N;
const int i3 = r % N;
const int prev_dist =
graph.Dist(orders[i0], orders[i1]) + graph.Dist(orders[i2], orders[i3]);
const int next_dist =
graph.Dist(orders[i0], orders[i2]) + graph.Dist(orders[i1], orders[i3]);
return next_dist - prev_dist;
}
void Swap(int l, int r) {
if (l > r) {
swap(l, r);
}
reverse(orders.begin() + l, orders.begin() + r);
}
vector<pair<int, int>> CalcAns() const {
int start_i = 0;
for (start_i = 0; start_i < N; start_i++) {
if (this->orders[start_i] == 0) {
break;
}
}
vector<pair<int, int>> ans;
ans.emplace_back(kPlanetType, this->orders[start_i]);
rep(add_i, N) {
const int i = (start_i + add_i) % N;
const int next_i = (i + 1) % N;
const Pos& from_planet = planets[this->orders[i]];
const Pos& to_planet = planets[this->orders[next_i]];
// stationを検討
int min_dist =
Solution::graph.Dist(this->orders[i], this->orders[next_i]);
vector<int> best_stations;
// i,jの順に経由
rep(i, M) {
const Pos& station_i = this->stations[i];
rep(j, M) {
if (i == j) continue;
const Pos& station_j = this->stations[j];
const int dist = from_planet.Euclid2(station_i) * alpha +
station_i.Euclid2(station_j) +
station_j.Euclid2(to_planet) * alpha;
if (dist < min_dist) {
min_dist = dist;
best_stations = {i, j};
}
}
}
// iのみ経由
rep(i, M) {
const Pos& station = this->stations[i];
const int dist = from_planet.Euclid2(station) * alpha +
station.Euclid2(to_planet) * alpha;
if (dist < min_dist) {
min_dist = dist;
best_stations = {i};
}
}
if (best_stations.size() == 0) {
// 惑星to惑星
int from = this->orders[i];
int to = this->orders[next_i];
vector<int> rev_planets;
while (from != to) {
rev_planets.emplace_back(to);
rep(i, N) {
if (i == to) continue;
if (graph.Dist(from, i) + graph.Dist(i, to) ==
graph.Dist(from, to)) {
to = i;
break;
}
}
}
REP(i, rev_planets.size()) {
ans.emplace_back(kPlanetType, rev_planets[i]);
}
} else {
for (const auto station : best_stations) {
ans.emplace_back(kStationType, station);
}
ans.emplace_back(kPlanetType, this->orders[next_i]);
}
}
return ans;
}
friend ostream& operator<<(ostream& os, const Solution& sol) {
const auto ans = sol.CalcAns();
rep(i, M) { os << sol.stations[i] << endl; }
os << ans.size() << endl;
rep(i, ans.size()) {
const auto& [t, r] = ans[i];
os << t << " " << r + 1;
if (i + 1 != (int)ans.size()) {
os << endl;
}
}
int dist_sum = 0;
rep(i, (int)ans.size() - 1) {
const auto& [t0, r0] = ans[i];
const auto& [t1, r1] = ans[(i + 1) % ans.size()];
const Pos& from = (t0 == kPlanetType ? planets[r0] : sol.stations[r0]);
const Pos& to = (t1 == kPlanetType ? planets[r1] : sol.stations[r1]);
dist_sum += from.Euclid2(to) * (t0 == kPlanetType ? alpha : 1) *
(t1 == kPlanetType ? alpha : 1);
}
cerr << "score: " << std::round(1e9 / (1000 + sqrt(dist_sum))) << endl;
return os;
}
static vector<Pos> planets;
static int N, M;
static Graph graph;
static constexpr int alpha = 5;
static constexpr int alpha2 = 25;
vector<int> orders;
vector<Pos> stations;
};
int Solution::N;
int Solution::M;
vector<Pos> Solution::planets;
Graph_T<int, int> Solution::graph;
/* start */
struct TSPSolution {
using Graph = Graph_T<int, int>;
TSPSolution() {
orders.resize(N);
std::iota(all(orders), 0);
}
TSPSolution(const vector<int>& _orders) : orders(_orders) {}
static void StaticUpdateGraph(const Graph& _graph) {
graph = _graph;
N = _graph.N();
}
// [l, r)のreverse
int DiffSwap(const int l, const int r) const {
// l-1, l, r-1, r
// l-1, r-1, l, r
const int i0 = (l - 1 + N) % N;
const int i1 = l % N;
const int i2 = (r - 1 + N) % N;
const int i3 = r % N;
const int prev_dist =
graph.Dist(orders[i0], orders[i1]) + graph.Dist(orders[i2], orders[i3]);
const int next_dist =
graph.Dist(orders[i0], orders[i2]) + graph.Dist(orders[i1], orders[i3]);
return next_dist - prev_dist;
}
void Swap(int l, int r) {
if (l > r) {
swap(l, r);
}
reverse(orders.begin() + l, orders.begin() + r);
}
// toの前にtargetを入れる
int DiffInsert(const int to, const int target) const {
// to-1, to, target-1, target, target, target+1
// to-1, target, target, to, target-1, target+1
const int i0 = (to - 1 + N) % N;
const int i1 = to;
const int i2 = (target - 1 + N) % N;
const int i3 = target;
const int i4 = target;
const int i5 = (target + 1) % N;
const int prev_dist = graph.Dist(orders[i0], orders[i1]) +
graph.Dist(orders[i2], orders[i3]) +
graph.Dist(orders[i4], orders[i5]);
const int next_dist = graph.Dist(orders[i0], orders[i3]) +
graph.Dist(orders[i4], orders[i1]) +
graph.Dist(orders[i2], orders[i5]);
return next_dist - prev_dist;
}
void Insert(const int to, const int target) {
orders.insert(orders.begin() + to, orders[target]);
orders.erase(orders.begin() + target + (to <= target ? 1 : 0));
}
static Graph graph;
static int N;
vector<int> orders;
};
int TSPSolution::N;
Graph_T<int, int> TSPSolution::graph;
/* start */
using Graph = Graph_T<int, int>;
struct NBD {
int Diff() const { return diff_; }
virtual void Update(TSPSolution*) const = 0;
void Restore(TSPSolution*) const {}
int diff_;
};
struct NBD_swap : public NBD {
NBD_swap() {}
NBD_swap(TSPSolution* const sol, const int l, const int r) {
l_ = l;
r_ = r;
diff_ = sol->DiffSwap(l, r);
}
void Update(TSPSolution* sol) const override { sol->Swap(l_, r_); }
int l_;
int r_;
};
struct NBD_insert : public NBD {
NBD_insert() {}
NBD_insert(TSPSolution* const sol, const int to, const int target) {
to_ = to;
target_ = target;
diff_ = sol->DiffInsert(to, target);
}
void Update(TSPSolution* sol) const override { sol->Insert(to_, target_); }
int to_;
int target_;
};
struct NBDGenerator {
NBDGenerator() {}
NBD* Generate(TSPSolution* const sol, const int g) {
if ((g & 0x4) == 0) {
// 1.5-opt
const int to = rand(0, TSPSolution::N);
const int target = rand(0, TSPSolution::N);
nbd_insert0 = NBD_insert(sol, to, target);
return &nbd_insert0;
} else {
// 2-opt
const int l = rand(0, TSPSolution::N);
const int r = rand_other_than(0, TSPSolution::N, l);
nbd_swap0 = NBD_swap(sol, l, r);
return &nbd_swap0;
}
}
NBD_swap nbd_swap0; // 2-opt
NBD_insert nbd_insert0; // 1.5-opt
};
struct State {
static State Dummy() {
State state;
state.sum_dist = INF;
return state;
}
static void StaticInit(const int _N, const int _M,
const vector<Pos>& _planets) {
N = _N;
M = _M;
planets = _planets;
}
bool IsDummy() const { return sum_dist == INF; }
void UpdateWarshall() {
if (stations.size() == 0) return;
const int total_N = N + stations.size();
const int k = total_N - 1;
const Pos last_station = stations.back();
// kにまっすぐ向かった時の距離
vector<int> dist_to_k_straight(total_N);
rep(i, N) {
dist_to_k_straight[i] = planets[i].Euclid2(last_station) * alpha;
}
for (int i = N; i < total_N; ++i) {
dist_to_k_straight[i] = stations[i - N].Euclid2(last_station);
}
// iからkへの最短距離をまずは求める
rep(i, total_N) {
int min_dist = INF;
rep(j, total_N - 1) {
const int dist = warshall[i][j] + dist_to_k_straight[j];
chmin(min_dist, dist);
}
assert(min_dist < INF);
warshall[i][k] = warshall[k][i] = min_dist;
}
// warshallを更新
warshall[k][k] = 0;
rep(i, total_N) {
rep(j, i) {
const int next_dist = warshall[i][k] + warshall[k][j];
if (next_dist < warshall[i][j]) {
warshall[i][j] = warshall[j][i] = next_dist;
}
}
}
}
int CalcSumDist() const {
int ret = 0;
rep(i, orders.size()) {
const int next_i = (i + 1) % orders.size();
ret += warshall[orders[i]][orders[next_i]];
}
return ret;
}
void TwoOpt(const double time_limit) {
SimulatedAnnealing<TSPSolution, NBDGenerator> sa(time_limit, 1.0, 0.03);
NBDGenerator nbd_generator;
Graph graph(N);
rep(i, N) {
rep(j, N) { graph.AddEdge(i, j, warshall[i][j]); }
}
TSPSolution::StaticUpdateGraph(graph);
assert(orders.size() == N);
TSPSolution initial_sol(orders);
const auto tsp_sol = sa.run(initial_sol, nbd_generator);
orders = tsp_sol.orders;
sum_dist = CalcSumDist();
}
vector<vector<int>> warshall;
vector<int> orders;
vector<Pos> stations;
int sum_dist;
int prev_j;
int prev_k;
static constexpr int alpha = 5;
static constexpr int alpha2 = 25;
static vector<Pos> planets;
static int N, M;
};
int State::N;
int State::M;
vector<Pos> State::planets;
class Solver {
public:
Solver(istream& is) {
is >> N >> M;
Pos::StaticInit(1001);
planets.reserve(N);
rep(i, N) {
int a, b;
is >> a >> b;
planets.emplace_back(a, b);
}
warshall_without_stations = CalcWarshallWithoutStations();
const auto graph = InitGraph(warshall_without_stations);
TSPSolution::StaticUpdateGraph(graph);
Solution::StaticInit(N, M, planets, graph);
State::StaticInit(N, M, planets);
}
vector<vector<int>> CalcWarshallWithoutStations() const {
// stationの分も確保しておく
vector<vector<int>> dists(N + M, vector<int>(N + M, INF));
rep(i, N + M) { dists[i][i] = 0; }
rep(i, N) {
rep(j, N) { dists[i][j] = alpha2 * planets[i].Euclid2(planets[j]); }
}
rep(k, N) {
rep(i, N) {
rep(j, N) { chmin(dists[i][j], dists[i][k] + dists[k][j]); }
}
}
return dists;
}
Graph InitGraph(const vector<vector<int>>& dists) const {
// dists.size() = N+M
// graph.size() = N
Graph graph(N);
rep(i, N) {
rep(j, N) {
const int dist = dists[i][j];
graph.AddEdge(i, j, dist);
}
}
return graph;
}
using BS = bitset<100>;
vector<pair<Pos, int>> GridSearchNextStation(const State& state,
const int delta) const {
vector<pair<Pos, int>> best_next_stations(N + 1, {Pos::Dummy(), INF});
assert(state.orders.size() == N);
for (int y = 0; y <= 1000; y += delta) {
for (int x = 0; x <= 1000; x += delta) {
const Pos station(x, y);
// このstationの経由を可能としたときの、各辺のdistを求める
int dist_sum = 0;
int used_station_count = 0;
rep(order_i, N) {
const int i = state.orders[order_i];
const int j = state.orders[(order_i + 1) % N];
const auto pos_i = planets[i];
const auto pos_j = planets[j];
const int used_dist =
(pos_i.Euclid2(station) + pos_j.Euclid2(station)) * alpha;
const int dist = min(state.warshall[i][j], used_dist);
dist_sum += dist;
used_station_count += (dist < warshall_without_stations[i][j]);
}
if (best_next_stations[used_station_count].second > dist_sum) {
best_next_stations[used_station_count] = {station, dist_sum};
}
}
}
return best_next_stations;
}
Solution BeamSearch(const vector<int>& initial_orders) const {
assert(initial_orders.size() == N);
assert(warshall_without_stations.size() == N + M);
vector dp(M + 1, vector(N + 1, State::Dummy()));
dp[0][0].warshall = warshall_without_stations;
dp[0][0].orders = initial_orders;
dp[0][0].sum_dist = dp[0][0].CalcSumDist();
rep(j, M + 1) {
rep(k, N + 1) {
auto& state = dp[j][k];
if (state.IsDummy()) continue;
// TSPを追加で1ms行い、既存の駅を有効活用できている状態にする
if (j != 0) {
state.warshall = dp[state.prev_j][state.prev_k].warshall;
}
state.UpdateWarshall();
state.TwoOpt(0.8);
if (j == M) continue;
// 訪問順は固定した状態で、次の駅を置く場所をグリッドサーチして、
//「何らかの駅を経由している辺の本数」ごとに「total_distが最小」となるような
// 駅の置き場所を求める
const auto best_next_stations = GridSearchNextStation(state, 33);
// それらを用いて遷移
rep(next_k, N + 1) {
const auto& [station, dist] = best_next_stations[next_k];
if (dist == INF) continue;
auto& next_state = dp[j + 1][next_k];
if (next_state.sum_dist <= dist) continue;
next_state.orders = state.orders;
next_state.stations = state.stations;
next_state.stations.emplace_back(station);
next_state.sum_dist = dist;
next_state.prev_j = j;
next_state.prev_k = k;
}
}
}
using IDX = tuple<int, int>;
IDX best_idx;
int min_dist = INF;
const int j = M;
rep(k, N + 1) {
auto& state = dp[j][k];
if (state.IsDummy()) continue;
if (state.sum_dist < min_dist) {
min_dist = state.sum_dist;
best_idx = {j, k};
}
}
#if 0
IDX idx = best_idx;
while (true) {
const auto [j, k] = idx;
const auto& state = dp[j][k];
Solution sol;
sol.orders = state.orders;
sol.stations = state.stations;
sol.stations.resize(M, Pos{0, 0});
cout << endl << endl << sol << endl;
if (j == 0) {
break;
}
idx = {state.prev_j, state.prev_k};
}
#endif
cerr << "min_dist: " << min_dist << endl;
const auto& [best_j, best_k] = best_idx;
Solution best_sol;
best_sol.orders = dp[best_j][best_k].orders;
best_sol.stations = dp[best_j][best_k].stations;
return best_sol;
}
Solution OptimizeStations(const Solution& base_sol) const {
auto sol = base_sol;
rep(_t, 3) {
const auto ans = base_sol.CalcAns();
// 各stationから伸びている辺のdistの総和を最小化する位置に各stationを移動する
// 重心に移動すればよい
// station->stationについては、重みを1/5にするのが自然だが、
// もう一つのstationもずれることが予想されるので、重みは1/10にする
vector<float> weight_sums(M);
vector<float> Xs(M), Ys(M);
rep(i, (int)ans.size() - 1) {
auto [t0, r0] = ans[i];
auto [t1, r1] = ans[i + 1];
Pos p0 = (t0 == kPlanetType ? planets[r0] : sol.stations[r0]);
Pos p1 = (t1 == kPlanetType ? planets[r1] : sol.stations[r1]);
rep(_t, 2) {
swap(t0, t1);
swap(r0, r1);
swap(p0, p1);
// 0について考える
if (t0 != kStationType) continue;
// station->planetなら1.0
// station->stationなら0.1
const float weight = (t1 == kPlanetType ? 1.0 : 0.1);
weight_sums[r0] += weight;
Xs[r0] += p1.X() * weight;
Ys[r0] += p1.Y() * weight;
}
}
rep(i, M) {
const auto weight = weight_sums[i];
assert(weight > 0.0);
float float_x = Xs[i] / weight;
float float_y = Ys[i] / weight;
const int x = std::clamp((int)std::round(float_x), 0, 1000);
const int y = std::clamp((int)std::round(float_y), 0, 1000);
sol.stations[i] = Pos{x, y};
}
// warshallを0から更新
vector<vector<int>> dists(N + M, vector<int>(N + M, INF));
vector<Pos> planets_and_stations(N + M);
rep(i, N) { planets_and_stations[i] = planets[i]; }
for (int i = N; i < N + M; i++) {
planets_and_stations[i] = sol.stations[i - N];
}
rep(i, N + M) { dists[i][i] = 0; }
rep(i, N + M) {
rep(j, i) {
dists[i][j] = dists[j][i] =
planets_and_stations[i].Euclid2(planets_and_stations[j]) *
(i < N ? alpha : 1) * (j < N ? alpha : 1);
}
}
rep(k, N + M) {
rep(i, N + M) {
rep(j, N + M) { chmin(dists[i][j], dists[i][k] + dists[k][j]); }
}
}
const auto graph = InitGraph(dists);
TSPSolution::StaticUpdateGraph(graph);
SimulatedAnnealing<TSPSolution, NBDGenerator> sa(10, 1.0, 0.03);
NBDGenerator nbd_generator;
TSPSolution initial_sol(sol.orders);
const auto tsp_sol = sa.run(initial_sol, nbd_generator);
sol.orders = tsp_sol.orders;
}
return sol;
}
Solution Solve(const int time_limit) const {
const TSPSolution initial_sol;
SimulatedAnnealing<TSPSolution, NBDGenerator> sa(100, 1.0, 0.03);
NBDGenerator nbd_generator;
const auto tsp_sol = sa.run(initial_sol, nbd_generator);
const auto beam_sol = BeamSearch(tsp_sol.orders);
const auto optimized_sol = OptimizeStations(beam_sol);
return optimized_sol;
}
private:
int N, M;
vector<vector<int>> warshall_without_stations;
vector<Pos> planets;
static constexpr int alpha = 5;
static constexpr int alpha2 = 25;
};
int main(int argc, char* argv[]) {
fast_io;
if (argc >= 2) {
int idx = 0;
for (int i = 1; i < argc; ++i) {
PARAMS[idx++] = std::stod(argv[i]);
}
}
Timer timer;
timer.start();
Solver solver(cin);
auto sol = solver.Solve(900 - timer.ms());
cerr << "ms: " << timer.ms() << endl;
cout << sol << endl;
return 0;
}