A - Stack to Match Pairs Editorial
by 
Kyo25
解法と考察など
今回コードを書いてくれたGeminiが解説も書いてくれました。
[AHC059] TSPバックボーンと利益ベースの包含探索
1. 解法の概要
本問題は、グリッド上のカードをペアで回収していく問題ですが、スタック (LIFO)の制約があるため、単純な巡回セールスマン問題 (TSP)よりも複雑です。しかし、基本的な戦略として以下の2段階のアプローチを採用することで高得点を得ることができました。
- バックボーンの構築 (Global Optimization): スタックの深さを1 (つまり、1ペアずつ確実に拾って捨てる)と仮定し、純粋なTSPとして最短経路を焼きなまし法 (SA)で求めます。
- 包含による短縮 (Local Optimization / Vacuum): 構築した経路（バックボーン）に従って移動する際、その道中で「ついでに拾えるペア」があれば、スタックに積み込んで（包含して）回収します。
2. アルゴリズム詳細
2.1. 問題のモデル化
カードのペア \(i\) (\(0 \le i < M\)) に対し、その位置を \(P_{i,1}, P_{i,2}\) とします。 ペア \(i\) を回収する手順は、以下の2通りの有向パスとみなせます。 * 順方向: \(P_{i,1} \to P_{i,2}\) (コスト: \(dist(Current, P_{i,1}) + dist(P_{i,1}, P_{i,2})\)) * 逆方向: \(P_{i,2} \to P_{i,1}\) (コスト: \(dist(Current, P_{i,2}) + dist(P_{i,2}, P_{i,1})\))
この「順序」と「方向」を決定する問題として捉えます。
2.2. 初期解の構築 (Clustering Strategy)
ランダムな初期解では収束に時間がかかるため、座標に基づくクラスタリングを行いました。 盤面を \(4 \times 4\) のグリッドに分割し、Z字やS字などの順序でクラスタを巡回する初期解を複数生成し、最もスコアが良いものを採用しました。
2.3. 焼きなまし法 (Simulated Annealing)
TSPとしての経路最適化を行います。評価関数は「バックボーンを順に辿った時の総移動距離」です。
近傍操作: 特に以下の 2-opt (Reverse) がスコア向上に大きく寄与しました。 * Reverse (2-opt): 区間 \([i, j]\) の訪問順序を反転させます。同時に、その区間内の各ペアの通過方向（\(dirs\)）も反転させることで、幾何学的な経路の整合性を保ち、交差（ねじれ）を解消します。 * Insert (Block Move): 連続する区間を別の場所に移動させます。 * Swap: 2つのペアの位置を交換します。 * Flip: ペアの通過方向 (\(P_1 \to P_2\) か \(P_2 \to P_1\) か)を反転させます。
高速化のため、距離計算は事前計算したテーブルを用い、また Reverse などの操作では差分計算（あるいは高速な \(O(M)\) 再計算）を行い、1.98秒ギリギリまで数千万回の試行を行っています。
2.4. 利益ベースの包含探索 (Gain-based Vacuum)
ここが本解法の肝です。SAで決定したバックボーン経路 \(A \to B \to C \dots\) に従って移動しますが、移動中に他のペアを「吸い込む（Vacuum）」判定を行います。
あるペア \(A\) の処理中（\(A_{start} \to A_{end}\)）に、まだ回収していないペア \(B\) を割り込ませるかどうかの判断基準として、純利益 (Profit) を導入しました。
\[ Profit(B) = Gain(B) - Detour(B) \]
- \(Gain(B)\): ペア \(B\) をバックボーンの本来の順序から削除した場合に短縮される距離（将来の得）。
- \(Detour(B)\): 今、現在地から \(B\) を拾って元の経路に戻るために増加する距離（現在の損）。
\(Profit(B) > - \epsilon\) （\(\epsilon\) は許容する微小な損）であれば、ペア \(B\) を割り込ませて回収します。これを再帰的に行うことで、\(A \supset B \supset C\) のような多重のマトリョーシカ構造（スタック活用）が自然に形成されます。
この戦略により、「近くにあるから拾う」だけでなく、「遠くにあっても、将来回収するのが大変な場所にあるなら、多少無理をしてでも今拾う」という高度な判断が可能になりました。
3. 実装上の工夫
- Xorshift128: 標準の乱数生成器より高速な乱数を使用。
- 動的閾値: Vacuum処理において、移動距離が長い場合は探索範囲（許容Detour）を広げる動的な閾値を採用。
- 計算量削減: 全点間距離およびペア内部距離を事前計算し、定数倍を高速化。
4. 考察と今後の課題
今回の解法は、「大域的なTSP」と「局所的なスタック最適化」を分離したハイブリッド手法でした。
- 分離のメリット: 探索空間を限定でき、短時間で良質な解に到達しやすい。
- 分離のデメリット: 「スタックに積むことを前提とした、TSP的には遠回りな経路」をSAが見つけにくい。
より上位を目指すには、SAの評価関数自体に「吸い込みによるコスト削減」を組み込む（Plan 2のようなアプローチ）か、ビームサーチ等で順序とスタック操作を同時に探索する方法が考えられます。しかし、計算資源が限られる中では、今回のような「強力なTSPソルバ + 賢い貪欲法」の組み合わせが非常にバランスの良い戦略であったと言えます。
// GCC 15.2.0 対応
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")
#pragma GCC target("avx2,bmi2,popcnt,lzcnt")
#include <bits/stdc++.h>
using namespace std;
// ============================================================================
// 超高速乱数生成器 (Xorshift128)
// ============================================================================
struct Xorshift128 {
uint64_t x = 123456789, y = 362436069, z = 521288629, w = 88675123;
inline uint64_t next() {
uint64_t t = x ^ (x << 11);
x = y; y = z; z = w;
return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
}
inline int next_int(int n) {
return (int)(next() % (uint64_t)n);
}
inline double next_double() {
return (next() >> 11) * (1.0 / 9007199254740992.0);
}
} rng;
// ============================================================================
// 定数・データ構造
// ============================================================================
const int N = 20;
const int M = N * N / 2; // 200 ペア
const int MAX_OPERATIONS = 16000;
struct Point {
int y, x;
};
// 距離テーブル: [ペアID * 2 + 向き][ペアID * 2 + 向き]
// 向き0: p1, 向き1: p2
int dist_table[400][400];
// ペア内の距離: internal_dist[g] = dist(p1, p2)
int internal_dist[200];
struct Group {
Point p1, p2;
};
// ============================================================================
// ソルバクラス
// ============================================================================
class Solver {
public:
int n;
vector<Group> groups;
vector<vector<Point>> card_positions;
Solver() : card_positions(N * N) {}
void read_input() {
cin >> n;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
int val; cin >> val;
card_positions[val].push_back({i, j});
}
}
// グループ生成
for (int val = 0; val < N * N; val++) {
if (!card_positions[val].empty()) {
groups.push_back({card_positions[val][0], card_positions[val][1]});
int g_idx = (int)groups.size() - 1;
Point p1 = groups.back().p1;
Point p2 = groups.back().p2;
internal_dist[g_idx] = abs(p1.y - p2.y) + abs(p1.x - p2.x);
}
}
// 全点間距離の事前計算
for (int i = 0; i < M; i++) {
Point ip[2] = {groups[i].p1, groups[i].p2};
for (int j = 0; j < M; j++) {
Point jp[2] = {groups[j].p1, groups[j].p2};
for (int di = 0; di < 2; di++) {
for (int dj = 0; dj < 2; dj++) {
int u = 2 * i + di;
int v = 2 * j + dj;
dist_table[u][v] = abs(ip[di].y - jp[dj].y) + abs(ip[di].x - jp[dj].x);
}
}
}
}
}
// ===========================
// 初期解構築 (Cluster Strategy)
// ===========================
void build_initial_solution(vector<int> &best_order, vector<int> &best_dirs,
int &best_init_score) {
vector<Point> center(M);
for (int g = 0; g < M; g++) {
center[g] = {
(groups[g].p1.y + groups[g].p2.y) / 2,
(groups[g].p1.x + groups[g].p2.x) / 2
};
}
const int G = 16;
array<vector<int>, G> base_clusters;
for (int g = 0; g < M; g++) {
int gr = center[g].y / 5;
int gc = center[g].x / 5;
if (gr < 0) gr = 0; if (gr > 3) gr = 3;
if (gc < 0) gc = 0; if (gc > 3) gc = 3;
int gid = gr * 4 + gc;
base_clusters[gid].push_back(g);
}
vector<array<int, G>> cluster_orders;
{ array<int, G> ord{}; iota(ord.begin(), ord.end(), 0); cluster_orders.push_back(ord); }
{
array<int, G> ord{}; int idx = 0;
for (int gr = 0; gr < 4; gr++) for (int gc = 0; gc < 4; gc++) ord[idx++] = gr * 4 + gc;
cluster_orders.push_back(ord);
}
{
array<int, G> ord{}; int idx = 0;
for (int gr = 0; gr < 4; gr++) {
if (gr % 2 == 0) for (int gc = 0; gc < 4; gc++) ord[idx++] = gr * 4 + gc;
else for (int gc = 3; gc >= 0; gc--) ord[idx++] = gr * 4 + gc;
}
cluster_orders.push_back(ord);
}
for (int rep = 0; rep < 2; rep++) {
array<int, G> ord{}; iota(ord.begin(), ord.end(), 0);
for (int i = G - 1; i > 0; --i) {
int j = rng.next_int(i + 1);
swap(ord[i], ord[j]);
}
cluster_orders.push_back(ord);
}
best_init_score = INT_MAX;
for (auto &cord : cluster_orders) {
vector<int> order(M);
vector<int> dirs(M);
vector<bool> used(M, false);
Point cur{0, 0};
int idx = 0;
for (int cid_k = 0; cid_k < G; cid_k++) {
int cid = cord[cid_k];
auto &vec = base_clusters[cid];
while (true) {
int best_g = -1, best_dist = INT_MAX, best_dir = 0;
for (int g : vec) {
if (used[g]) continue;
const Group &group = groups[g];
int d1 = abs(cur.y - group.p1.y) + abs(cur.x - group.p1.x);
int d2 = abs(cur.y - group.p2.y) + abs(cur.x - group.p2.x);
if (d1 < best_dist) { best_dist = d1; best_g = g; best_dir = 0; }
if (d2 < best_dist) { best_dist = d2; best_g = g; best_dir = 1; }
}
if (best_g == -1) break;
order[idx] = best_g;
dirs[idx] = best_dir;
used[best_g] = true;
idx++;
cur = (best_dir == 0) ? groups[best_g].p2 : groups[best_g].p1;
}
}
while (idx < M) {
int best_g = -1, best_dist = INT_MAX, best_dir = 0;
for (int g = 0; g < M; g++) {
if (used[g]) continue;
const Group &group = groups[g];
int d1 = abs(cur.y - group.p1.y) + abs(cur.x - group.p1.x);
int d2 = abs(cur.y - group.p2.y) + abs(cur.x - group.p2.x);
if (d1 < best_dist) { best_dist = d1; best_g = g; best_dir = 0; }
if (d2 < best_dist) { best_dist = d2; best_g = g; best_dir = 1; }
}
if (best_g == -1) break;
order[idx] = best_g;
dirs[idx] = best_dir;
used[best_g] = true;
idx++;
cur = (best_dir == 0) ? groups[best_g].p2 : groups[best_g].p1;
}
// 初期スコア計算
int score = 0;
int prev_node = -1;
Point start_p{0, 0};
for (int i = 0; i < M; i++) {
int g = order[i];
int d = dirs[i];
int entry_node = 2 * g + d;
int exit_node = 2 * g + (1 - d);
int dist_to_entry;
if (prev_node == -1) {
Point entry_p = (d == 0) ? groups[g].p1 : groups[g].p2;
dist_to_entry = abs(start_p.y - entry_p.y) + abs(start_p.x - entry_p.x);
} else {
dist_to_entry = dist_table[prev_node][entry_node];
}
score += dist_to_entry + internal_dist[g];
prev_node = exit_node;
}
if (score < best_init_score) {
best_init_score = score;
best_order = order;
best_dirs = dirs;
}
}
}
// --- 高速SA ---
void solve() {
vector<int> order(M), dirs(M);
int init_score;
build_initial_solution(order, dirs, init_score);
auto start_time = chrono::high_resolution_clock::now();
const double TIME_LIMIT = 1.98; // ギリギリまで
const double TEMP_START = 2500.0;
const double TEMP_END = 0.01;
int current_score = init_score;
vector<int> best_order = order;
vector<int> best_dirs = dirs;
int best_score = current_score;
int iter_count = 0;
double temp = TEMP_START;
// 作業用バッファ
vector<int> next_order; next_order.reserve(M);
vector<int> next_dirs; next_dirs.reserve(M);
while (true) {
iter_count++;
if ((iter_count & 1023) == 0) {
auto now = chrono::high_resolution_clock::now();
double elapsed = chrono::duration_cast<chrono::duration<double>>(now - start_time).count();
if (elapsed > TIME_LIMIT) break;
double progress = elapsed / TIME_LIMIT;
temp = TEMP_START * pow(TEMP_END / TEMP_START, progress);
}
int r = rng.next_int(100);
int delta = 0;
// ★改良点: 2-opt (Reverse) を50%の確率で投入
if (r < 50) {
// ================= Reverse (2-opt) =================
// 区間 [i, j] を反転
// 幾何学的整合性のため、区間内の各要素の向き(dirs)も反転する
int i = rng.next_int(M - 1);
int len = 1 + rng.next_int(min(M - 1 - i, 30));
int j = i + len;
int idx_i = i;
int idx_j = j;
int u_prev = (idx_i == 0) ? -1 : (2 * order[idx_i - 1] + (1 - dirs[idx_i - 1]));
int u_next = (idx_j == M - 1) ? -2 : (2 * order[idx_j + 1] + dirs[idx_j + 1]);
// Old connections
// prev -> i_in ... j_out -> next
int u_in = 2 * order[idx_i] + dirs[idx_i];
int v_out = 2 * order[idx_j] + (1 - dirs[idx_j]);
int loss = 0;
if (u_prev == -1) {
Point p = (dirs[idx_i] == 0) ? groups[order[idx_i]].p1 : groups[order[idx_i]].p2;
loss += abs(p.y) + abs(p.x);
} else loss += dist_table[u_prev][u_in];
if (u_next != -2) loss += dist_table[v_out][u_next];
// New connections (after reverse and flip)
// prev -> j_out ... i_in -> next
// (Note: j's old 'out' becomes new 'in' when reversed, etc.)
int v_in_new = v_out;
int u_out_new = u_in;
int gain = 0;
if (u_prev == -1) {
Point p = (dirs[idx_j] == 0) ? groups[order[idx_j]].p2 : groups[order[idx_j]].p1; // reversed entrance
gain += abs(p.y) + abs(p.x);
} else gain += dist_table[u_prev][v_in_new];
if (u_next != -2) gain += dist_table[u_out_new][u_next];
delta = gain - loss;
if (delta < 0 || rng.next_double() < exp(-delta / temp)) {
reverse(order.begin() + i, order.begin() + j + 1);
reverse(dirs.begin() + i, dirs.begin() + j + 1);
for(int k=i; k<=j; ++k) dirs[k] = 1 - dirs[k];
current_score += delta;
if (current_score < best_score) {
best_score = current_score;
best_order = order;
best_dirs = dirs;
}
}
} else if (r < 80) {
// ================= Insert (Block Move) =================
int len = 1 + rng.next_int(6);
if (len > M) len = 1;
int i = rng.next_int(M);
if (i + len > M) i = M - len;
int j = rng.next_int(M);
if (j >= i && j <= i + len) continue;
next_order.clear(); next_dirs.clear();
for (int k = 0; k < M; k++) {
if (k >= i && k < i + len) continue;
next_order.push_back(order[k]);
next_dirs.push_back(dirs[k]);
}
if (j > i) j -= len;
if (j < 0) j = 0;
if (j > (int)next_order.size()) j = (int)next_order.size();
vector<int> block_o(order.begin() + i, order.begin() + i + len);
vector<int> block_d(dirs.begin() + i, dirs.begin() + i + len);
next_order.insert(next_order.begin() + j, block_o.begin(), block_o.end());
next_dirs.insert(next_dirs.begin() + j, block_d.begin(), block_d.end());
// Full Recalculation (Safe & Fast enough)
int next_score = 0;
int p_node = -1;
Point s{0, 0};
for (int k = 0; k < M; k++) {
int g = next_order[k], d = next_dirs[k];
int u_in = 2*g+d, u_out=2*g+(1-d);
if (p_node == -1) {
Point ep = (d==0)?groups[g].p1:groups[g].p2;
next_score += abs(s.y-ep.y)+abs(s.x-ep.x);
} else next_score += dist_table[p_node][u_in];
next_score += internal_dist[g];
p_node = u_out;
}
delta = next_score - current_score;
if (delta < 0 || rng.next_double() < exp(-delta / temp)) {
order = move(next_order);
dirs = move(next_dirs);
current_score = next_score;
if (current_score < best_score) {
best_score = current_score;
best_order = order;
best_dirs = dirs;
}
}
} else {
// ================= Swap / Flip =================
int i = rng.next_int(M);
int j = rng.next_int(M);
if (i == j) { // Flip
dirs[i] = 1 - dirs[i];
} else { // Swap
swap(order[i], order[j]);
swap(dirs[i], dirs[j]);
}
// Full Recalc
int next_score = 0;
int p_node = -1;
Point s{0, 0};
for (int k = 0; k < M; k++) {
int g = order[k], d = dirs[k];
int u_in = 2*g+d, u_out=2*g+(1-d);
if (p_node == -1) {
Point ep = (d==0)?groups[g].p1:groups[g].p2;
next_score += abs(s.y-ep.y)+abs(s.x-ep.x);
} else next_score += dist_table[p_node][u_in];
next_score += internal_dist[g];
p_node = u_out;
}
delta = next_score - current_score;
if (delta < 0 || rng.next_double() < exp(-delta / temp)) {
current_score = next_score;
if (current_score < best_score) {
best_score = current_score;
best_order = order;
best_dirs = dirs;
}
} else {
if (i == j) dirs[i] = 1 - dirs[i];
else { swap(order[i], order[j]); swap(dirs[i], dirs[j]); }
}
}
}
// --- Vacuum (Gain-based) ---
generate_output_with_gain(best_order, best_dirs);
}
// ===========================
// Vacuum 部分 (Gain-based)
// ===========================
vector<char> operations;
vector<int> removal_gain; // 各ペアをスキップしたときの短縮距離
void calc_gains(const vector<int>& order, const vector<int>& dirs) {
removal_gain.assign(M * 2, 0); // safe size
// removal_gain[g] = (Original Connect Cost) - (Shortcut Connect Cost)
// Original: prev -> g_in ... g_out -> next
// Shortcut: prev -> next
for (int i = 0; i < M; i++) {
int g = order[i];
int d = dirs[i];
int u_in = 2 * g + d;
int u_out = 2 * g + (1 - d);
int u_prev = (i == 0) ? -1 : (2 * order[i-1] + (1 - dirs[i-1]));
int u_next = (i == M - 1) ? -2 : (2 * order[i+1] + dirs[i+1]);
int orig_cost = 0;
if (u_prev == -1) {
Point p = (d==0)?groups[g].p1:groups[g].p2;
orig_cost += abs(p.y)+abs(p.x);
} else orig_cost += dist_table[u_prev][u_in];
if (u_next != -2) orig_cost += dist_table[u_out][u_next];
int short_cost = 0;
if (u_prev == -1) {
if (u_next != -2) {
int ng = u_next/2, nd = u_next%2;
Point p = (nd==0)?groups[ng].p1:groups[ng].p2;
short_cost += abs(p.y)+abs(p.x);
}
} else {
if (u_next != -2) short_cost += dist_table[u_prev][u_next];
}
// Gain cannot be negative (triangle inequality) but for safety
removal_gain[g] = max(0, orig_cost - short_cost);
}
}
void add_move(Point from, Point to) {
int dy = to.y - from.y;
int dx = to.x - from.x;
if (dy > 0) for (int k = 0; k < dy; ++k) operations.push_back('D');
else for (int k = 0; k < -dy; ++k) operations.push_back('U');
if (dx > 0) for (int k = 0; k < dx; ++k) operations.push_back('R');
else for (int k = 0; k < -dx; ++k) operations.push_back('L');
}
void process_recursive(int g_idx, int dir_idx, Point& current_pos,
unordered_set<int>& remaining_set) {
const Group& gA = groups[g_idx];
Point A_start = (dir_idx == 0) ? gA.p1 : gA.p2;
Point A_end = (dir_idx == 0) ? gA.p2 : gA.p1;
add_move(current_pos, A_start);
operations.push_back('Z');
current_pos = A_start;
while (true) {
int best_B = -1;
int best_B_dir = 0;
int max_profit = -1e9; // Profit = Gain - Detour
int base_dist = abs(current_pos.y - A_end.y) + abs(current_pos.x - A_end.x);
for (int cand_idx : remaining_set) {
int gain = removal_gain[cand_idx];
const Group& gB = groups[cand_idx];
for (int d = 0; d < 2; ++d) {
Point B_in = (d == 0) ? gB.p1 : gB.p2;
Point B_out = (d == 0) ? gB.p2 : gB.p1;
int d_in = abs(current_pos.y - B_in.y) + abs(current_pos.x - B_in.x);
int d_out = abs(B_out.y - A_end.y) + abs(B_out.x - A_end.x);
int detour = (d_in + d_out) - base_dist;
// Profit: 将来浮く距離 - 今増える距離
int profit = gain - detour;
// Profit > -2: 微小な損は許容 (スタック活用の機会損失を防ぐ)
// かつ、detour自体が極端に大きくないこと (リスク管理)
if (profit > -2 && detour < 20 && profit > max_profit) {
max_profit = profit;
best_B = cand_idx;
best_B_dir = d;
}
}
}
if (best_B != -1) {
remaining_set.erase(best_B);
process_recursive(best_B, best_B_dir, current_pos, remaining_set);
} else {
break;
}
}
add_move(current_pos, A_end);
operations.push_back('Z');
current_pos = A_end;
}
void generate_output_with_gain(const vector<int>& best_order, const vector<int>& best_dirs) {
calc_gains(best_order, best_dirs); // Gain事前計算
operations.clear();
operations.reserve(MAX_OPERATIONS);
unordered_set<int> remaining_set;
remaining_set.reserve(M * 2);
for (int g : best_order) remaining_set.insert(g);
Point robot_pos{0, 0};
for (int i = 0; i < M; i++) {
int g = best_order[i];
if (remaining_set.count(g)) {
remaining_set.erase(g);
process_recursive(g, best_dirs[i], robot_pos, remaining_set);
}
}
if ((int)operations.size() > MAX_OPERATIONS) operations.resize(MAX_OPERATIONS);
for (char op : operations) cout << op << '\n';
}
};
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
Solver solver;
solver.read_input();
solver.solve();
return 0;
}
私からすると何をやっているか全く理解できませんでしたが、Visualizerやスコアを観察すると改善が進んでいくことがわかりました。次回は黄色パフォーマンスが出せるように邁進したいと思います！
posted:
last update: