結果

問題 No.5007 Steiner Space Travel
ユーザー Jiro_tech15Jiro_tech15
提出日時 2023-04-30 12:18:19
言語 C++17
(gcc 12.3.0 + boost 1.83.0)
結果
TLE  
実行時間 -
コード長 21,689 bytes
コンパイル時間 4,706 ms
コンパイル使用メモリ 215,116 KB
実行使用メモリ 215,276 KB
スコア 411,212
最終ジャッジ日時 2023-04-30 12:19:04
合計ジャッジ時間 39,370 ms
ジャッジサーバーID
(参考情報)
judge14 / judge11
このコードへのチャレンジ
(要ログイン)

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 TLE -
testcase_01 TLE -
testcase_02 TLE -
testcase_03 TLE -
testcase_04 TLE -
testcase_05 TLE -
testcase_06 TLE -
testcase_07 TLE -
testcase_08 AC 987 ms
198,216 KB
testcase_09 TLE -
testcase_10 TLE -
testcase_11 TLE -
testcase_12 TLE -
testcase_13 AC 983 ms
205,964 KB
testcase_14 TLE -
testcase_15 TLE -
testcase_16 TLE -
testcase_17 TLE -
testcase_18 AC 994 ms
197,212 KB
testcase_19 TLE -
testcase_20 TLE -
testcase_21 TLE -
testcase_22 TLE -
testcase_23 TLE -
testcase_24 TLE -
testcase_25 TLE -
testcase_26 TLE -
testcase_27 TLE -
testcase_28 TLE -
testcase_29 TLE -
権限があれば一括ダウンロードができます
コンパイルメッセージ
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;
      |                 ^~~

ソースコード

diff #

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),
        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_.ms(); }
  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);
  }

  friend ostream& operator<<(ostream& os, const Solution& sol) {
    int start_i = 0;
    for (start_i = 0; start_i < N; start_i++) {
      if (sol.orders[start_i] == 0) {
        break;
      }
    }

#if 0
    // その順序固定だとして、10^6通りの駅について、「その駅を通ったほうが嬉しい」辺のbitsetは何通りあるか?
    using BS = bitset<100>;
    unordered_set<BS> S;
    rep(y, 1001) {
      rep(x, 1001) {
        const Pos station(x, y);
        BS bs = 0;
        rep(i, N) {
          const int next_i = (i + 1) % N;
          const auto p0 = planets[sol.orders[i]];
          const auto p1 = planets[sol.orders[next_i]];

          const int default_dist =
              graph.Dist(sol.orders[i], sol.orders[next_i]);
          const int station_dist =
              p0.Euclid2(station) * alpha + station.Euclid2(p1) * alpha;
          bs[i] = station_dist <= default_dist;
        }
        S.insert(bs);
      }
    }
    cerr << endl;
    cerr << "size: " << S.size() << endl;
#endif

    vector<pair<int, int>> ans;
    ans.emplace_back(kPlanetType, sol.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[sol.orders[i]];
      const Pos& to_planet = planets[sol.orders[next_i]];

      // stationを検討
      int min_dist = Solution::graph.Dist(sol.orders[i], sol.orders[next_i]);
      vector<int> best_stations;

      // i,jの順に経由
      rep(i, M) {
        const Pos& station_i = sol.stations[i];
        rep(j, M) {
          if (i == j) continue;
          const Pos& station_j = sol.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 = sol.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 = sol.orders[i];
        int to = sol.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, sol.orders[next_i]);
      }
    }

    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 NBD {
  int Diff() const { return diff_; }
  virtual void Update(Solution*) const = 0;
  void Restore(Solution*) const {}

  int diff_;
};

struct NBD_swap : public NBD {
  NBD_swap() {}
  NBD_swap(Solution* const sol, const int l, const int r) {
    l_ = l;
    r_ = r;
    diff_ = sol->DiffSwap(l, r);
  }
  void Update(Solution* sol) const override { sol->Swap(l_, r_); }

  int l_;
  int r_;
};

struct NBDGenerator {
  NBDGenerator() {}
  NBD* Generate(Solution* const sol, const int g) {
    // swap
    const int l = rand(0, Solution::N);
    const int r = rand_other_than(0, Solution::N, l);
    nbd_swap0 = NBD_swap(sol, l, r);
    return &nbd_swap0;
  }

  NBD_swap nbd_swap0;
};

class Solver {
 public:
  using Graph = Graph_T<int, int>;
  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);
    }

    graph = InitGraph();
    Solution::StaticInit(N, M, planets, graph);
  }

  Graph InitGraph() const {
    vector<vector<int>> dists(N, vector<int>(N, INF));
    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]); }
      }
    }
    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>;

  unordered_map<BS, pair<Pos, ll>> CalcBestBSs(
      const vector<int>& orders) const {
    // その順序固定だとして、10^6通りの駅について、「その駅を通ったほうが嬉しい」辺のbitsetを求めて、
    // 各bitsetごとに、距離の改善量が最も大きいものを保持する
    // TODO: 駅to駅を考慮
    // TODO: intにするかも
    unordered_map<BS, pair<Pos, ll>> mp;
    rep(y, 1001) {
      rep(x, 1001) {
        const Pos station(x, y);
        BS bs = 0;
        ll diff_sum = 0;

        rep(i, N) {
          const int next_i = (i + 1) % N;
          const auto p0 = planets[orders[i]];
          const auto p1 = planets[orders[next_i]];

          const int default_dist = graph.Dist(orders[i], orders[next_i]);
          const int station_dist =
              p0.Euclid2(station) * alpha + station.Euclid2(p1) * alpha;
          const ll diff = default_dist - station_dist;
          if (diff > 0) {
            bs[i] = true;
            diff_sum += diff;
          }
        }
        if (diff_sum == 0) continue;
        if (mp[bs].second < diff_sum) {
          mp[bs] = {Pos{x, y}, diff_sum};
        }
      }
    }
    cerr << "size: " << mp.size() << endl;
    return mp;
  }

  Solution DPStations(const Solution& order_sol) const {
    const auto orders = order_sol.orders;
    const auto mp = CalcBestBSs(orders);

    // dp[j][k] := i個目まで見て、j個駅を置いて、bit_countがk個
    using P = pair<BS, ll>;
    using IDX = tuple<int, int, int>;
    vector dp(M + 1, vector(N + 1, P{BS{0}, (int)-INF}));
    vector prev(mp.size() + 1, vector(M + 1, vector(101, IDX{0, 0, 0})));
    dp[0][0].second = 0;

    vector<Pos> pos_candidates;
    pos_candidates.reserve(mp.size());
    int i = -1;
    for (const auto& [bs, pos_diff] : mp) {
      ++i;
      const int bs_count = bs.count();
      const auto& [pos, diff] = pos_diff;
      pos_candidates.emplace_back(pos);

      prev[i + 1] = prev[i];

      REP(j, M) {
        const auto next_j = j + 1;
        REP(k, N + 1 - bs_count) {
          if (dp[j][k].second < 0) continue;
          if ((dp[j][k].first & bs) != 0) continue;
          const auto next_bs = (dp[j][k].first | bs);
          const auto next_k = next_bs.count();
          const auto next_val = dp[j][k].second + diff;
          if (dp[next_j][next_k].second < next_val) {
            dp[next_j][next_k].first = next_bs;
            dp[next_j][next_k].second = next_val;
            prev[i + 1][next_j][next_k] = make_tuple(i, j, k);
          }
        }
      }
    }

    IDX best_idx{0, 0, 0};
    ll best_diff = 0;
    rep(j, M + 1) {
      rep(k, N + 1) {
        const auto diff = dp[j][k].second;
        if (best_diff < diff) {
          best_diff = diff;
          best_idx = make_tuple((int)mp.size(), j, k);
        }
      }
    }

    // 復元
    vector<Pos> stations;
    IDX idx = best_idx;
    while (idx != make_tuple(0, 0, 0)) {
      const auto& [i, j, k] = idx;
      const auto& prev_idx = prev[i][j][k];
      const auto& [prev_i, prev_j, prev_k] = prev_idx;
      if (prev_j != j) {
        stations.emplace_back(pos_candidates[prev_i]);
      }
      idx = prev_idx;
    }

    cerr << "station size: " << stations.size() << endl;
    while (stations.size() < M) {
      stations.emplace_back(Pos{500, 500});
    }

    auto sol = order_sol;
    sol.stations = stations;

#ifndef LOCAL
    cout << sol << endl;
    std::quick_exit(0);
#endif
    return sol;
  }

  Solution Solve(const int time_limit) const {
    const Solution initial_sol;
    SimulatedAnnealing<Solution, NBDGenerator> sa(time_limit - 700, 1, 0.01);
    NBDGenerator nbd_generator;
    const auto order_sol = sa.run(initial_sol, nbd_generator);
    const auto final_sol = DPStations(order_sol);

    return final_sol;
  }

 private:
  int N, M;
  vector<Pos> planets;
  Graph graph;
  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());
  cout << sol << endl;
  return 0;
}
0