#pragma GCC optimize("O3") #pragma GCC optimize("unroll-loops") #pragma GCC optimize("Ofast") #include using namespace std; constexpr bool DEBUG = false; // ========================================== // SA Parameters (調整箇所) // ========================================== using score_t = int64_t; // スコアの型 constexpr bool MAXIMIZE = true; // true: 最大化, false: 最小化 constexpr double TIME_LIMIT = 0.98; constexpr double START_TEMP = 1000.0; constexpr double END_TEMP = 10.0; constexpr int TIME_CHECK_INTERVAL = 1; // 時間チェック頻度 (ビットマスク) constexpr int STATS_INTERVAL = 10 * (TIME_CHECK_INTERVAL + 1); // 統計出力頻度 constexpr bool USE_EXPONENTIAL_DECAY = false; // true: 指数減衰, false: べき乗減衰 constexpr double POWER_DECAY_EXP = 1.0; // べき乗減衰の指数 (1.0で線形, >1で序盤高温維持, <1で序盤急冷) constexpr bool ALLOW_WORSE_MOVES = false; // true: SA, false: 山登り // 近傍操作の定義 enum MoveType { DESTROY_RECONNECT, NUM_MOVES }; // 近傍の名前(enum順に対応) vector MOVE_NAMES = {"DESTROY_RECONNECT"}; // 近傍の重み(整数、合計値は任意) vector MOVE_WEIGHTS = {1.0}; // 0: そのまま, 1: sqrt, 2: 二乗 #ifndef WEIGHT_STYLE #define WEIGHT_STYLE 0 #endif // Macros #define el '\n' #define all(v) (v).begin(), (v).end() using i8 = int8_t; using u8 = uint8_t; using i16 = int16_t; using u16 = uint16_t; using i32 = int32_t; using u32 = uint32_t; using i64 = int64_t; using u64 = uint64_t; using f32 = float; using f64 = double; #define rep(i, n) for (i64 i = 0; i < (i64)(n); i++) template using min_queue = priority_queue, greater>; template using max_queue = priority_queue; // Constant const double pi = 3.141592653589793238; const i32 inf32 = 1073741823; const i64 inf64 = 1LL << 60; const string ABC = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"; const string abc = "abcdefghijklmnopqrstuvwxyz"; const int MOD = 998244353; const array dx = {0, 0, -1, 1, -1, -1, 1, 1}; const array dy = {-1, 1, 0, 0, -1, 1, -1, 1}; // Printing template void print_collection(ostream &out, T const &x); template void print_tuple(ostream &out, T const &a, index_sequence); namespace std { template ostream &operator<<(ostream &out, tuple const &x) { print_tuple(out, x, index_sequence_for{}); return out; } template ostream &operator<<(ostream &out, pair const &x) { print_tuple(out, x, index_sequence_for{}); return out; } template ostream &operator<<(ostream &out, array const &x) { print_collection(out, x); return out; } template ostream &operator<<(ostream &out, vector const &x) { print_collection(out, x); return out; } template ostream &operator<<(ostream &out, deque const &x) { print_collection(out, x); return out; } template ostream &operator<<(ostream &out, multiset const &x) { print_collection(out, x); return out; } template ostream &operator<<(ostream &out, multimap const &x) { print_collection(out, x); return out; } template ostream &operator<<(ostream &out, set const &x) { print_collection(out, x); return out; } template ostream &operator<<(ostream &out, map const &x) { print_collection(out, x); return out; } template ostream &operator<<(ostream &out, unordered_set const &x) { print_collection(out, x); return out; } } // namespace std template void print_tuple(ostream &out, T const &a, index_sequence) { using swallow = int[]; out << '('; (void)swallow{0, (void(out << (I == 0 ? "" : ", ") << get(a)), 0)...}; out << ')'; } template void print_collection(ostream &out, T const &x) { int f = 0; out << '['; for (auto const &i : x) { out << (f++ ? "," : ""); out << i; } out << "]"; } // Random struct RNG { uint64_t s[2]; RNG(u64 seed) { reset(seed); } RNG() { reset(time(0)); } using result_type = u32; static constexpr u32 min() { return numeric_limits::min(); } static constexpr u32 max() { return numeric_limits::max(); } u32 operator()() { return randomInt32(); } static __attribute__((always_inline)) inline uint64_t rotl(const uint64_t x, int k) { return (x << k) | (x >> (64 - k)); } inline void reset(u64 seed) { struct splitmix64_state { u64 s; u64 splitmix64() { u64 result = (s += 0x9E3779B97f4A7C15); result = (result ^ (result >> 30)) * 0xBF58476D1CE4E5B9; result = (result ^ (result >> 27)) * 0x94D049BB133111EB; return result ^ (result >> 31); } }; splitmix64_state sm{seed}; s[0] = sm.splitmix64(); s[1] = sm.splitmix64(); } uint64_t next() { const uint64_t s0 = s[0]; uint64_t s1 = s[1]; const uint64_t result = rotl(s0 * 5, 7) * 9; s1 ^= s0; s[0] = rotl(s0, 24) ^ s1 ^ (s1 << 16); // a, b s[1] = rotl(s1, 37); // c return result; } inline u32 randomInt32() { return next(); } inline u64 randomInt64() { return next(); } inline u32 random32(u32 r) { return (((u64)randomInt32()) * r) >> 32; } inline u64 random64(u64 r) { return randomInt64() % r; } inline u32 randomRange32(u32 l, u32 r) { return l + random32(r - l + 1); } inline u64 randomRange64(u64 l, u64 r) { return l + random64(r - l + 1); } inline double randomDouble() { return (double)randomInt32() / 4294967296.0; } inline double randomDoubleOpen01() { return (randomInt32() + 1.0) / 4294967297.0; } inline float randomFloat() { return (float)randomInt32() / 4294967296.0; } inline double randomRangeDouble(double l, double r) { return l + randomDouble() * (r - l); } inline double randomGaussian(double mean = 0.0, double stddev = 1.0) { double u1 = randomDoubleOpen01(); double u2 = randomDouble(); double z = sqrt(-2.0 * log(u1)) * cos(2.0 * 3.141592653589793238 * u2); return mean + stddev * z; } template void shuffle(vector &v) { i32 sz = v.size(); for (i32 i = sz; i > 1; i--) { i32 p = random32(i); swap(v[i - 1], v[p]); } } template void shuffle(T *fr, T *to) { i32 sz = distance(fr, to); for (int i = sz; i > 1; i--) { int p = random32(i); swap(fr[i - 1], fr[p]); } } template inline int sample_index(vector const &v) { return random32(v.size()); } template inline T sample(vector const &v) { return v[sample_index(v)]; } } rng; class FenwickWeightedSampler { private: int n_ = 0; vector bit_; vector weight_; double total_ = 0.0; void add_internal(int idx0, double delta) { int i = idx0 + 1; while (i <= n_) { bit_[i] += delta; i += i & -i; } } public: // Complexity: O(n) void init(int n, double initial_weight = 0.0) { n_ = n; bit_.assign(n_ + 1, 0.0); weight_.assign(n_, initial_weight); total_ = 0.0; if (initial_weight != 0.0) { for (int i = 0; i < n_; i++) { add_internal(i, initial_weight); } total_ = initial_weight * n_; } } // Complexity: O(n log n) void build(const vector &weights) { init((int)weights.size(), 0.0); for (int i = 0; i < n_; i++) { weight_[i] = weights[i]; add_internal(i, weights[i]); total_ += weights[i]; } } // Complexity: O(log n) void set_weight(int idx, double new_weight) { assert(0 <= idx && idx < n_); const double delta = new_weight - weight_[idx]; weight_[idx] = new_weight; add_internal(idx, delta); total_ += delta; } // Complexity: O(log n) void add_weight(int idx, double delta) { assert(0 <= idx && idx < n_); weight_[idx] += delta; add_internal(idx, delta); total_ += delta; } // Complexity: O(log n) int sample() const { assert(n_ > 0); assert(total_ > 0.0); const double r = rng.randomDouble() * total_; int idx = 0; double prefix = 0.0; int step = 1; while ((step << 1) <= n_) step <<= 1; for (int k = step; k > 0; k >>= 1) { const int nxt = idx + k; if (nxt <= n_ && prefix + bit_[nxt] <= r) { idx = nxt; prefix += bit_[nxt]; } } if (idx >= n_) idx = n_ - 1; return idx; } int size() const { return n_; } double total_weight() const { return total_; } double weight(int idx) const { assert(0 <= idx && idx < n_); return weight_[idx]; } }; class AliasWeightedSampler { private: int n_ = 0; vector weight_; double total_ = 0.0; vector prob_; vector alias_; bool dirty_ = true; // Complexity: O(n) void rebuild_alias_table() { prob_.assign(n_, 0.0); alias_.assign(n_, 0); total_ = 0.0; for (double w : weight_) total_ += w; assert(total_ > 0.0); vector scaled(n_); vector small; vector large; small.reserve(n_); large.reserve(n_); for (int i = 0; i < n_; i++) { scaled[i] = weight_[i] * n_ / total_; if (scaled[i] < 1.0) { small.push_back(i); } else { large.push_back(i); } } while (!small.empty() && !large.empty()) { const int s = small.back(); small.pop_back(); const int l = large.back(); large.pop_back(); prob_[s] = scaled[s]; alias_[s] = l; scaled[l] -= (1.0 - scaled[s]); if (scaled[l] < 1.0) { small.push_back(l); } else { large.push_back(l); } } while (!large.empty()) { const int i = large.back(); large.pop_back(); prob_[i] = 1.0; alias_[i] = i; } while (!small.empty()) { const int i = small.back(); small.pop_back(); prob_[i] = 1.0; alias_[i] = i; } dirty_ = false; } public: // Complexity: O(n) void init(int n, double initial_weight = 0.0) { n_ = n; weight_.assign(n_, initial_weight); total_ = initial_weight * n_; prob_.clear(); alias_.clear(); dirty_ = true; } // Complexity: O(n) void build(const vector &weights) { n_ = (int)weights.size(); weight_ = weights; dirty_ = true; // Build lazily on sample() to allow batch updates before first use. } // Complexity: O(1), table becomes dirty void set_weight(int idx, double new_weight) { assert(0 <= idx && idx < n_); weight_[idx] = new_weight; dirty_ = true; } // Complexity: O(1), table becomes dirty void add_weight(int idx, double delta) { assert(0 <= idx && idx < n_); weight_[idx] += delta; dirty_ = true; } // Complexity: O(1) after table is ready, O(n) if rebuild is needed int sample() { assert(n_ > 0); if (dirty_) rebuild_alias_table(); const int i = (int)rng.random32((uint32_t)n_); return (rng.randomDouble() < prob_[i]) ? i : alias_[i]; } int size() const { return n_; } double total_weight() { if (dirty_) rebuild_alias_table(); return total_; } double weight(int idx) const { assert(0 <= idx && idx < n_); return weight_[idx]; } }; // Timer struct Timer { chrono::steady_clock::time_point t_begin; Timer() { t_begin = chrono::steady_clock::now(); } void reset() { t_begin = chrono::steady_clock::now(); } float elapsed() const { return chrono::duration(chrono::steady_clock::now() - t_begin) .count(); } } timer; // Util template T &smin(T &x, T const &y) { x = min(x, y); return x; } template T &smax(T &x, T const &y) { x = max(x, y); return x; } template int sgn(T val) { if (val < 0) return -1; if (val > 0) return 1; return 0; } static inline string int_to_string(int val, int digits = 0) { string s = to_string(val); reverse(begin(s), end(s)); while ((int)s.size() < digits) s.push_back('0'); reverse(begin(s), end(s)); return s; } // Debug static inline void debug() { cerr << "\n"; } template void debug(T const &t, V const &...v) { if constexpr (DEBUG) { cerr << t; if (sizeof...(v)) { cerr << ", "; } debug(v...); } } // Bits __attribute__((always_inline)) inline u64 bit(u64 x) { return 1ull << x; } __attribute__((always_inline)) inline void setbit(u64 &a, u32 b, u64 value = 1) { a = (a & ~bit(b)) | (value << b); } __attribute__((always_inline)) inline u64 getbit(u64 a, u32 b) { return (a >> b) & 1; } __attribute__((always_inline)) inline u64 lsb(u64 a) { return __builtin_ctzll(a); } __attribute__((always_inline)) inline int msb(uint64_t bb) { return __builtin_clzll(bb) ^ 63; } // 入出力高速化 struct Init { Init() { ios::sync_with_stdio(0); cin.tie(0); } } init; // テーブルサイズ(温度・logで共通) constexpr int TABLE_SIZE = 1024; // 温度減衰テーブル double temp_table[TABLE_SIZE + 1]; inline void init_temp_table() { for (int i = 0; i <= TABLE_SIZE; i++) { double progress = (double)i / TABLE_SIZE; if constexpr (USE_EXPONENTIAL_DECAY) { // 指数減衰 double ratio = END_TEMP / START_TEMP; temp_table[i] = START_TEMP * pow(ratio, progress); } else { // べき乗減衰: start - (start - end) * progress^x (x=1で線形) double p = pow(progress, POWER_DECAY_EXP); temp_table[i] = START_TEMP - (START_TEMP - END_TEMP) * p; } } } inline double temp_decay(double progress) { // テーブルルックアップ + 線形補間 double index = progress * TABLE_SIZE; int idx = (int)index; if (idx >= TABLE_SIZE) return temp_table[TABLE_SIZE]; double frac = index - idx; return temp_table[idx] * (1.0 - frac) + temp_table[idx + 1] * frac; } // logテーブル double log_table[TABLE_SIZE + 1]; inline void init_log_table() { for (int i = 0; i <= TABLE_SIZE; i++) { double x = (double)i / TABLE_SIZE; if (x == 0.0) { log_table[i] = -10.0; // log(0) = -inf を有限値で近似 } else { log_table[i] = log(x); } } } inline double log_fast(double x) { if constexpr (DEBUG) assert(0.0 < x && x < 1.0); double index = x * TABLE_SIZE; int idx = (int)index; if (idx >= TABLE_SIZE) return log_table[TABLE_SIZE]; double frac = index - idx; return log_table[idx] * (1.0 - frac) + log_table[idx + 1] * frac; } // ガウシアンテーブル(Box-Mullerの事前計算) constexpr int GAUSS_TABLE_SIZE = 131072; double gauss_table[GAUSS_TABLE_SIZE]; inline void init_gauss_table() { for (int i = 0; i < GAUSS_TABLE_SIZE; i++) { double u1 = rng.randomDoubleOpen01(); double u2 = rng.randomDouble(); gauss_table[i] = sqrt(-2.0 * log(u1)) * cos(2.0 * 3.141592653589793238 * u2); } } inline double fast_gaussian(double mean, double stddev) { return mean + stddev * gauss_table[rng.random32(GAUSS_TABLE_SIZE)]; } // ========================================== // Input & State // ========================================== struct Flight { i16 from = -1; i16 to = -1; i16 dep = -1; i16 arr = -1; }; struct PairInfo { i16 src = -1; i16 dst = -1; i64 gain = 0; }; double g_sa_progress = 0.0; constexpr i32 START_MIN = 6 * 60; constexpr i32 END_MIN = 21 * 60; struct Input { i32 N = 0, R = 0, M = 0, K = 0; vector x, y; vector w; vector sq_flights; static i32 parse_time(string const &s) { i32 hh = stoi(s.substr(0, 2)); i32 mm = stoi(s.substr(3, 2)); return hh * 60 + mm; } void input() { string head; cin >> head; if (!head.empty() && (unsigned char)head[0] == 0xEF) { if (head.size() >= 3 && (unsigned char)head[1] == 0xBB && (unsigned char)head[2] == 0xBF) { head = head.substr(3); } } N = stoi(head); cin >> R; x.resize(N); y.resize(N); w.resize(N); rep(i, N) { cin >> x[i] >> y[i] >> w[i]; } cin >> M; sq_flights.resize(M); rep(i, M) { i32 a, b; string s, t; cin >> a >> s >> b >> t; --a; --b; sq_flights[i] = Flight{(i16)a, (i16)b, (i16)parse_time(s), (i16)parse_time(t)}; } cin >> K; } } in; namespace GlobalData { struct EvalTerm { i32 idx = -1; i64 gain = 0; }; vector targets; vector> duration; vector pairs; vector dp_sq_flat; vector dp_ci_flat_buf; vector ci_flights_buf; vector eval_terms; __int128 total_gain_all = 0; } // namespace GlobalData struct State { score_t score = -1; vector> planes; struct RollbackInfo { i32 plane_id = -1; i32 cut_pos = 0; vector tail; }; static inline i32 dp_index(i32 dest, i32 tk, i32 src) { return (dest * (i32)GlobalData::targets.size() + tk) * in.N + src; } static i16 flight_duration(double d) { double raw = 60.0 * d / 800.0 + 40.0; return (i16)(ceil(raw / 5.0) * 5); } static void sort_desc(vector &flights) { const int SLOTS = (END_MIN - START_MIN) / 5 + 1; const int n = (int)flights.size(); if (n <= 1) return; static vector cnt, start, cur; static vector tmp; if ((int)cnt.size() != SLOTS) { cnt.resize(SLOTS); start.resize(SLOTS + 1); cur.resize(SLOTS); } memset(cnt.data(), 0, SLOTS * sizeof(int)); start[0] = 0; // start[s+1] はループ内で cnt から計算、cur もループ内で設定 tmp.resize(n); for (auto const &f : flights) { int s = (f.dep - START_MIN) / 5; cnt[s]++; } for (int s = 0; s < SLOTS; s++) { start[s + 1] = start[s] + cnt[s]; cur[s] = start[s]; } for (auto const &f : flights) { int s = (f.dep - START_MIN) / 5; tmp[cur[s]++] = f; } int out = 0; for (int s = SLOTS - 1; s >= 0; s--) { int L = start[s], R = start[s + 1]; if (R - L > 1) { sort(tmp.begin() + L, tmp.begin() + R, [](Flight const &a, Flight const &b) { return a.arr > b.arr; }); } for (int i = L; i < R; i++) flights[out++] = tmp[i]; } } static void compute_all_dp_flat_sorted(vector const &flights_sorted, vector &out_dp) { const i32 N = in.N; const i32 TK = (i32)GlobalData::targets.size(); const i32 need = N * TK * N; if ((i32)out_dp.size() != need) out_dp.resize(need); fill(out_dp.begin(), out_dp.end(), (i16)-1); rep(dest, N) { rep(k, TK) { i16 *cur = out_dp.data() + dp_index((i32)dest, (i32)k, 0); cur[dest] = GlobalData::targets[k]; for (auto const &f : flights_sorted) { if (cur[f.to] >= f.arr && f.dep > cur[f.from]) { cur[f.from] = f.dep; } } } } } static vector> build_random_plan() { vector> plan(in.K); // 人口重みの累積和 vector cumw(in.N + 1, 0); for (int i = 0; i < in.N; i++) cumw[i + 1] = cumw[i] + city_weight(i); u64 total_w = cumw[in.N]; // 全都市から重みづけランダム選択 auto pick_city = [&]() -> i32 { u64 r = rng.random64(total_w); return (i32)(upper_bound(cumw.begin(), cumw.end(), r) - cumw.begin()) - 1; }; // except を除いた都市から重みづけランダム選択 auto pick_city_except = [&](i32 except) -> i32 { u64 ew = city_weight(except); u64 total_ex = total_w - ew; if (total_ex == 0) return (except + 1) % in.N; u64 r = rng.random64(total_ex); if (r < cumw[except]) { return (i32)(upper_bound(cumw.begin(), cumw.begin() + except + 1, r) - cumw.begin()) - 1; } else { r += ew; return (i32)(upper_bound(cumw.begin() + except + 1, cumw.end(), r) - cumw.begin()) - 1; } }; rep(p, in.K) { i32 cur_city = pick_city(); i32 cur_time = START_MIN; plan[p].reserve(16); while (true) { i32 dst = pick_city_except(cur_city); i32 dep = cur_time; i32 arr = dep + GlobalData::duration[cur_city][dst]; if (arr > END_MIN) break; plan[p].push_back( Flight{(i16)cur_city, (i16)dst, (i16)dep, (i16)arr}); cur_city = dst; cur_time = arr; } } return plan; } static inline u64 city_weight(i16 c) { u64 w = (u64)in.w[c]; #if WEIGHT_STYLE == 1 return (u64)llround(sqrt((long double)w)); #elif WEIGHT_STYLE == 2 return w * w; #else return w; #endif } State() { planes = build_random_plan(); } score_t calc_score_full() { auto &ci = GlobalData::ci_flights_buf; ci.clear(); for (auto const &pl : planes) { for (auto const &f : pl) ci.push_back(f); } sort_desc(ci); compute_all_dp_flat_sorted(ci, GlobalData::dp_ci_flat_buf); __int128 v_ci = 0; for (auto const &term : GlobalData::eval_terms) { if (GlobalData::dp_sq_flat[term.idx] < GlobalData::dp_ci_flat_buf[term.idx]) { v_ci += (__int128)term.gain; } } if (GlobalData::total_gain_all == 0) return 0; return (score_t)((1000000 * v_ci) / GlobalData::total_gain_all); } bool apply_destroy_reconnect(RollbackInfo &rb) { if (in.K <= 0) return false; rb.plane_id = (i32)rng.random32(in.K); vector &fl = planes[rb.plane_id]; const int n = (int)fl.size(); if (n <= 1) return false; // 部分破壊: [l, r] を壊す(rは末尾でも可) int l = (int)rng.random32((u32)n); int max_len = rng.randomRange32(2, 5); int r = l + (int)rng.random32((u32)min(max_len, n - l)); // r <= n-1 i16 cur_city = fl[l].from; i16 cur_time = (l == 0) ? (i16)START_MIN : fl[l - 1].arr; bool has_suffix = (r + 1 < n); i16 target_city = -1; i16 target_dep = -1; if (has_suffix) { target_city = fl[r + 1].from; target_dep = fl[r + 1].dep; } rb.cut_pos = l; rb.tail.assign(fl.begin() + l, fl.end()); fl.resize(l); // 破壊区間のみを再接続 // suffixがある場合: // - 寄り道候補は「その便を経由しても target_dep までに target_city // に戻れる」もののみ // - 寄り道不能なら、追加せず直接 target_city に戻る // suffixがない場合: 当日中に到着可能な候補のみ // 候補バッファ: (累積重み, 都市番号) を1回のスキャンで構築 static vector> cands; int guard = 0; while (guard < in.N * 8) { if (has_suffix && cur_city == target_city) break; guard++; if (has_suffix) { i16 direct_arr = cur_time + GlobalData::duration[cur_city][target_city]; bool can_direct = (direct_arr <= target_dep); cands.clear(); u64 acc = 0; for (i16 dst = 0; dst < in.N; dst++) { if (dst == cur_city || dst == target_city) continue; i16 arr = cur_time + GlobalData::duration[cur_city][dst]; if (arr > target_dep) continue; i16 back_arr = arr + GlobalData::duration[dst][target_city]; if (back_arr > target_dep) continue; acc += city_weight(dst); cands.push_back({acc, dst}); } if (acc == 0) { if (!can_direct) break; // 追加せず直接戻る fl.push_back( Flight{cur_city, target_city, cur_time, direct_arr}); cur_city = target_city; cur_time = direct_arr; continue; } u64 rsel = rng.random64(acc); auto it = lower_bound(cands.begin(), cands.end(), make_pair(rsel + 1, (i16)0)); if (it == cands.end()) break; i16 dst = it->second; i16 dep = cur_time; i16 arr = dep + GlobalData::duration[cur_city][dst]; fl.push_back(Flight{cur_city, dst, dep, arr}); cur_city = dst; cur_time = arr; } else { cands.clear(); u64 acc = 0; for (i16 dst = 0; dst < in.N; dst++) { if (dst == cur_city) continue; i16 arr = cur_time + GlobalData::duration[cur_city][dst]; if (arr > END_MIN) continue; acc += city_weight(dst); cands.push_back({acc, dst}); } if (acc == 0) break; u64 rsel = rng.random64(acc); auto it = lower_bound(cands.begin(), cands.end(), make_pair(rsel + 1, (i16)0)); if (it == cands.end()) break; i16 dst = it->second; i16 dep = cur_time; i16 arr = dep + GlobalData::duration[cur_city][dst]; fl.push_back(Flight{cur_city, dst, dep, arr}); cur_city = dst; cur_time = arr; } } // suffixがあるのに接続できなければロールバック if (has_suffix && (cur_city != target_city || cur_time > target_dep)) { fl.resize(rb.cut_pos); fl.insert(fl.end(), rb.tail.begin(), rb.tail.end()); return false; } // 壊していないsuffixを復元 if (has_suffix) { int suffix_begin_in_tail = r + 1 - l; fl.insert(fl.end(), rb.tail.begin() + suffix_begin_in_tail, rb.tail.end()); } return true; } }; // ========================================== // SA // ========================================== // 多点スタートに対応するためグローバル変数で持つ struct MoveStats { i64 try_cnt = 0; i64 evaluated_cnt = 0; i64 accept_cnt = 0; i64 improve_cnt = 0; i64 update_best_cnt = 0; } stats[NUM_MOVES]; i64 iter_count = 0; class SA { const double SA_TIME_LIMIT; const State &init_state; AliasWeightedSampler &move_selector; inline void print_stats(i64 iter, f64 time, f64 temp, score_t score, score_t best) { if constexpr (!DEBUG) return; cerr << "VISUALIZE: {"; cerr << "\"iter\": " << iter << ", "; cerr << "\"time\": " << fixed << setprecision(3) << time << ", "; cerr << "\"temp\": " << fixed << setprecision(3) << temp << ", "; cerr << "\"score\": " << score << ", "; cerr << "\"best_score\": " << best << ", "; cerr << "\"moves\": ["; for (int i = 0; i < NUM_MOVES; i++) { if (i > 0) cerr << ", "; cerr << "{\"name\": \"" << MOVE_NAMES[i] << "\", \"tried\": " << stats[i].try_cnt << ", \"evaluated\": " << stats[i].evaluated_cnt << ", \"accepted\": " << stats[i].accept_cnt << ", \"improved\": " << stats[i].improve_cnt << ", \"updated_best\": " << stats[i].update_best_cnt << "}"; } cerr << "]}" << el; } public: SA(const double SA_TIME_LIMIT, const State &init_state, AliasWeightedSampler &move_selector) : SA_TIME_LIMIT(SA_TIME_LIMIT), init_state(init_state), move_selector(move_selector) {} State run() { State state = init_state; state.score = state.calc_score_full(); State best_state = state; f64 current_temp = START_TEMP; f64 time_start = timer.elapsed(); f64 time_now = time_start; State::RollbackInfo rb; // ループ外で確保してtailのメモリを再利用 while (true) { iter_count++; // 時間チェック if ((iter_count & TIME_CHECK_INTERVAL) == 0) { time_now = timer.elapsed(); if (time_now > SA_TIME_LIMIT) break; g_sa_progress = (time_now - time_start) / (SA_TIME_LIMIT - time_start); if constexpr (ALLOW_WORSE_MOVES) { f64 progress = (time_now - time_start) / (SA_TIME_LIMIT - time_start); current_temp = temp_decay(progress); } if (DEBUG && (iter_count % STATS_INTERVAL == 0)) { if constexpr (ALLOW_WORSE_MOVES) print_stats(iter_count, time_now, current_temp, state.score, best_state.score); else print_stats(iter_count, time_now, current_temp, state.score, state.score); } } // 採用閾値の事前計算 f64 threshold; if constexpr (!ALLOW_WORSE_MOVES) { threshold = 0.0; } else { f64 rand_val = rng.randomDoubleOpen01(); threshold = current_temp * log_fast(rand_val); } int move_type = move_selector.sample(); if constexpr (DEBUG) stats[move_type].try_cnt++; // 遷移適用と差分計算 score_t next_score = state.score; bool moved = false; moved = state.apply_destroy_reconnect(rb); if (!moved) { continue; } next_score = state.calc_score_full(); if constexpr (DEBUG) stats[move_type].evaluated_cnt++; // 採用判定(閾値ベース) f64 delta = (f64)(next_score - state.score); if constexpr (!MAXIMIZE) delta = -delta; bool accept = (delta >= threshold); if constexpr (DEBUG) { if (delta > 0) stats[move_type].improve_cnt++; } // 更新 or ロールバック if (accept) { if constexpr (DEBUG) stats[move_type].accept_cnt++; bool is_better; if constexpr (ALLOW_WORSE_MOVES) { is_better = MAXIMIZE ? (next_score > best_state.score) : (next_score < best_state.score); } else { is_better = MAXIMIZE ? (next_score > state.score) : (next_score < state.score); } state.score = next_score; if (is_better) { if constexpr (ALLOW_WORSE_MOVES) { best_state = state; } if constexpr (DEBUG) stats[move_type].update_best_cnt++; } } else { if (0 <= rb.plane_id && rb.plane_id < (int)state.planes.size()) { auto &fl = state.planes[rb.plane_id]; fl.resize(rb.cut_pos); fl.insert(fl.end(), rb.tail.begin(), rb.tail.end()); } } } if constexpr (!ALLOW_WORSE_MOVES) best_state = state; // DEBUGに関わらず出力 cerr << "Time: " << fixed << setprecision(3) << (timer.elapsed()) << " Iter: " << iter_count << " Score: " << best_state.score << el; return best_state; } }; // ========================================== // Solver // ========================================== class Solver { optional ans; static string format_time(i16 m) { i32 hh = m / 60; i32 mm = m % 60; return int_to_string(hh, 2) + ":" + int_to_string(mm, 2); } public: Solver() {} void solve() { using namespace GlobalData; targets.clear(); for (i16 t = 11 * 60; t <= 21 * 60; t += 30) targets.push_back(t); duration.assign(in.N, vector(in.N, 0)); pairs.clear(); rep(i, in.N) { rep(j, in.N) { if (i == j) continue; double dx = (double)in.x[i] - (double)in.x[j]; double dy = (double)in.y[i] - (double)in.y[j]; double dist = sqrt(dx * dx + dy * dy); duration[i][j] = State::flight_duration(dist); if (dist >= 0.25 * (double)in.R) { pairs.push_back( PairInfo{(i16)i, (i16)j, in.w[i] * in.w[j]}); } } } total_gain_all = 0; for (auto const &pi : pairs) { total_gain_all += (__int128)pi.gain * (__int128)targets.size(); } vector sq_sorted = in.sq_flights; State::sort_desc(sq_sorted); State::compute_all_dp_flat_sorted(sq_sorted, dp_sq_flat); eval_terms.clear(); eval_terms.reserve((size_t)pairs.size() * targets.size()); for (auto const &pi : pairs) { rep(k, targets.size()) { eval_terms.push_back( EvalTerm{State::dp_index(pi.dst, (i32)k, pi.src), pi.gain}); } } ci_flights_buf.clear(); ci_flights_buf.reserve((size_t)in.K * 32); dp_ci_flat_buf.clear(); dp_ci_flat_buf.reserve((size_t)in.N * targets.size() * in.N); AliasWeightedSampler move_selector; if constexpr (DEBUG) { assert((int)MOVE_NAMES.size() == NUM_MOVES); assert((int)MOVE_WEIGHTS.size() == NUM_MOVES); } move_selector.build(MOVE_WEIGHTS); // 初期解: 元のランダム生成を50回繰り返して最良を採用 State best_init; best_init.score = best_init.calc_score_full(); for (int i = 1; i < 50; i++) { State cand; cand.score = cand.calc_score_full(); if (cand.score > best_init.score) best_init = move(cand); } SA sa(TIME_LIMIT, best_init, move_selector); ans = sa.run(); } void print() { if (!ans.has_value()) return; State const &ret = ans.value(); rep(p, in.K) { cout << ret.planes[p].size() << el; for (auto const &f : ret.planes[p]) { cout << (f.from + 1) << ' ' << format_time(f.dep) << ' ' << (f.to + 1) << ' ' << format_time(f.arr) << el; } } } }; int main() { init_temp_table(); // 温度減衰テーブルの初期化 init_log_table(); // logテーブルの初期化 // init_gauss_table(); // ガウシアンテーブルの初期化 in.input(); Solver solver; solver.solve(); solver.print(); // 高速に終了する cout.flush(); cerr.flush(); _Exit(0); }