eijirouの競プロ参加記

主に AtCoder Heuristic Contest の解説記事を書く予定です。

AHC030 参加記

公式ビジュアライザ (Seed=19, Cost=100.837210)

THIRD プログラミングコンテスト2023(AtCoder Heuristic Contest 030) お疲れ様でした。

システムテストの得点が 2,050,795,126,975 で 3 位でした。

最終順位表

問題概要

クエリを通して、石油がある場所を特定する問題です。

油田の個数と形は入力で与えられますが、各油田がどこにあるかは与えられません。

詳しくは公式の問題文を参考にしてください。

atcoder.jp

方針

毎ターン次のことを行いました。

  1. 焼きなまし法を用いて、油田配置の中で対数尤度が大きいものを 32 個ほど得る。
  2. 得られた油田配置の中に実際の配置が含まれていると仮定し、各油田配置の事後確率を求める。
  3. 得られた事後分布を利用して占うマスを選ぶ。

他の上位者と大まかな流れは一緒だと思います。焼きなまし法の部分が私は得意だったようなので、この記事では焼きなまし法で油田配置を推定する部分を中心に説明します。

ベイズ推定

占う場所を選ぶときに油田配置の事後分布が欲しくなるので、油田配置の事後分布を求めることを目標とします。

占いの履歴を  Q とし、各油田配置を  x と表記します。

 Q が得られたときに実際の油田配置が  x である事後確率  P(x | Q) を求めたいです。入力生成方法を読むと事前確率  P(x) が全て等しいことがわかります。ベイズの定理を踏まえると次の式が成り立ちます。

 P(x | Q) = \frac{P(Q | x)}{\sum_{x}{P(Q | x)}}

 P(x | Q) を求める代わりに  P(Q | x) を求めればよいことがわかりました。

 P(Q | x) = \prod_{q \in Q}{P(q | x)} です。

 P(q | x)正規分布の一部を積分した値であり、正規分布の累積分布関数が  \frac{1}{2} \left( 1 + \mathrm{erf} \frac{x - \mu}{\sqrt{2 \sigma^{2}}} \right) であることから計算できます。

油田の個数  M がある程度大きいとき、全ての油田配置  x について  P(Q | x) を計算する時間がありません。そこで、 P(Q | x) が比較的大きいものをいくつか集めた集合を  X' とし、 X' に実際の油田配置が含まれていると仮定して

 P(x | Q) = \frac{P(Q | x)}{\sum_{x \in X'}{P(Q | x)}}

で分布を近似することにしました。

 |X'| と過去の占いの情報量が十分大きければ、高い精度で近似できると思います。

焼きなまし法

 X' を求める手法はいくつか考えられ、大きく 2 つに分類できると思います。1 種類目は枝刈り全探索、MCTS、ビームサーチなどの全探索を改善したもので、2 種類目は山登り法や焼きなまし法などの局所探索です。

結論として、今回は焼きなまし法を使用しました。 M が大きいときに探索空間が広すぎて全探索系の手法が弱そうで、また、局所性がある近傍を設定して焼きなますことができそうだったからです。

焼きなまし法の評価値として、対数尤度  \log_{2}{P(Q | x)} を使用しました。

評価値が高い状態を 32 個ほど優先度付きキューで管理し、最後に保持していたものを  X' としました。同じ状態やクエリで推測して間違えた配置は省くようにします。

近傍は 3 種類です。

近傍 1: 1 マス移動(選択率: 9 / 16)

ランダムに 1 つの油田を選び、前後左右のランダムな方向に 1 マスずらします。

右への移動で揃う例

近傍 2: ランダムな場所に移動(選択率: 3 / 16)

ランダムに 1 つの油田を選び、ランダムに選んだ場所に移動させます。

ランダムな場所への移動で揃う例

近傍 3: 2 つの油田の入れ替え(選択率: 4 / 16)

ランダムに 2 つの油田を選び、2 つの油田の位置を入れ替えます。

石油埋蔵量が変化するマスの数が最も少なくなるような入れ替えのみを近傍として設定しました。複数通りある場合はそのいずれかをランダムに選びます。入れ替え方は入力を読み込んだ直後に事前計算します。

油田の入れ替えで揃う例

状態の持ち方

状態は、各油田の左上の座標と、各占いに対する石油埋蔵量を保持するようにしました。 N \times N の二次元グリッドは使いません。

油田が  M 個あって、それぞれの置き方はたかだか  N^{2} 通りしかありません。すなわち、1 つの油田を選んで配置する場合の数は  O(N^{2} M) です。そこで、占いの結果を得るたびに  O(N^{2} M) 通りそれぞれに対して占った場所の個数を計算するようにしました。この処理によって石油埋蔵量を  O(|Q|) で更新できます。

また、占いの結果を得るたびに、実際の石油埋蔵量が  k だったときにその占いの観測値が得られる確率も計算しておきます。焼きなまし法の最中に誤差関数を呼ぶ必要がなくなり、尤度の計算を  O(|Q|) で実現できます。

実装では石油埋蔵量が変わる占いのみを取り出して更新していますが、石油埋蔵量が変わる占いの割合が高そうなので高速化されているかはわかりません。

多点スタート

8 つの状態を並行して焼きなまし、少しずつ状態を減らすようにしました。

状態数の変化

状態を減らすとき、基本的にはその時点での評価値が最も低いものが削除されます。ただし、評価値の悪化を許容した直後に削除されるようなケースもあるため、評価値の指数移動平均をとったり、遷移の受容確率を評価したりしました。改善の寄与は小さいです。

初期状態

焼きなまし法の 1 ターン目の初期状態はランダムに生成しますが、2 ターン目以降の初期状態として、1 つ前のターンで得られた油田配置の中から尤度が高いものを選ぶようにしました。

初期解の尤度が高いと焼きなまし法の収束が速くなり、探索性能が改善されるようでした。

8 つの状態から 32 個の状態を生成し、次のターンは 32 個の状態から評価値が高い 8 つの状態を選ぶため、焼きなまし法を遷移としたビームサーチと捉えることもできます。

温度の自動調節

基本的には、初期温度を 6.0、最終温度を 0.6 にして温度を指数的に減少させます。

ただし、温度を上記のように固定すると、占いの回数が増えるにつれて焼きなまし法の遷移採択率が低くなるようだったので、採択率が 0.04 を下回ったときに初期温度と最終温度を 1.05 倍しました。採択率が 0.04 以上になるように温度を自動調節しているつもりです。難しいケースだと初期温度が 10.0 程度まで上がるようでした。

温度の自動調節

近傍における油田の選び方

近傍の説明でランダムに油田を選ぶと書きましたが、油田を等しい確率で選ぶのではなく、過去の遷移における採択率が高かったものを優先して選びました。

改善の寄与は小さいです。

対数尤度を整数で表現する

浮動小数点演算よりも整数演算のほうが速いです。そこで、対数尤度に  2^{32} をかけたものを 64 ビット整数で管理しました。焼きなまし法の中で浮動小数点数は扱いません。

速度を測っていないので確実なことはわからないですが、あまり速くなっていない気がします。

占うマスの決め方

最尤油田配置とそうでない油田配置を分離することを考えました。

焼きなまし法で求めた擬似的な事後分布から、各マスの石油埋蔵量の期待値  E \left[ v(i, j) \right] を求め、最尤油田配置における  v(i, j) E \left[ v(i, j) \right] 以上なら占うマスとして暫定的に採用します。

他の候補同士も分離したかったので、適当に評価関数を作って山登り法で占うマスの集合を改善しているつもりですが、評価関数の質が悪かったからか、思うようには改善されませんでした。

感想戦エントロピー最小化(占いとしては相互情報量の最大化)をすればよいことを知って納得しました。情報理論が身についていないのが敗因です。

序盤の占い方

最初の数ターンはルールベースで占うマスを決めました。事後分布は使いません。

 M が小さいときには損だったかもしれません。

時間管理

必要なターン数を雑に見積もったうえで、序盤のターンほど時間が長くなるように調整しました。

非常に難しいケースでは各ターンが使う時間を等しくしました。

提出コード

atcoder.jp

最後に

天才的なアイディアが必要というよりは、理論とヒューリスティック最適化の基礎的な理解と応用力が問われる問題だと思っていて、かなり自分好みの問題でした。

優勝できなかったのは悔しいですが、情報理論の大切さや MCMC との関連性など、多くのことを学べてよかったです。

コンテストを主催してくださった THIRD 様、writer の wata さん、コンテスト参加者の皆様、最後まで読んでくださった読者の皆様、ありがとうございました。

差分更新ビームサーチライブラリの実装 (C++)

この記事について

  • 差分更新型ビームサーチライブラリの実装例について説明します。
  • 差分更新型のビームサーチについては 高速なビームサーチが欲しい!!! などで既に解説されているため、被る部分については詳しく説明しません。
  • この記事に書いたソースコードプログラミングコンテストで自由に使っていただいて構いませんが、当ブログは損害などに対する責任を負いかねますのでご了承ください。

前提知識

差分更新型のビームサーチとは

木上のビームサーチ、Euler Tour ビームサーチとも呼ばれています1

ビームサーチは、幅優先探索に枝刈りを取り入れた手法です。探索木における深さ  d の頂点集合から深さ  d + 1 の頂点集合を生成し、その中から評価値が高い上位  W 個を選択します。 W のことをビーム幅と呼びます。

ビームサーチの例( W = 3、赤色が探索中のノード、橙色が採用されたノード、灰色が不採用のノード)

愚直なビームサーチでは、頂点ごとに独立した状態を作成します。実装は比較的楽だと思いますが、状態や履歴をコピーする必要が出てくるため、実行速度が遅くなりやすいです2

そこで、探索木を明示的に作成し、Euler Tour の順序で1つの状態を更新するようにしたものが差分更新型のビームサーチです。探索木の葉を訪れたときに新しい葉の候補を生成します。

Euler Tour の例

状態遷移の履歴は探索木から復元できるため、状態のコピーだけでなく履歴のコピーも省略できます。

深さが  d の1つの葉を探索するのに最悪で  2 d 回の遷移を必要としますが、実際には多くの葉で近い先祖を共有するため、遷移回数は平均して  2 d 回よりもかなり少なくなります。

差分更新型ビームサーチは愚直なビームサーチよりも高速に動作する場合が多く3、実装がやや複雑なため、ライブラリを作成しました。

実装

ライブラリ全体を折りたたみに記載しました。コードを読む場合は下の方にある beam_search 関数を最初に読むと分かりやすいかもしれません。


ビームサーチライブラリ

#include <bits/stdc++.h>
#include <atcoder/all>

using namespace std;

namespace beam_search {

// メモリの再利用を行いつつ集合を管理するクラス
template<class T>
class ObjectPool {
    public:
        // 配列と同じようにアクセスできる
        T& operator[](int i) {
            return data_[i];
        }

        // 配列の長さを変更せずにメモリを確保する
        void reserve(size_t capacity) {
            data_.reserve(capacity);
        }

        // 要素を追加し、追加されたインデックスを返す
        int push(const T& x) {
            if (garbage_.empty()) {
                data_.push_back(x);
                return data_.size() - 1;
            } else {
                int i = garbage_.top();
                garbage_.pop();
                data_[i] = x;
                return i;
            }
        }

        // 要素を(見かけ上)削除する
        void pop(int i) {
            garbage_.push(i);
        }

        // 使用した最大のインデックス(+1)を得る
        // この値より少し大きい値をreserveすることでメモリの再割り当てがなくなる
        size_t size() {
            return data_.size();
        }

    private:
        vector<T> data_;
        stack<int> garbage_;
};

// 連想配列
// Keyにハッシュ関数を適用しない
// open addressing with linear probing
// unordered_mapよりも速い
// nは格納する要素数よりも4~16倍ほど大きくする
template <class Key, class T>
struct HashMap {
    public:
        explicit HashMap(uint32_t n) {
            n_ = n;
            valid_.resize(n_, false);
            data_.resize(n_);
        }

        // 戻り値
        // - 存在するならtrue、存在しないならfalse
        // - index
        pair<bool,int> get_index(Key key) const {
            Key i = key % n_;
            while (valid_[i]) {
                if (data_[i].first == key) {
                    return {true, i};
                }
                if (++i == n_) {
                    i = 0;
                }
            }
            return {false, i};
        }

        // 指定したindexにkeyとvalueを格納する
        void set(int i, Key key, T value) {
            valid_[i] = true;
            data_[i] = {key, value};
        }

        // 指定したindexのvalueを返す
        T get(int i) const {
            assert(valid_[i]);
            return data_[i].second;
        }

        void clear() {
            fill(valid_.begin(), valid_.end(), false);
        }

    private:
        uint32_t n_;
        vector<bool> valid_;
        vector<pair<Key,T>> data_;
};

using Hash = uint32_t; // TODO

// 状態遷移を行うために必要な情報
// メモリ使用量をできるだけ小さくしてください
struct Action {
    // TODO

    Action() {
        // TODO
    }
};

using Cost = int; // TODO

// 状態のコストを評価するための構造体
// メモリ使用量をできるだけ小さくしてください
struct Evaluator {
    // TODO

    Evaluator() {
        // TODO
    }

    // 低いほどよい
    Cost evaluate() const {
        // TODO
    }
};

// 展開するノードの候補を表す構造体
struct Candidate {
    Action action;
    Evaluator evaluator;
    Hash hash;
    int parent;
    Cost cost;

    Candidate(Action action, Evaluator evaluator, Hash hash, int parent, Cost cost) :
        action(action),
        evaluator(evaluator),
        hash(hash),
        parent(parent),
        cost(cost) {}
};

// ビームサーチの設定
struct Config {
    int max_turn;
    size_t beam_width;
    size_t nodes_capacity;
    uint32_t hash_map_capacity;
};

// 削除可能な優先度付きキュー
using MaxSegtree = atcoder::segtree<
    pair<Cost,int>,
    [](pair<Cost,int> a, pair<Cost,int> b){
        if (a.first >= b.first) {
            return a;
        } else {
            return b;
        }
    },
    []() { return make_pair(numeric_limits<Cost>::min(), -1); }
>;

// ノードの候補から実際に追加するものを選ぶクラス
// ビーム幅の個数だけ、評価がよいものを選ぶ
// ハッシュ値が一致したものについては、評価がよいほうのみを残す
class Selector {
    public:
        explicit Selector(const Config& config) :
            hash_to_index_(config.hash_map_capacity)
        {
            beam_width = config.beam_width;
            candidates_.reserve(beam_width);
            full_ = false;
            st_original_.resize(beam_width);
        }

        // 候補を追加する
        // ターン数最小化型の問題で、candidateによって実行可能解が得られる場合にのみ finished = true とする
        // ビーム幅分の候補をCandidateを追加したときにsegment treeを構築する
        void push(Action action, const Evaluator& evaluator, Hash hash, int parent, bool finished) {
            Cost cost = evaluator.evaluate();
            if (finished) {
                finished_candidates_.emplace_back(Candidate(action, evaluator, hash, parent, cost));
                return;
            }
            if (full_ && cost >= st_.all_prod().first) {
                // 保持しているどの候補よりもコストが小さくないとき
                return;
            }
            auto [valid, i] = hash_to_index_.get_index(hash);

            if (valid) {
                int j = hash_to_index_.get(i);
                if (hash == candidates_[j].hash) {
                    // ハッシュ値が等しいものが存在しているとき
                    if (cost < candidates_[j].cost) {
                        // 更新する場合
                        candidates_[j] = Candidate(action, evaluator, hash, parent, cost);
                        if (full_) {
                            st_.set(j, {cost, j});
                        }
                    }
                    return;
                }
            }
            if (full_) {
                // segment treeが構築されている場合
                int j = st_.all_prod().second;
                hash_to_index_.set(i, hash, j);
                candidates_[j] = Candidate(action, evaluator, hash, parent, cost);
                st_.set(j, {cost, j});
            } else {
                // segment treeが構築されていない場合
                hash_to_index_.set(i, hash, candidates_.size());
                candidates_.emplace_back(Candidate(action, evaluator, hash, parent, cost));

                if (candidates_.size() == beam_width) {
                    // 保持している候補がビーム幅分になったとき
                    construct_segment_tree();
                }
            }
        }

        // 選んだ候補を返す
        const vector<Candidate>& select() const {
            return candidates_;
        }

        // 実行可能解が見つかったか
        bool have_finished() const {
            return !finished_candidates_.empty();
        }

        // 実行可能解に到達する「候補」を返す
        vector<Candidate> get_finished_candidates() const {
            return finished_candidates_;
        }

        void clear() {
            candidates_.clear();
            hash_to_index_.clear();
            full_ = false;
        }

    private:
        size_t beam_width;
        vector<Candidate> candidates_;
        HashMap<Hash,int> hash_to_index_;
        bool full_;
        vector<pair<Cost,int>> st_original_;
        MaxSegtree st_;
        vector<Candidate> finished_candidates_;

        void construct_segment_tree() {
            full_ = true;
            for (size_t i = 0; i < beam_width; ++i) {
                st_original_[i] = {candidates_[i].cost, i};
            }
            st_ = MaxSegtree(st_original_);
        }
};

// 深さ優先探索に沿って更新する情報をまとめたクラス
class State {
    public:
        explicit State(/* const Input& input */) {
            // TODO
        }

        // 次の状態候補を全てselectorに追加する
        // 引数
        //   evaluator : 今の評価器
        //   hash      : 今のハッシュ値
        //   parent    : 今のノードID(次のノードにとって親となる)
        void expand(const Evaluator& evaluator, Hash hash, int parent, Selector& selector) {
            // TODO
        }

        // actionを実行して次の状態に遷移する
        void move_forward(Action action) {
            // TODO
        }

        // actionを実行する前の状態に遷移する
        // 今の状態は、親からactionを実行して遷移した状態である
        void move_backward(Action action) {
            // TODO
        }

    private:
        // TODO
};

// 探索木(二重連鎖木)のノード
struct Node {
    Action action;
    Evaluator evaluator;
    Hash hash;
    int parent, child, left, right;

    // 根のコンストラクタ
    Node(Action action, const Evaluator& evaluator, Hash hash) :
        action(action),
        evaluator(evaluator),
        hash(hash),
        parent(-1),
        child(-1),
        left(-1),
        right(-1) {}

    // 通常のコンストラクタ
    Node(const Candidate& candidate, int right) :
        action(candidate.action),
        evaluator(candidate.evaluator),
        hash(candidate.hash),
        parent(candidate.parent),
        child(-1),
        left(-1),
        right(right) {}
};

// 二重連鎖木に対する操作をまとめたクラス
class Tree {
    public:
        explicit Tree(const State& state, size_t nodes_capacity, const Node& root) :
            state_(state)
        {
            nodes_.reserve(nodes_capacity);
            root_ = nodes_.push(root);
        }

        // 状態を更新しながら深さ優先探索を行い、次のノードの候補を全てselectorに追加する
        void dfs(Selector& selector) {
            update_root();

            int v = root_;
            while (true) {
                v = move_to_leaf(v);
                state_.expand(nodes_[v].evaluator, nodes_[v].hash, v, selector);
                v = move_to_ancestor(v);
                if (v == root_) {
                    break;
                }
                v = move_to_right(v);
            }
        }

        // 根からノードvまでのパスを取得する
        vector<Action> get_path(int v) {
            // cerr << nodes_.size() << endl;

            vector<Action> path;
            while (nodes_[v].parent != -1) {
                path.push_back(nodes_[v].action);
                v = nodes_[v].parent;
            }
            reverse(path.begin(), path.end());
            return path;
        }

        // 新しいノードを追加する
        int add_leaf(const Candidate& candidate) {
            int parent = candidate.parent;
            int sibling = nodes_[parent].child;
            int v = nodes_.push(Node(candidate, sibling));

            nodes_[parent].child = v;

            if (sibling != -1) {
                nodes_[sibling].left = v;
            }
            return v;
        }

        // ノードvに子がいなかった場合、vと不要な先祖を削除する
        void remove_if_leaf(int v) {
            if (nodes_[v].child == -1) {
                remove_leaf(v);
            }
        }

        // 最も評価がよいノードを返す
        int get_best_leaf(const vector<int>& last_nodes) {
            assert(!last_nodes.empty());
            int ret = last_nodes[0];
            for (int v : last_nodes) {
                if (nodes_[v].evaluator.evaluate() < nodes_[ret].evaluator.evaluate()) {
                    ret = v;
                }
            }
            return ret;
        }

    private:
        State state_;
        ObjectPool<Node> nodes_;
        int root_;

        // 根から一本道の部分は往復しないようにする
        void update_root() {
            int child = nodes_[root_].child;
            while (child != -1 && nodes_[child].right == -1) {
                root_ = child;
                state_.move_forward(nodes_[child].action);
                child = nodes_[child].child;
            }
        }

        // ノードvの子孫で、最も左にある葉に移動する
        int move_to_leaf(int v) {
            int child = nodes_[v].child;
            while (child != -1) {
                v = child;
                state_.move_forward(nodes_[child].action);
                child = nodes_[child].child;
            }
            return v;
        }

        // ノードvの先祖で、右への分岐があるところまで移動する
        int move_to_ancestor(int v) {
            while (v != root_ && nodes_[v].right == -1) {
                state_.move_backward(nodes_[v].action);
                v = nodes_[v].parent;
            }
            return v;
        }

        // ノードvの右のノードに移動する
        int move_to_right(int v) {
            state_.move_backward(nodes_[v].action);
            v = nodes_[v].right;
            state_.move_forward(nodes_[v].action);
            return v;
        }

        // 不要になった葉を再帰的に削除する
        void remove_leaf(int v) {
            while (true) {
                int left = nodes_[v].left;
                int right = nodes_[v].right;
                if (left == -1) {
                    int parent = nodes_[v].parent;

                    if (parent == -1) {
                        cerr << "ERROR: root is removed" << endl;
                        exit(-1);
                    }
                    nodes_.pop(v);
                    nodes_[parent].child = right;
                    if (right != -1) {
                        nodes_[right].left = -1;
                        return;
                    }
                    v = parent;
                } else {
                    nodes_.pop(v);
                    nodes_[left].right = right;
                    if (right != -1) {
                        nodes_[right].left = left;
                    }
                    return;
                }
            }
        }
};

// ビームサーチを行う関数
vector<Action> beam_search(const Config& config, State state, Node root) {
    Tree tree(state, config.nodes_capacity, root);

    // 探索中のノード集合
    vector<int> curr_nodes;
    curr_nodes.reserve(config.beam_width);
    // 本来は curr_nodes = {state.root_} とすべきだが, 省略しても問題ない

    // 新しいノードの集合
    vector<int> next_nodes;
    next_nodes.reserve(config.beam_width);

    // 新しいノード候補の集合
    Selector selector(config);

    for (int turn = 0; turn < config.max_turn; ++turn) {
        // Euler Tour で selector に候補を追加する
        tree.dfs(selector);

        if (selector.have_finished()) {
            // ターン数最小化型の問題で実行可能解が見つかったとき
            Candidate candidate = selector.get_finished_candidates()[0];
            vector<Action> ret = tree.get_path(candidate.parent);
            ret.push_back(candidate.action);
            return ret;
        }
        // 新しいノードを追加する
        for (const Candidate& candidate : selector.select()) {
            next_nodes.push_back(tree.add_leaf(candidate));
        }
        if (next_nodes.empty()) {
            // 新しいノードがないとき
            cerr << "ERROR: Failed to find any valid solution" << endl;
            return {};
        }
        // 不要なノードを再帰的に削除する
        for (int v : curr_nodes) {
            tree.remove_if_leaf(v);
        }
        // ダブルバッファリングで配列を使い回す
        swap(curr_nodes, next_nodes);
        next_nodes.clear();

        selector.clear();
    }
    // ターン数固定型の問題で全ターンが終了したとき
    int best_leaf = tree.get_best_leaf(curr_nodes);
    return tree.get_path(best_leaf);
}

} // namespace beam_search


それぞれの構造体やクラスについて見ていきます。

Object Pool

ビームサーチと直接の関係がないクラスです。

配列にオブジェクトを保存し、削除したオブジェクトの場所を再利用します。

また、std::vector などと同様に reserve でメモリを確保できるようにしました4


例を用いた説明

最初は長さ4の空の配列とします。

a[0] a[1] a[2] a[3]

3, 1, 4 を順に追加します。追加した場所である 0, 1, 2 を順に返します。

a[0] a[1] a[2] a[3]
3 1 4

a[1] を削除します。

a[0] a[1] a[2] a[3]
3 4

5 を追加します。a[1] と a[3] のどちらに追加してもよいのですが、最後に使用した場所を優先的に使用することにします。今回の例だと a[1] です。追加した場所である 1 を返します。

a[0] a[1] a[2] a[3]
3 5 4

削除した場所のインデックスを保持することでこのような挙動を実装することができます。


差分更新型のビームサーチでは、ノードの追加と削除を頻繁に繰り返します。Object Pool を使うことで次のようなメリットを享受できます。

  • メモリの再利用により空間計算量を削減できる。
  • メモリ上の連続した領域を使用するので、データがキャッシュに乗りやすくなる。
  • 十分な大きさのメモリを最初に確保することで、メモリの再割り当てをなくすことができる。

HashMap

Selector でハッシュ値が重複した候補を除去するところで使います。

open addressing を使用し、インデックスが衝突したときは linear probing を行いました。(ハッシュ値に対する)ハッシュ関数がなく、挿入と一括削除しか行わないので、std::unordered_set よりも単純で高速に動作します。

Hash が整数型でないときは整数型に変換するか、HashMap を std::unordered_set に変更する必要があります。

ちなみに、元々私は std::unordered_set を使っていて、saharan さんのツイート を読んで参考にしたら速度が少し上がりました5

Action

状態遷移に必要な情報をまとめます。ここでいう状態遷移は、親からの移動と親への移動の両方を指しています。

Evaluator

評価値を計算するための構造体です。

基本的な実装は次のようなものです。

struct Evaluator {
    Cost cost;

    Evaluator(Cost cost) : cost(cost) {}

    Cost evaluate() const {
        return cost;
    }
};

Evaluator::evaluate でコストを返します。コストが低いほど採用されやすくなります。

補足: Action や Evaluator で何を保持すべきか

既に述べたように、状態遷移に必要な情報を Action で保持し、評価結果の比較に必要な情報を Evaluator で保持することを想定しています。

一方で、Action と Evaluator に、より多くの情報を持たせることもできます。

状態が変数 x をメンバとして保持し、状態を更新するときに毎度 x を更新するものとします。このとき、x を状態ではなく Action や Evaluator が保持するようにすれば、Euler Tour における x の更新を省略できます。状態の関数内で x を使用したいときには、Action や Evaluator のメンバにアクセスすればよいです。x の後退処理(Euler Tour における子から親への遷移)を実装する必要がなくなるというメリットもあります。

例えば、複数の評価項目があるときに、Evaluator で各評価項目の値を保持するということが考えられます。また、後退処理がなくなるため、浮動小数点数も扱いやすくなります6

一方で、Action や Evaluator の使用メモリが小さいほどよいという側面もあるため、全ての変数を Action や Evaluator に保持すればいいわけではありません。更新が面倒で使用メモリが少ない変数だけを Action や Evaluator に保持するとよいと思います。

例えば、状態内で探索木の深さを管理する場合、深さを更新するときにインクリメントやデクリメントという非常に軽い処理しか行われないため、深さは状態に保持すればよいと思います。

ちなみに、Action と State を空にして Evaluator で全ての情報を保持するようにすると、愚直なビームサーチらしくなります。

Candidate

新しいノードの候補を表現します。

Selector

Candidate の中から新しいノードとして採用するものを選びます。採用された Candidate を元に新しいノードが追加されます。

新しいノードを生成してから不要なものを削除するという実装も考えられますが、木に対する操作は定数倍が重いため、Candidate 構造体を経由する実装になっています。

コードのコメントに書かれているように、ハッシュ値が一致した候補については評価が最もよいものに限定し、その中から評価がよい候補をビーム幅の個数だけ選びます。ハッシュ値の重複を検出するために自作の HashMap を使ったり、ソートをなくすために atcoder::segtree を削除可能な優先度付きキューとして使ったり、高速化を意識して実装しました。

ハッシュ値以外で多様性を確保したい場合、Selector を実装し直す必要があります。例えばスライドパズルなどで「現在のマスの座標が一致するものの中から上位  k 個を選ぶ」という場合には書き直す必要があります。

State

Euler Tour に沿って更新する情報をまとめたクラスです。問題ごとに各メソッドを実装する必要があります。

ビームサーチの最中に State がコピーされることはないため、空間計算量は大きくても構いません。一方で、各メソッドは速いほどよいです。

Tree

探索木を二重連鎖木で表現し、木に対する操作をまとめたクラスです。

重連鎖木のノードは次のノードへのポインタ7を持ちます。

  • 1つ上の兄
  • 1つ下の弟
  • 最も上の子供

状態の更新順序は Euler Tour と一緒です。

Euler Tour の例

一方で、二重連鎖木上では兄弟間を直接移動するようにします。

重連鎖木の遷移

重連鎖木のノードは配列を使用しないため、Euler Tour もノードの追加や削除も簡潔に実装できます。

不要なノードは全て削除します。不要なノードというのは、子ができなかった、あるいは子が全て削除されたノードのことです。

高速なビームサーチが欲しい!!! で紹介されているように、根から一本道の部分は反復しないようにします。

枠で囲った範囲で状態遷移を行う

上図では探索したノードが全て描かれていますが、実際には灰色のノードは作成されず、さらに赤色のノードを子孫として持たない6つの橙色のノードは削除されていることに注意してください。

ビームサーチを実行する関数です。ライブラリの外からこの関数を呼び出します。

使用例

TOYOTA Programming Contest 2023 Summer(AtCoder Heuristic Contest 021) の実装例を紹介します。

大まかな方針

番号が小さいボールから揃えます。

番号が最小のボール、あるいはその左上または右上のボールを、左上か右上に移動させます。

紫色のスワップを遷移の候補とする

評価関数は次のように設定しました。小さいほうがよいです。

 \sum_{i = 0}^{464}{(i \times (ボールiの高さ))} - 600 \times (揃えたボールの数)

ハッシュ値は、揃えているボールの位置と、既に揃えたボールの位置の集合の2つから生成しました。

ライブラリの使い方を紹介することが目的なので、考察などは省略します。

ハッシュ関数

inline int get_pyramid_index(int x, int y) {
    return x * (x - 1) / 2 + y;
}

using Hash = uint32_t;

constexpr Hash hash_mask = ((1U << 23) - 1U) << 9;

inline Hash update_target_position(Hash hash, int x, int y) {
    return (hash & hash_mask) | get_pyramid_index(x, y);
}

inline Hash update_sorted_position(Hash hash, int x, int y) {
    Hash zobrist_hash = get_pyramid_index(x, y);
    zobrist_hash |= 512U; // 10-bit
    zobrist_hash *= zobrist_hash * zobrist_hash; // 30-bit
    return hash ^ (zobrist_hash & hash_mask);
}

下位9ビットで揃えているボールの位置を保持し、残りのビットで既に揃えたボールの位置の集合の Zobrist hashing を保持しました。

Action

struct Action {
    int xyxy;

    Action(int x1, int y1, int x2, int y2) {
        xyxy = x1 | (y1 << 8) | (x2 << 16) | (y2 << 24);
    }

    tuple<int,int,int,int> decode() const {
        return {xyxy & 255, (xyxy >> 8) & 255, (xyxy >> 16) & 255, xyxy >> 24};
    }
};

AHC021の場合、スワップ位置が分かれば盤面の更新が行えるため、スワップする2つの位置を保持しました。

4つの整数を1つの intエンコードし、メモリ使用量を減らしています8

Evaluator

using Cost = int;

constexpr int target_coefficient = 600;

struct Evaluator {
    int target_ball;
    int potential;

    Evaluator(int target_ball, int potential) :
        target_ball(target_ball),
        potential(potential) {}

    Cost evaluate() const {
        return potential - target_coefficient * target_ball;
    }
};

2つの評価項目を保持しました。

State

更新部分だけ説明します。他の部分を読みたい場合は提出コードをご覧ください。

class State {
    public:
        void move_forward(Action action) {
            auto [x1, y1, x2, y2] = action.decode();
            swap_balls(x1, y1, x2, y2);
        }

        void move_backward(Action action) {
            auto [x1, y1, x2, y2] = action.decode();
            swap_balls(x1, y1, x2, y2);
        }

    private:
        vector<vector<int>> b_;
        array<pair<int,int>,m> positions_;

        void swap_balls(int x1, int y1, int x2, int y2) {
            int b1 = b_[x1][y1];
            int b2 = b_[x2][y2];
            b_[x1][y1] = b2;
            b_[x2][y2] = b1;
            positions_[b2] = {x1, y1};
            positions_[b1] = {x2, y2};
        }
};

b_ は盤面を表しています。positions_[i] はボール i の位置を表します。 positions_[b_[x][y]] = {x, y} が成り立ちます。positions_ は、あるボールを揃えた後に次のボールの位置を得るときなどに使われます。

Evaluator が揃えたボールの数を保持するため、現在揃えているボールの番号を保持・更新する必要がないことに注意してください。

今回は2つのボールを入れ替えるだけなので move_forwardmove_backward が等しくなっていますが、一般的には異なります。

提出コード

ビーム幅を3500に設定しました。

提出言語は C++ 20 (Clang 16.0.6) です。私のビームサーチライブラリ(の主に Selector)は GCC よりも Clang のほうが高速に動作するようでした。

atcoder.jp

他の使用例

ゲーム実況者Xの挑戦 の提出コードです。解説などは省略します。

atcoder.jp

Euler Tour の辺を保持する実装 (2024/02/07 追記)

重連鎖木や Object Pool を使うのではなく、探索木の有向辺を Euler Tour の順序で保持するほうが速いという話があり、実装してみました。

状態を差分計算するときは Euler Tour の配列に前から順にアクセスします。

探索木を更新するときは、辺の追加や削除をしながら、別の配列に Euler Tour をコピーしています。


ビームサーチライブラリ

#include <bits/stdc++.h>

using namespace std;

namespace beam_search {

// ビームサーチの設定
struct Config {
    int max_turn;
    size_t beam_width;
    size_t tour_capacity;
    uint32_t hash_map_capacity;
};

// 連想配列
// Keyにハッシュ関数を適用しない
// open addressing with linear probing
// unordered_mapよりも速い
// nは格納する要素数よりも16倍ほど大きくする
template <class Key, class T>
struct HashMap {
    public:
        explicit HashMap(uint32_t n) {
            if (n % 2 == 0) {
                ++n;
            }
            n_ = n;
            valid_.resize(n_, false);
            data_.resize(n_);
        }

        // 戻り値
        // - 存在するならtrue、存在しないならfalse
        // - index
        pair<bool,int> get_index(Key key) const {
            Key i = key % n_;
            while (valid_[i]) {
                if (data_[i].first == key) {
                    return {true, i};
                }
                if (++i == n_) {
                    i = 0;
                }
            }
            return {false, i};
        }

        // 指定したindexにkeyとvalueを格納する
        void set(int i, Key key, T value) {
            valid_[i] = true;
            data_[i] = {key, value};
        }

        // 指定したindexのvalueを返す
        T get(int i) const {
            assert(valid_[i]);
            return data_[i].second;
        }

        void clear() {
            fill(valid_.begin(), valid_.end(), false);
        }

    private:
        uint32_t n_;
        vector<bool> valid_;
        vector<pair<Key,T>> data_;
};

using Hash = uint32_t; // TODO

// 状態遷移を行うために必要な情報
// メモリ使用量をできるだけ小さくしてください
struct Action {
    // TODO

    Action() {
        // TODO
    }

    bool operator==(const Action& other) const {
        // TODO
    }
};

using Cost = int;

// 状態のコストを評価するための構造体
// メモリ使用量をできるだけ小さくしてください
struct Evaluator {
    // TODO

    Evaluator() {
        // TODO
    }

    // 低いほどよい
    Cost evaluate() const {
        // TODO
    }
};

// 展開するノードの候補を表す構造体
struct Candidate {
    Action action;
    Evaluator evaluator;
    Hash hash;
    int parent;

    Candidate(Action action, Evaluator evaluator, Hash hash, int parent) :
        action(action),
        evaluator(evaluator),
        hash(hash),
        parent(parent) {}
};

// ノードの候補から実際に追加するものを選ぶクラス
// ビーム幅の個数だけ、評価がよいものを選ぶ
// ハッシュ値が一致したものについては、評価がよいほうのみを残す
class Selector {
    public:
        explicit Selector(const Config& config) :
            hash_to_index_(config.hash_map_capacity)
        {
            beam_width = config.beam_width;
            candidates_.reserve(beam_width);
            full_ = false;

            costs_.resize(beam_width);
            for (size_t i = 0; i < beam_width; ++i) {
                costs_[i] = {0, i};
            }
        }

        // 候補を追加する
        // ターン数最小化型の問題で、candidateによって実行可能解が得られる場合にのみ finished = true とする
        // ビーム幅分の候補をCandidateを追加したときにsegment treeを構築する
        void push(const Candidate& candidate, bool finished) {
            if (finished) {
                finished_candidates_.emplace_back(candidate);
                return;
            }
            Cost cost = candidate.evaluator.evaluate();
            if (full_ && cost >= st_.all_prod().first) {
                // 保持しているどの候補よりもコストが小さくないとき
                return;
            }
            auto [valid, i] = hash_to_index_.get_index(candidate.hash);

            if (valid) {
                int j = hash_to_index_.get(i);
                if (candidate.hash == candidates_[j].hash) {
                    // ハッシュ値が等しいものが存在しているとき
                    if (full_) {
                        // segment treeが構築されている場合
                        if (cost < st_.get(j).first) {
                            candidates_[j] = candidate;
                            st_.set(j, {cost, j});
                        }
                    } else {
                        // segment treeが構築されていない場合
                        if (cost < costs_[j].first) {
                            candidates_[j] = candidate;
                            costs_[j].first = cost;
                        }
                    }
                    return;
                }
            }
            if (full_) {
                // segment treeが構築されている場合
                int j = st_.all_prod().second;
                hash_to_index_.set(i, candidate.hash, j);
                candidates_[j] = candidate;
                st_.set(j, {cost, j});
            } else {
                // segment treeが構築されていない場合
                int j = candidates_.size();
                hash_to_index_.set(i, candidate.hash, j);
                candidates_.emplace_back(candidate);
                costs_[j].first = cost;

                if (candidates_.size() == beam_width) {
                    // 保持している候補がビーム幅分になったときにsegment treeを構築する
                    full_ = true;
                    st_ = MaxSegtree(costs_);
                }
            }
        }

        // 選んだ候補を返す
        const vector<Candidate>& select() const {
            return candidates_;
        }

        // 実行可能解が見つかったか
        bool have_finished() const {
            return !finished_candidates_.empty();
        }

        // 実行可能解に到達するCandidateを返す
        vector<Candidate> get_finished_candidates() const {
            return finished_candidates_;
        }

        // 最もよいCandidateを返す
        Candidate calculate_best_candidate() const {
            if (full_) {
                size_t best = 0;
                for (size_t i = 0; i < beam_width; ++i) {
                    if (st_.get(i).first < st_.get(best).first) {
                        best = i;
                    }
                }
                return candidates_[best];
            } else {
                size_t best = 0;
                for (size_t i = 0; i < candidates_.size(); ++i) {
                    if (costs_[i].first < costs_[best].first) {
                        best = i;
                    }
                }
                return candidates_[best];
            }
        }

        void clear() {
            candidates_.clear();
            hash_to_index_.clear();
            full_ = false;
        }

    private:
        // 削除可能な優先度付きキュー
        using MaxSegtree = atcoder::segtree<
            pair<Cost,int>,
            [](pair<Cost,int> a, pair<Cost,int> b){
                if (a.first >= b.first) {
                    return a;
                } else {
                    return b;
                }
            },
            []() { return make_pair(numeric_limits<Cost>::min(), -1); }
        >;

        size_t beam_width;
        vector<Candidate> candidates_;
        HashMap<Hash,int> hash_to_index_;
        bool full_;
        vector<pair<Cost,int>> costs_;
        MaxSegtree st_;
        vector<Candidate> finished_candidates_;
};

// 深さ優先探索に沿って更新する情報をまとめたクラス
class State {
    public:
        explicit State() {
            // TODO
        }

        // EvaluatorとHashの初期値を返す
        pair<Evaluator,Hash> make_initial_node() {
            // TODO
        }

        // 次の状態候補を全てselectorに追加する
        // 引数
        //   evaluator : 今の評価器
        //   hash      : 今のハッシュ値
        //   parent    : 今のノードID(次のノードにとって親となる)
        void expand(const Evaluator& evaluator, Hash hash, int parent, Selector& selector) {
            // TODO
        }

        // actionを実行して次の状態に遷移する
        void move_forward(Action action) {
            // TODO
        }

        // actionを実行する前の状態に遷移する
        // 今の状態は、親からactionを実行して遷移した状態である
        void move_backward(Action action) {
            // TODO
        }

    private:
        // TODO
};

// Euler Tourを管理するためのクラス
class Tree {
    public:
        explicit Tree(const State& state, const Config& config) :
            state_(state)
        {
            curr_tour_.reserve(config.tour_capacity);
            next_tour_.reserve(config.tour_capacity);
            leaves_.reserve(config.beam_width);
            buckets_.assign(config.beam_width, {});
        }

        // 状態を更新しながら深さ優先探索を行い、次のノードの候補を全てselectorに追加する
        void dfs(Selector& selector) {
            if (curr_tour_.empty()) {
                // 最初のターン
                auto [evaluator, hash] = state_.make_initial_node();
                state_.expand(evaluator, hash, 0, selector);
                return;
            }

            for (auto [leaf_index, action] : curr_tour_) {
                if (leaf_index >= 0) {
                    // 葉
                    state_.move_forward(action);
                    auto& [evaluator, hash] = leaves_[leaf_index];
                    state_.expand(evaluator, hash, leaf_index, selector);
                    state_.move_backward(action);
                } else if (leaf_index == -1) {
                    // 前進辺
                    state_.move_forward(action);
                } else {
                    // 後退辺
                    state_.move_backward(action);
                }
            }
        }

        // 木を更新する
        void update(const vector<Candidate>& candidates) {
            leaves_.clear();

            if (curr_tour_.empty()) {
                // 最初のターン
                for (const Candidate& candidate : candidates) {
                    curr_tour_.push_back({(int)leaves_.size(), candidate.action});
                    leaves_.push_back({candidate.evaluator, candidate.hash});
                }
                return;
            }

            for (const Candidate& candidate : candidates) {
                buckets_[candidate.parent].push_back({candidate.action, candidate.evaluator, candidate.hash});
            }

            auto it = curr_tour_.begin();

            // 一本道を反復しないようにする
            while (it->first == -1 && it->second == curr_tour_.back().second) {
                Action action = (it++)->second;
                state_.move_forward(action);
                direct_road_.push_back(action);
                curr_tour_.pop_back();
            }

            // 葉の追加や不要な辺の削除をする
            while (it != curr_tour_.end()) {
                auto [leaf_index, action] = *(it++);
                if (leaf_index >= 0) {
                    // 葉
                    if (buckets_[leaf_index].empty()) {
                        continue;
                    }
                    next_tour_.push_back({-1, action});
                    for (auto [new_action, evaluator, hash] : buckets_[leaf_index]) {
                        int new_leaf_index = leaves_.size();
                        next_tour_.push_back({new_leaf_index, new_action});
                        leaves_.push_back({evaluator, hash});
                    }
                    buckets_[leaf_index].clear();
                    next_tour_.push_back({-2, action});
                } else if (leaf_index == -1) {
                    // 前進辺
                    next_tour_.push_back({-1, action});
                } else {
                    // 後退辺
                    auto [old_leaf_index, old_action] = next_tour_.back();
                    if (old_leaf_index == -1) {
                        next_tour_.pop_back();
                    } else {
                        next_tour_.push_back({-2, action});
                    }
                }
            }
            swap(curr_tour_, next_tour_);
            next_tour_.clear();
        }

        // 根からのパスを取得する
        vector<Action> calculate_path(int parent, int turn) const {
            // cerr << curr_tour_.size() << endl;

            vector<Action> ret = direct_road_;
            ret.reserve(turn);
            for (auto [leaf_index, action] : curr_tour_) {
                if (leaf_index >= 0) {
                    if (leaf_index == parent) {
                        ret.push_back(action);
                        return ret;
                    }
                } else if (leaf_index == -1) {
                    ret.push_back(action);
                } else {
                    ret.pop_back();
                }
            }

            unreachable();
        }

    private:
        State state_;
        vector<pair<int,Action>> curr_tour_;
        vector<pair<int,Action>> next_tour_;
        vector<pair<Evaluator,Hash>> leaves_;
        vector<vector<tuple<Action,Evaluator,Hash>>> buckets_;
        vector<Action> direct_road_;
};

// ビームサーチを行う関数
vector<Action> beam_search(const Config& config, const State& state) {
    Tree tree(state, config);

    // 新しいノード候補の集合
    Selector selector(config);

    for (int turn = 0; turn < config.max_turn; ++turn) {
        // Euler Tourでselectorに候補を追加する
        tree.dfs(selector);

        if (selector.have_finished()) {
            // ターン数最小化型の問題で実行可能解が見つかったとき
            Candidate candidate = selector.get_finished_candidates()[0];
            vector<Action> ret = tree.calculate_path(candidate.parent, turn + 1);
            ret.push_back(candidate.action);
            return ret;
        }

        assert(!selector.select().empty());

        if (turn == config.max_turn - 1) {
            // ターン数固定型の問題で全ターンが終了したとき
            Candidate best_candidate = selector.calculate_best_candidate();
            vector<Action> ret = tree.calculate_path(best_candidate.parent, turn + 1);
            ret.push_back(best_candidate.action);
            return ret;
        }

        // 木を更新する
        tree.update(selector.select());

        selector.clear();
    }

    unreachable();
}

} // namespace beam_search


元の二重連鎖木を使った実装よりも高速に動作するようでした。

atcoder.jp

atcoder.jp

最後に

この記事では私の実装のみを紹介しました。実装した人によって異なる部分があるので調べてみると面白いかもしれません。

最後まで読んでくださりありがとうございました。


  1. 若干意味合いが異なるかもしれません。私は同一視しています。
  2. 履歴は永続配列を使用することで高速化できます。
  3. 状態遷移の計算量が状態をコピーする計算量と同程度の場合には愚直なビームサーチでよいと思います。過去のAHCを見る限りだと、状態遷移が  O(1) の場合が多く、愚直なビームサーチが強い問題は少ないです。
  4. 内部的には std::vector::reserve を呼び出しているだけです。
  5. open addressing と linear probing について参考にさせていただきました。私の実装では、遅延削除は行っていません。
  6. 浮動小数点演算を含む関数に対して、正確な逆関数を定義することは難しいことが多いです。例えば  y = 1.1 \times x として  y から  x を復元するときに、 y / 1.1 を行うと思いますが、数値誤差を考慮すると  y / 1.1 x と等しいとは限りません。状態の更新に浮動小数点演算が絡むと、数値誤差によって正しく復元されないリスクがあります。
  7. 正確には ObjectPool におけるインデックスです。
  8. 16ビットにエンコードすることも可能ですが、個人的には限界まで減らすモチベーションがなかったです。

AHC029 参加記

公式ビジュアライザ (seed=2)

RECRUIT 日本橋ハーフマラソン 2024冬(AtCoder Heuristic Contest 029) お疲れ様でした。

システムテストの得点率が 73.4 % で優勝しました!

順位表

問題概要

所持金の最大化を行うカードゲームです。

詳細は公式の問題文を参考にしてください。

atcoder.jp

方針

貪欲法をモンテカルロ法によって改善しました。

コンテスト上位者の中で私の貪欲法の性能は悪かったため、貪欲法についてはあまり参考にならないかもしれません。

貪欲法の準備

貪欲法を説明する前に、貪欲法の説明に必要なものについて述べます。

各ターンの流れ

各ターンの流れについて、問題文では

  1. 手札からカードを1枚選んで使う。
  2. 必要に応じてプロジェクが補充される。
  3. 補充するカードの候補が提示される。
  4. 補充するカードを1枚選んで補充する。

の順序で書かれています。

しかし、補充するカードを選んだ直後に手札のカードを使うため、これら2つの処理をまとめて考えることができます。

そこで、各ターンの流れを

  1. 補充するカードの候補が提示される。
  2. 補充するカードを1枚選んで補充し、手札からカードを1枚選んで使う。
  3. 必要に応じてプロジェクトが補充される。

のように捉えることにします。

1ターン目は、手札の1枚目のカードが欠けていて、補充するカードの候補として手札の1枚目のカードだけが提示されていることにすると実装しやすいです。また、問題文における最終ターンでは、補充するカードとして1枚目を選ぶのが最適なので、最後のカードを選ぶ手順は省略して考えることができます。

行動の定義

(補充するカード、使用するカードと対象のプロジェクト)という組を行動と定義します。問題文中の記号で表すと  (r, c, m) です。

貪欲法

行動価値関数を定義し、行動価値が最大になる行動を選択します。

直後に取りうる行動を全探索するため、ループ内部の時間計算量を  O(1) とすると、全体の計算量は  O(KNM) になります。

行動価値関数は、後述する状態価値関数とは直接的な関係がないことに注意してください。(強化学習における行動価値関数とは定義が異なると思います。)

行動価値関数は、以下の項目の変化を評価します。

  • 所持金  money
  • プロジェクト
  • 手札のカード
  • 使用した投資カードの枚数  L

単位は所持金に合わせることにします。すなわち、所持金以外の項目の変化をお金に換算します。所持金ベースで考えることで、定量的な考察がしやすくなると思います。

所持金の変化

カードのコスト  p と、カードを使うことによって完了するプロジェクトの価値  v を用いて計算することができます。

プロジェクトの残務量を減らす価値

プロジェクトが完了する場合  (h \leq w)、所持金の変化のみを考えることにします。すなわち、新しいプロジェクトの評価は0とします。

プロジェクトが完了しない場合  (h \gt w) v w / h をベースにして評価しました。減らした残務量の割合に応じてプロジェクトの価値を得られるイメージです。

実際にはお金を稼げないことや、累積和が  v よりも大きくなることを踏まえると、 v w / h を下方修正したほうがよさそうです。

また、カード生成方法を見ると、所持金と減らせる残務量が大体同じであることがわかります。そのため、所持金よりも大きい残務量のプロジェクトを完了するには時間がかかることが予想されます。

そこで、 \max{(0, h - w - money)} が大きいほど、 v w / h を下方修正するようにしました。


具体的な修正方法とその高速化

 v w / h に以下の値をかけました。

 1.1 / (1 + 0.5 \times (1 + \lfloor \max{(0, h - w - money)} / 2^{L} \rfloor )^{0.4})

0.4乗の計算が重そうに見えますが、0.4乗する前の値が257以下の自然数であることから、事前計算によって配列アクセスに変更できます。0.4乗以降の計算も事前計算に置き換えることにより、浮動小数点演算の回数を大幅に減らすことができます。


プロジェクトを取りやめる価値

 h - v をベースとしました。

ただし、 h が小さいほど、取りやめる価値を下げました。以下の2つの理由があります。

  •  h が小さければすぐに完了できるから。
  • 途中まで進めたプロジェクトを業務転換カードで換えないようにするため。

カードの変化

カードの価値を計算すればよいです。

ゲーム終了時に強いカードを残すのは無駄なので、最後20ターンほどからカードの価値を下げるようにしています。

通常労働カードの価値

 0.9 \times w としました。プロジェクトの残務量を0未満にする(overkill と呼ぶことにします)ことがあるため、 w よりも少しだけ小さくしています。

全力労働カードの価値

 0.6 \times M w としました。overkill が通常労働よりも多発することや、効率の悪いプロジェクトも進めるため、 M w よりも小さくしています。

キャンセルカードと業務転換カードの価値

 2^{L} としました。

増資について

900ターンぐらいまで、増資カードは買えるときに必ず買うようにし、買ったらすぐに使用しました。評価値が無限大と考えることもできます。

実は増資カードが高いときには買わないべきで、増資ストックも検討するべきなのですが、コンテスト中にうまく実装できませんでした。今後暇なときに実装したいです。


増資ストックとは

増資カードを購入し、しばらく使わずに手札に残すことを増資ストックと呼ぶことにします。

増資すると多くの変数が2倍になりますが、所持金は2倍になりません。相対的に所持金が半分になると考えることができます。

所持金が少ないと効率が悪くなりやすいため、所持金が少ない期間は短いほうがよいです。

以上より、所持金が少ないときに増資カードを一気に使いたいというモチベーションが生まれます。

一方で、所持金を一時的でも極端に減らしていいのか、増資ストック中に使える手札のカード枚数を減らしていいのかなど、増資ストックのデメリットに関しては自分の中で結論が出てないです。

(2024/01/08 追記: 増資ストックではなく増資スタックでした。同じアイテムを複数まとめることをスタックというらしいです。)


貪欲法の高速化

モンテカルロ法におけるボトルネックは貪欲法になります。そのため、貪欲法は速ければ速いほどよいです。

処理をまとめる

基本的に可能な行動全てについて行動価値を計算しますが、補充したカードの価値など、いくつかの行動で共通する要素があるので、それらはなるべく1回しか計算しないようにしました。

行動価値関数として実装しないため、実装が汚くなりました。

貪欲法における枝刈り

行動を全探索すると書きましたが、明らかに無駄な探索は枝刈りします。

枝刈りするのは以下のケースです。

  • 下位互換のカードは補充しない。すなわち、補充するカード候補の中に、同じ種類でコスト  p が小さくて効能がそれ以上のものがあるようなカードを補充することを禁止する。
  • 手札に存在する効能が同じカードを同一視する。

枝刈りは貪欲法を高速化するだけでなく、モンテカルロ法でプレイアウトを行う行動の候補を選ぶときにも効果的に働きます。

状態評価関数

モンテカルロ法で状態を評価するときに使います。

状態評価値の比較は同じターン同士でしか行われないため、同じターンで評価基準が揃っていればよいものとします。

ゲームが終了している場合、所持金を評価しました。以降、ゲームが終了していない場合について説明します。

基本的には行動価値関数と同じ考え方で状態を評価します。

所持金の評価

所持金が  200 \times 2^{L} より大きいときは増資カードを買える可能性があるので所持金に応じて高く評価しました。

カードの評価

行動価値関数とほぼ一緒です。

プロジェクトの評価

 v - h をベースとしました。

ただし、キャンセルカードや業務転換カードを保持しているときは  v \lt h の場合のペナルティを緩和しました。

増資カード使用回数の評価

 425 \times 2^{L} としました。

その他

最後に評価値を 0.95 乗して極端な値の影響を下げました。スコアはあまり変わらなかったと思います。

モンテカルロ法

可能な行動について、貪欲法によるプレイアウトで評価値の期待値を求め、期待値が最大の行動を選択します。

処理の流れは以下の通りです。

  1. 行動価値が大きい7つの行動を候補として取得する。
  2. およそ2msec経過するまで以下の処理を繰り返す。
    1. プレイアウトに必要な情報を生成する。
    2. 各候補について9ターンの貪欲プレイアウトを行う。
    3. 一定時間が経過した場合、期待値が最小の候補を削除する。
  3. 期待値が最大の候補を選択する。

重要そうな部分について説明します。

シミュレーション回数

より正確な期待値を求めるにはプレイアウトの回数を増やすことが大切です。

そこで、プレイアウト回数を増やす工夫をいくつか行いました。

プレイアウトを行う行動の選択

可能な行動が  O(KNM) 通りあり、理想的にはそれら全てについて十分な回数のプレイアウトを行いたいです。

しかし、実際には実行時間が非常に限られているため、可能な行動全てに対してプレイアウトを行うとプレイアウト回数が減ってしまいます。

そこで、行動価値が大きい行動についてのみプレイアウトを行うようにしました。行動価値関数は完璧ではありませんが、上位何個か選べば最良のものが含まれているだろうという感覚です。

また、プレイアウトを行う行動候補を選ぶときに、補充するカードと使用するカードの組が等しいものは1個までとしました。すなわち、候補の中にカードを適用するプロジェクトだけが異なるような行動がないようにしました。おそらく候補の多様性が重要なのだと思います。

データ生成

基本的には問題文と同様の方法でプロジェクトとカードについてのデータを生成します。

プロジェクトの生成には隠しパラメータが使用されていないため、一度生成したプロジェクトを別のターンでも使いまわしました。

一方で、カードの生成には隠しパラメータが使用されているため、今までに出現したカードの統計データをとり、毎ターンその分布にしたがってカードを生成しました。

統計データは各カードの種類について出現した回数を数えるだけですが、出現した回数の初期値は 21, 11, 11, 6, 4 としました。序盤はデータが少ないので無情報の期待値を参考にできるということです。スコアへの寄与は小さいと思います。

プレイアウト

狭義的には最終状態までシミュレーションすることをプレイアウトと言いますが、途中までシミュレーションすることもプレイアウトと呼ぶことにします。

最終状態までシミュレーションした場合、序盤や中盤では非常に長い時間がかかるうえ、ランダム性が大きくなって評価したい行動のよさが計りにくくなります。そこで、シミュレーションを9ターン後まで行い、9ターン後の状態価値を最大化することを目標としました。ただし、9ターン後までにゲームが終了する場合はゲームを終了させ、そのときの評価は所持金になります。

候補の削除

7つの行動候補を段階的に減らし、最終的に2つの候補についてプレイアウトを行うようにしました。一定時間が経過するごとに最も期待値が低い候補を1つ削除します。

最良でなさそうな候補を途中で削除することにより、よさそうな候補のシミュレーション回数を増やすことができます。

感想

全体

モンテカルロ法について、AHCの過去問で何度も練習していて、今回結果を出せて非常に嬉しく思います。

問題文を読み終えた時点で、数ターンのプレイアウトを行うモンテカルロ法が強いことを確信し、自分の得意ジャンルであることから優勝するチャンスだと思いました。

大まかな方針はすぐに決まりましたが、「数ターンのプレイアウト」の「数ターン」は20ターンぐらいだと予想していましたし、モンテカルロ法によって、スコアはせいぜい2倍ぐらいにしかならないだろうと思っていたので、その辺りの感覚はずれていました。(スコアは平均して32倍程度、順位表のスコアは12倍程度でした。)

優勝できたのはよかったですが、感想戦において貪欲法で大きく負けていたことがわかり、貪欲法については反省点が多そうです。貪欲法(と状態評価関数)が基本で、モンテカルロ法を行う場合でも、最終的なスコアは貪欲法に大きく依存すると思っていたのでかなり意外な結果でした。

相対評価システム

自分がコンテストを荒らしている感じがして楽しかったです。

今回の相対評価が気に入らない人がいるみたいなので個人的な意見を少し書きます。

順位スコアはテストケースにおける圧倒的1位の人にとって嬉しくないので基本的にはよくないと思っています。元のスコアが指数的な分布なら、対数をとって絶対スコアか相対スコアだとバランスがいいかもしれないです。指数的な分布でないとき(一様分布など)に対数をとると、僅差になって非ACのペナルティが大きくなるのでよくなさそうです。

バグ

自分の実装において、直すとスコアが悪化するバグが2つあったので紹介します。

状態評価関数におけるカードの評価

私の状態評価関数は、カードを使ってプロジェクトが補充された直後の状態に適用するようになっています。カードが1枚欠けた状態です。

バグで欠けた1枚のカード、すなわち最後に使用したカードも評価に入れていることに気づき、バグを直したところスコアが少し下がりました。

原因に心当たりはありません。単なるスコアのぶれかもしれないです。

カードの種類の統計の取り方

カードの種類の統計をとるとき、提示されたカードの1枚目は必ず通常労働カードで統計データに含めてはいけないのですが、バグで統計に含まれていました。バグを直したところ、スコアが少し下がりました。

 K が小さいときに大きく影響しそうですが、通常労働カードは基本的なカードなので、少し出やすくしたほうがいいのでしょうか。納得はしていませんが、致命的なバグではないような気がしています。

最後に

コンテストを主催してくださったリクルート様、学生賞金があって大変感謝しております。長期AHCとしては珍しいタイプの問題で、とても面白かったです。

AtCoderの皆様、参加者の皆様、最後まで読んでくださった読者の皆様、ありがとうございました。

ALGO ARTIS プログラミングコンテスト2023 冬(AtCoder Heuristic Contest 028) も参加する予定なので対戦よろしくお願いします。

AHC027 参加記

Seed = 0, Score = 1317173

HACK TO THE FUTURE 2024 (AtCoder Heuristic Contest 027) お疲れ様でした。

システムテストの得点率が 98.44 % で準優勝しました!

順位表

前提知識

  • 木上のビームサーチ
  • バックトラック法
  • low-link

問題概要

 N \times N マスの二次元グリッドがあり、各マス  (i, j) の汚れは毎ターン  d_{i, j} だけ増加します。ロボット掃除機が訪れたマスは汚れが 0 になります。ロボットの掃除ルートで、平均汚れができるだけ小さいものを求めてください。

詳細は公式の問題文を参考にしてください。

atcoder.jp

制約

 20 \leq N \leq 40

方針

ビームサーチで初期解を生成し、焼きなまし法を適用しました。

私が調べた限りでは、bowwowforeach さんや siman さんと大まかな方針は一緒で、他の上位の人は方針が違いそうでした。

ビームサーチ

深さ

深さ1 t_{max} := 2.7 \times N^{2} に設定しました。

ただし、深さが偶数になるように必要に応じてインクリメントしました。深さを偶数にすることで始点と終点が同じものを得ることができます。

多様性の確保

(現在のマス, 直前にいたマス) をハッシュ値の代わりに使用しました。ビーム幅は  2 N^{2} 弱ということになります。

次のノードを選ぶときにソートする必要がなくなるため、ビームサーチが高速になります。

評価関数

4つの項目を評価しました。

影響が大きそうなものから順に説明します。

評価項目1: 汚れの累積和

後述する汚れの総和を毎ターン加算したものです。

評価項目2: 汚れの総和

現在の各マスの汚れの総和です。初期状態における汚れは  4000 \times d_{i, j} で初期化しました。

汚れの総和は、直前の汚れの総和と、そのターンに回収した汚れの量から計算することができます。回収した汚れの量を求めるには、各マスについて最後に訪れたターンをメモしておけばよいです。

 d の代わりに  \max(d, 5)^{0.85} を使うとスコアがよくなりました。 d をそのまま使うと  d が大きいマスばかり訪問するのに対し、修正したものを使うと比較的バランスよく訪問していました。

評価項目3: 移動平均

指数移動平均を計算し、絶対値が大きいほどよい評価をしました。

ここでいう指数移動平均とは、 t ターン目の座標の変化量を  Y_{t} として

 \sum_{t=0}^{T}{\alpha (1 - \alpha)^{T - t} Y_t}

を表すものとします。平滑化係数は  \alpha = 0.05 としました。

更新は

 X_{t} = \alpha Y_{t} + (1 - \alpha) X_{t - 1}

のように行います。

 d が大きい狭い領域に居続けるのがよくないので、移動平均の絶対値が大きいものを優先するとよいと考えられます。

例えば Seed = 25 は、 d が大きい領域が非常に限られており、評価関数によっては中心やや下の領域から出なくなります。

Seed = 25, d の分布

評価項目4: 周期解としての生スコア

汚れの累積和では、掃除ルートが繰り返されることが考慮されていません。そこで、 t_{max} ターン目まで何もしなかった場合の生スコアを評価しました。

各マスを訪問するたびに、訪問時刻を挿入するような感じで差分計算を行うことができます。

未訪問のマスは  2 t_{max} ターンの周期で訪問されるものとして計算しました。

差分更新

木上のビームサーチでは、評価の更新に必要な情報を深さ優先探索に沿って更新します。私の実装だと、評価値はノードに保存されるので深さ優先探索では更新する必要がありません。

評価の更新に必要な情報で動的なものは、現在のターンと各マスを最後に訪れた時刻の2つです。これらは  O(1) で更新できます。

簡易的な C++ のコードを載せておきます。座標が一次元に平坦化されていることに注意してください。

struct State {
    // 経過したターン数
    int turn;

    // timestamps[i] = マスiを訪問した時刻の配列
    vector<vector<int>> timestamps;

    // 深さ優先探索における子孫への移動
    void move_forward(int position) {
        timestamps[position].push_back(++turn)
    }

    // 深さ優先探索における親への移動
    void move_backward(int position) {
        assert(timestamps[position].back() == turn)
        timestamps[position].pop_back();
        --turn;
    }
}

前回訪問した時刻をノードに保存することで、一次元配列のみで実装することも可能です。2

未訪問頂点の対応

未訪問頂点が減るように評価関数を工夫しているつもりですが、それでも未訪問の頂点が残る場合があります。

なるべく低コストで未訪問頂点を挿入したいので、2マス拡張(後述)を基本として、2マス拡張ができないときは1マス拡張(後述)を行いました。

元々

元の経路

だった経路を

1マス拡張の例

のように拡張するのが1マス拡張で

2マス拡張の例

のように拡張するのが2マス拡張です。

ビームサーチの結果

コンテスト後にビームサーチのみのコードを提出したところ、44Gほどのスコアが出ていました。

焼きなまし法

3種類の近傍で焼きなまします。

近傍1: L字型の変更

a, b, c を a, d, c のように変更します。L字型の部分を反対側に曲げる感じです。

L字型の変更の例

経路長が変わらないので高速に計算できます。

近傍2: 十字型の解消

掃除ルートが同じマスでクロスしていた場合に、区間の訪問順序を反転させてクロスを解消します。2-optをイメージするとよいかもしれません。

クロスを解消した後にL字型の変更が行われると1マス分得をすることができます。

例えば

クロスが解消される前

区間を反転させて、クロスしていたところにL字型の変更を施すと

クロスが解消されてL字型の変更が行われた後

になります。

反転区間の長さが約150以上の場合にはクロスの解消を行いませんでした。区間が長すぎると訪問時間の変更によるデメリットのほうが大きくなります。

近傍3: バックトラック法による破壊再構築

掃除ルートの連続する 4 ~ 8 マスを破壊し、バックトラック法で生成した経路で置き換えます。

バックトラック法の工夫1: 同じマスを通らない

よくあるバックトラック法と同様に、同じマスを複数回通らないようにしました。

ただし、関節点については2回まで同じマスを通ってもよいものとしました。関節点も1回しか通れないようにした場合、関節点の削除はできても挿入はできないので、操作が不可逆的になります。

関節点への侵入

関節点については、low-link で事前に列挙しました。

バックトラック法の工夫2: 経路長の制限

経路長を無制限にすると探索空間が膨大すぎて計算が終わりません。そこで、経路長は「元の経路長 + 2」までとしました。

バックトラック法の工夫3: 経路長による枝刈り

経路長の最大値を設定したため、経路長による枝刈りができます。

具体的には、現在のパスの長さと現在位置から終点までの距離を足したものが、経路長の最大値よりも大きいならば枝刈りできます。

全頂点から幅優先探索をするのは計算量的に少し重そうだったので、距離はマンハッタン距離で代用しました。

バックトラック法の工夫4: よさそうな頂点から探索する

次の頂点を選ぶときに、訪問したときに得られるスコアの改善量が大きい頂点を優先しました。

具体的にはスコアの改善量を  \delta として、 \delta \times \rm{randint}(10, 50) が大きいものから順に探索しました。

工夫できなかったこと

経路長が変わる場合、変更部分以降の時刻が全てずれ、更新の計算量が大きくなりました。

ただ更新するだけだともったいなかったので、経路全体を逆順にしたり、掃除ルートのスタート地点を変更したりしました。ときどきこの処理が行われることを前提とし、各近傍では掃除ルートの終わりと始まりをまたぐような遷移は考えませんでした。

各近傍の選択比率

近傍1, 2, 3 の選択比率は 6 : 3 : 1 としました。選択比率と実行時間は異なることに注意してください。

途中で掃除ルートを2倍にする

焼きなましの実行時間の8割が経過したときに、同じ掃除ルートを2回繰り返したものに変更しました。

掃除ルートが長くなって焼きなましにくくなるデメリットがある一方で、全頂点を2回訪れていると制約が緩くなって焼きなましやすくなったり、汚れやすさが極端に低い頂点の掃除頻度を半減できたりするメリットがありそうです。

試行回数

Seed = 0 ( N = 20) で約 2,500,000 回のループを回せていました。

日記

暇な人向けです。

12/1

  • 問題を読む。
  • 適当に初期解を作り、破壊再構築の焼きなましをしたい。
    • 経路長が変わると差分計算が面倒。
  • 遺伝的アルゴリズムとか使えないかな。無理か。
  • 移動が自由なら Introduction to Heuristic Contest に似ているな。新ジャッジコン で使った貪欲とか活かせるといいな。
  • 掃除ルートの周期性を考えなくても、経路が十分長ければ接続部分は無視できそう。

12/2

  • BFSに沿ったDPをベースとした貪欲とかがいいのかな。
  • 掃除できた汚れの量を評価関数にしてビームサーチを打つほうがいいか。
    • 今のマスをハッシュ値として多様性を確保できる。
    • 差分更新ビームサーチライブラリの出番が来た!
  • ビームサーチを実装する。
  • 初提出で38G。悪くない。
  • 最初と最後の接続部分を何とかしたい

12/3

  • 汚れの累積値という生スコアを評価していないことに気づいて評価する。
  • その他にも最適化をして41G。
  • 局所探索で改善できそうなことを確認する。

12/4

  • 未訪問箇所の挿入を改善する。
  • 44G。

12/5

  • 平日で唯一の休み。
  • 焼きなましの遷移をいくつか考える。
  • L字型変更と十字型解消を実装して45G。
  • バックトラック法といくつか最適化をして49Gで暫定2位。よし!

12/6, 12/7

  • ビームサーチで子孫の数を制限して悪化したのでやめる。
  • 過去改変機能をビームサーチに組み込んでも悪化するのでやめる。
    • 過去改変した情報を保持しておけば木上のビームサーチでも過去改変できる。
    • アイディアとしては面白いと思うが、うまくいかないので諦める。

12/8

  • プレテストの  N の分布を雑に確認したところ、 N が小さいものが多く、しかも  N が小さいものが得意なことに気づく。
  • 移動平均をビームサーチの評価関数に組み込んでスコアが上がる。

12/9

  • 朝から微熱が出ていて嫌な気分。
  • バックトラック法で関節点を2回通れるようにしたらスコアが上がった。
  • 午後に40度の熱が出て寝込む。発熱してもコンテストに参加するつもりだったが、40度までいくと自分には無理だった。
    • 調子いいときに発熱するのはやめてほしい。
    • まだパラメータ調整してないので明日はコンテストに取り組みたい。

12/10

  • 38度弱だったので問題なくコンテストに参加できた。
  • リファクタリングとパラメータ調整をやって暫定トップに立つ。
  • おそらく潜伏していた cuthbert さんに抜かれ、2位でフィニッシュ。
  • 最後やりきれない感じだったが悪くはなさそう。
  • 解説放送を見ずに寝る。

12/11

  • システムテストが終わって、ギリギリ bowwowforeach さんに勝っていて嬉しかった。
  • エラーが起きなかったのもよかった。
  • 解説放送を見て、 \sqrt{d} に比例した割合で訪問するなどの話を聞いて勉強になった。
  • 検査によりインフルエンザに罹患していたことがわかる。

感想

木上のビームサーチも、バックトラック法を遷移とする焼きなましも、過去のAHCで何度か実装していたので、非常にスムーズに実装することができました。過去問の復習を活かせたという実感があって非常に嬉しいです。

一方で、優勝するにはもう少し adhoc な考察や実装をする必要があった気がしていて反省もしています。

最後に

コンテストを主催してくださったフューチャー様、writer の wata さん、コンテスト参加者の皆様、そして最後まで読んでくださった読者の皆様、ありがとうございました。

日本橋ハーフマラソンでも対戦よろしくお願いします。


  1. ビームサーチの深さとは、シミュレーションを行うターン数を表すものとします。
  2. コンテスト終了後に気づきました。

AHC025 参加記

公式ビジュアライザ(seed = 5)

AHC025で優勝しました。長期AHCの過去最高順位は10位だったということもあり、非常に嬉しく思います。

順位表

はじめに

  • 解法だけでなく、考察や余談についても記載しました。
  • 重要度が低いと思われるものは折りたたみにしました。必要に応じてクリック(タップ)してください。
  • ベイズ推定やMCMCについては私自身が詳しくないため、コンテスト後に少し勉強して理解したことを書きました。「お気持ち」だけ汲み取っていただければ幸いです。
  • アイテム集合  A に対して、それらの重さ(の推定値)の和を  sum(A) と表記しました。行儀が悪い表現ですがお許しください。
  • 環境によっては数式が読みにくいと思います。下の記事などを参考にすると読みやすくなるかもしれません。

syleir.hatenablog.com

問題

公式の問題文と入力生成方法を参考にしてください。

atcoder.jp

この記事では  N, D, Q, \lambda などを問題文と同様に定義して使用します。

処理の流れ

  1.  (Q - N \log_2{N}) / N が大きいとき、クエリを消費しながらマージソートを行い、アイテムを重さ順にソートする。
  2.  D 個の集合を初期化する。
  3. 1~500個ほどの推定器(各アイテムの重さを推定するもの)を初期化する。
  4. 山登り法で改善する。具体的には以下の処理を繰り返す。
    1. 各推定器について、今までのクエリとの矛盾がなくなるように、反復法によって改善を行う。
    2. 時間に余裕がある場合、推定器を大きく更新し、再び反復法を適用する。
    3. 推定器の平均値を求める。
    4. 推定器の平均値を元に、2つの集合の間でいくつかのアイテムを入れ替える遷移でよさそうなものを列挙する。
    5. 推定器をサンプルとみなして各遷移のスコア改善量の期待値を求め、期待値が最も高いものを選択する。
    6. 選ばれた遷移によってスコアが改善するか否かをクエリで判定し、改善する場合にはその遷移を採用する。
    7. 残りの実行時間が十分短いとき、適当なクエリを投げてクエリを無駄に消費する。(TLE緊急回避モードです。基本的には実行されません。)
    8. クエリを合計で  Q 回行った場合、山登り法を終了し、解を出力する。

各処理について説明します。

アイテムのソート

 (Q - N \log_2{N}) / N がおよそ3.0以上のとき、マージソートでアイテムをソートしました。


閾値の設定について
マージソート N \log_2{N} 回弱のクエリを消費します。

ソートした後、山登り法のために  N の定数倍ほどのクエリを残したいという理由から、このような閾値を設定しました。

実際には  D が大きいときに閾値を上げていますが、改善の度合いはごくわずかです。


ソートを行う理由

後で重さを推定するときに、各アイテム間の順序関係がわかっていると非常に推定の精度が上がります。

例えば 0-indexed で  i 番目に軽いアイテムの重さを  - \frac{\log{(1 - (i + 0.5) / N)}}{\lambda} で推定するということが考えられます。


式の導出
指数分布について検索すると、累積分布関数が

 F(x; \lambda) = 1 - e^{-\lambda x}

であることがわかります。 F(x; \lambda)逆関数を求めます。

 y := F(x; \lambda)

 e^{-\lambda x} = 1 - y

 -\lambda x = \log{(1 - y)}

 x = - \frac{\log{(1 - y)}}{\lambda}

[0.0, 1.0] を  N 等分し、中央の値で代表すると  - \frac{\log{(1 - (i + 0.5) / N)}}{\lambda} が得られます。

 x = - \frac{\log{(1 - y)}}{\lambda} は一様分布から指数分布を生成することにも使用することができます。


seed = 1 の場合、次のような推定結果になりました。

ソートによる推定(seed = 1)

グラフにおいて、アイテムの真の重さの昇順にソートされ、重さの和が1になるように正規化されています。以降の推定に関するグラフにも同様のことが行われています。

比較回数の最小化

なるべく比較回数が少ないソートを選ぶことにより、山登り法でより多くのクエリを使用することができます。

私は比較回数が少なそうなマージソートを使用しました。


ソートの種類について
ソートの性能を評価する指標として次のようなものがあります。

  • 今回のコンテストで重要な指標
    • 比較回数
    • ソートを途中で止めたときの情報
  • 一般的な指標
    • 実行速度
    • 安定性

代表的なソートアルゴリズム

バブルソート

実装は簡単ですが、比較回数が  O(N^{2}) で、実行速度は遅いです。

クイックソート

一般的に速いとされていますが、安定性はありません。(安定性をもたせる実装だと遅くなると思います。)

理想的な場合の比較回数はあまり多くないですが、実際には理想的に動くとは限らないため、比較回数はやや多くなります。また、ピボットの選択が乱択でない場合、計算量が  O(N^{2}) になるケースを意図的に生成できます。

ソートを途中で止めた場合、いくつかのグループに分かれ、グループ間の大小関係がわかる状態になります。

軽いアイテムのソートを優先したり、逆に重いアイテムのソートを優先したりすることが可能です。

ヒープソート

小さいもの(あるいは大きいもの)から得ることができます。

ソートを途中で止めた場合、軽いアイテム(あるいは重いアイテム)の順序関係を正しく把握できます。

マージソート

安定性があり、比較回数は少ないです。

トップダウン方式とボトムアップ方式があります。

トップダウン方式は全体を半分ずつに区切る方式で、再帰関数を用いることで実装しやすくなります。

ボトムアップ方式は隣り合うグループをマージしていくような方式(Segment Tree をイメージするとわかりやすいかもしれません)で、通常は非再帰で実装します。

ボトムアップ方式はダブルバッファリングとの相性がよく、トップダウン方式よりも若干速い場合があります。

一方でトップダウン方式のほうが比較回数は少ないです。ボトムアップ方式だと  N = 33 のときに、最後に 32 のグループと 1 のグループに分かれ、単純にマージすると遅い可能性が高いです。

私が今回のコンテストで使用したのはトップダウン方式です。

途中で止めると、いくつかのグループに分かれ、グループ内でソートされた状態になります。クイックソートと対照的ですね。

merge-insertion sort

比較回数の最悪ケースが少ないらしいです。よく知りません。

ソートを途中で止める場合に、軽いアイテムと重いアイテムのどちらの情報が重要か

よくわからないというのが私の結論です。

直感的には重いアイテムのほうが影響が大きいため、重いアイテムの情報が重要そうに思えます。

一方で、山登り法で頻繁に動かすのは軽いアイテムであり、軽いアイテムの大小関係だけ把握するという考え方もあります。

ソートするかしないかの二択が離散的すぎて気持ち悪いという感覚があり、途中までソートすることを試したのですが、結局のところうまく使えませんでした。

ソート方法のまとめ

今回はマージソートか merge-insertion sort を使うのがいいと思います。


集合の初期化

ソートを行ったかどうかにかかわらず、アイテム  i をグループ  i \mod D に割り当てました

マージソートを行った場合、各アイテムの重さに関してある程度正しい推定ができ、初期解のコストを下げることが可能です。しかし、初期解をよくすると最終的なコストが増加しやすかったため、初期解は常に  i \mod D で設定しました。

局所探索の初期解はランダム性が高いほうがよいのでしょうか。私には理由がわかりませんでした。

(2023/10/30 追記: 初期解がよいとアイテムの入れ替え回数が減り、同一集合内のアイテム同士の比較が行われず、推定の精度が下がるからだと思います。)

推定

反復法

まず、今までのクエリ結果と矛盾しない(このことを制約と呼ぶことにします)重さを求めるための反復法について説明します。

推定の初期値は指数分布で生成し、ソートした場合には推定値もソートします。

以下に反復の流れを示します。

  • 次の処理を最大30回行う。
    1. 今までの全てのクエリ  l, r について次の処理を行う。
      1. 集合  l, r の和をそれぞれ求める。
      2. 求めた和とクエリの結果が矛盾している場合、重さの和が入れ替わるように  l, r の各要素の重さをスケーリングする(修正量をアイテムの重さに比例させる)。
      3. 重さが1.0未満のものは1.0にし、重さが 105 N / D 以上のものは 105 N / D にする。
    2. 更新が行われなかった場合、反復法を終了する。

推定値を更新する例

クエリの結果が等号だった場合、重さの和が等しくなるようにスケーリングしました。稀にしか起こらないと思います。


重さのスワップを行った経緯
元々は最急降下法によって重さを推定しようと考えていました。

不等式が満たされない場合に、重さの和が等しくなるように推定値を更新するときの修正量を「学習率が1.0である」と定義することにします。すなわち、学習率を  \alpha とすると、 sum(l) sum(l) + \alpha \frac{sum(r) - sum(l)}{2} になるように更新します。

α = 1.0 の例

学習率の値をいくらぐらいに設定するのがよいか考えます。

学習率が1.0未満の場合、満たされない制約が1つである場合などに更新を無限に繰り返すことになるため、不適切であると考えられます。(機械学習の文脈では基本的に学習率を1.0未満にしますが、今回は学習率の定義が異なることやデータが正確であることに注意してください。)

学習率が2.0より大きい場合、重さが負になる場合が考えられ、発散する場合もあるかもしれません。不適切です。

学習率は1.0以上2.0以下がよさそうです。パラメータ調整を行ったところ、学習率は2.0になりました。学習率2.0というのは重さの和を入れ替えるということになります。学習率を2.0にした場合、学習率を1.9にした場合よりもかなり性能がよく、学習率は敏感なパラメータだったと言えます。

連立一次方程式を解く際に、修正量を拡大して収束を加速させるという手法として、SOR法というものがあるらしいです。今回は制約が不等式のためSOR法ではありませんが、修正量を大きめに設定するという考え方は似ているかもしれません。ちなみにSOR法(かそれに近い何か)はトヨタコン決勝の懇親会でsaharanさんがAHC022(リクルートハーフマラソン)の文脈で教えてくださいました。ありがとうございます。


重さの和をスワップすることの利点

 sum(l) + sum(r) が不変のため、アイテムの重さの和が保たれます。重さの和が保たれない場合、クエリに含まれないアイテムの値を相対的に変更することになり、あまりよくない感じがします。

クエリの比較対象が要素1個同士の場合、単純に2つのアイテムの重さを入れ替えるだけなので、初期値が指数分布であれば指数分布が保たれることになります。分布が保たれるというのは非常によい性質です。実際には比較対象が1個同士とは限らないため、分布は少しずつ崩れると思います。

ベイズ推定の考え方

ソートを行った場合に、全ての制約を満たす重さとして次の2つの解を見つけたとします。

一様分布と指数分布

事前分布が指数分布であることを踏まえると、右側の推定のほうが確率が高そうです。

制約を満たす解の中でも、指数分布に近いものを見つけたいというモチベーションが生まれます。

理想的には解析的に事後分布を求めたいですが、おそらくそのような計算は困難です。そのためサンプリングを行って分布を近似することを考えます。

MCMCの考え方

サンプリングを行う度に指数分布から値を生成し(ソートした場合はソートも行い)、全ての制約を満たすもののみをサンプルとして扱うという方針が考えられますが、制約が増えるにつれて受理率が非常に低くなることが予想されます。

そこで、既に得られたサンプル(特に直前のサンプル)から新たなサンプルを生成して使用することにします。既存のサンプルから新たなサンプルを生成する方法として多くのバリエーションが考えられています。

サンプリング

並列に  \max{\left( 1, \frac{2 \times 10^{7} \times D}{N Q^{2}} \right)} 個ほどのサンプル(推定器)を用意します。


推定器の個数について
実行時間が2秒に収まる範囲で、なるべく推定器の個数を増やしたいです。

プログラムのボトルネックは、反復法で集合の和を求め、推定値を更新するところです。

集合の平均要素数はおよそ  N / D です。

クエリの数は  Q 回です。

制約の数は平均して  Q / 2 回程度です。

以上より、推定器の個数を  \frac{N Q^{2}}{D} に反比例させるとよさそうだということがわかります。


実行時間に余裕があるとき、直前のサンプルから新たなサンプルを生成します。

具体的には次のような手順です。

  1. 指数分布から  N 個の値を生成する。
  2. アイテムの大小関係が元のサンプルと等しくなるように生成した値を並べ替える。
  3. 反復法を適用し、制約を満たす解へ収束させることを目指す。

ソートを行っていた場合、アイテムの大小関係はサンプルによらないため、元のサンプルを無視して新たにサンプルを生成することと等価になります。

なお、実装上は反復法が収束したか否かはほとんど区別しませんでした。処理を収束したものに限定するよりもサンプル数を減らさないことのほうが重要だったようです。

ちなみに私よりも高得点だったwriterのwataさんは、ギブスサンプリングでサンプルを更新したそうです。ここで差をつけられた可能性もありますね。

サンプリング結果

サンプリング結果

ソートしない場合に、いくつかのアイテムの推定を大きく間違えていました。

山登り法

初期解を山登り法で改善します。

2つの集合の間でいくつかのアイテムを入れ替え、コストが減少したなら採用し、そうでないなら元に戻します。

コストが減少したか否かを判定する方法について述べます。

元の集合を  A \cup B, C \cup D とし、新たな集合を  A \cup D, C \cup B とします。 A, B, C, D の間に共通するアイテムはないものとします。

 sum(A \cup B) \lt sum(C \cup D) の場合、 sum(A) \lt sum(C) かつ  sum(B) \lt sum(D) ならばコストが減少します。 sum(A \cup B) \gt sum(C \cup D) の場合は不等号が逆になります。


最初の条件は必要か
コストが減少する必要十分条件

 sum(A) \lt sum(C) かつ  sum(B) \lt sum(D)

または

 sum(A) \gt sum(C) かつ  sum(B) \gt sum(D)

です。2回のクエリで判定することができ、 sum(A \cup B) \lt sum(C \cup D) で場合分けする方法よりも効率的に思えるかもしれません。

しかし、実際には場合分けを行ったほうがよいと思います。元のグループの大小関係はわかっている場合が多く、そのような場合は  sum(A \cup B) sum(C \cup D) の比較を省略できます。 sum(A \cup B) \lt sum(C \cup D) sum(A) \gt sum(C) の場合に、 sum(B) sum(D) の比較を省略でき、比較が1回で終わることがあります。

最大3回の比較を要しますが、集合間の大小関係を把握できるメリットもあるため、集合同士の比較も行ったほうがよいと考えられます。

不採用のときに比較を1回で済ませるために、 sum(A) \lt sum(C) sum(B) \lt sum(D) のうち、成り立ちにくそうなほうから比較を行いました。


クエリの省略

結果がわかる場合はクエリを省略します。省略するパターンを以下に示します。

  • 片方の集合が空のとき
  • 既に同じクエリを投げているとき
  • 集合の比較で、推移律から結果がわかるとき

集合の大小関係が  sum(A) \lt sum(L) \lt sum(B) \lt sum(R) \lt sum(C) となっていて、集合  L, R を変更してコストが下がったとします。このとき、 sum(L), sum(R) に関して以下の関係が保存されます。

  •  sum(A) \lt sum(L)
  •  sum(A) \lt sum(R)
  •  sum(L) \lt sum(C)
  •  sum(R) \lt sum(C)

コンテスト上位の方の中には、より複雑な比較の省略を行った人もいたようでした。

入れ替え候補の列挙

2つのグループを選択し、重さの推定値が与えられるものとすると、partition problem に帰着されます。partition problem はNP完全で、擬多項式時間アルゴリズムヒューリスティックな手法が存在します。

私は半分全列挙をベースとして partition problem を解きました。理想的には厳密解が得られますが、集合の要素数が多いときは長い時間がかかります。そこで、各集合の入れ替え可能なアイテムを、重さが小さいものを中心に 2 ~ 10 個に限定しました。


「重さが小さいものを中心に」の補足
より正確には、元の集合における和の差の推定値の 1/4 に近い重さのアイテムを優先して選びました。集合の差はすぐに小さくなるため、重さが小さいものを優先しているのとあまり変わりません。

優先度には乱数をかけました。

具体的な優先度の付け方を以下に示します。

  • アイテムの重さの推定値を  w とする。
  • 優先する重さを  \Delta := \frac{|sum(l) - sum(r)|}{4} とする。
  • 乱数の大きさを調節する変数(温度)を  t:= \frac{N Q'}{23000 \times D} とする。ただし  Q' は残りのクエリ数を表すものとする。
  • 優先度を  rand\_double(0, t) + | \log_2{\frac{w}{\Delta}} | で定義し、優先度が小さいものを優先して選ぶ。


入れ替え可能なアイテム数について
 Q が大きいほど入れ替え可能なアイテム数を増やしました。

 Q が小さいときは推定の精度が悪かったからか、変更を小さくしたほうがスコアが上がりました。


半分全列挙の実装上の工夫
全列挙パートでは要素を1つずつ追加します。その際、マージ手法を使うことによりソート関数を呼ばずにソートできます。また、復元用に全ての配列は削除せずに保持します。

2つの全列挙された配列はソートされているため、尺取法の要領で最適な値を計算できます。要するに  O(2^{n / 2} n) ではなく、 O(2^{n / 2}) で計算できます。

解の復元は、全列挙パートで追加したときと逆の順序で行います。和が前の配列にあればその要素を使う必要がなく、そうでないときは使う必要があります。配列がソート済みであることに注意すると、値が含まれるか否かは二分探索でわかるため、復元は  O(n^{2}) で行えます。

数値誤差に対処するのは面倒なため、重さは全て整数として扱いました。

関連する AtCoder Beginner Contest の問題 atcoder.jp

私の提出コード atcoder.jp


3-partition problem
2-partition problem へ帰着すると解きやすい問題になりますが、本来解くべき問題は D-partition problem です。

2-partition problem しか解かない場合、山登り法の局所解に陥るケースがありそうだったため、3-partition などのより多様な近傍を使用することをコンテスト中に考えました。しかし、私は 3-partition の近傍でスコアを上げることができませんでした。

うまくいかなかった理由として以下のことが考えられます。

  • そもそも 2-partition の近傍で十分よい解が得られる。
  • 半分全列挙をベースとする方法で 3-partition problem を解くことに無理があった。Largest Differencing Method などの別の手法を使うべきだった。(Largest Differencing Method はコンテスト後に知りました。)

入れ替え候補の選択

以上のことを踏まえて、次のような流れで入れ替え候補の選択を行います。

  1. サンプルの平均値を計算する。
  2. 2つの集合の組を全列挙し、次の処理を行う。ただし、時間がないときは集合の和の差が大きいペアについてのみ処理を行う。
    1. サンプルの平均値を推定値とし、半分全列挙ベースの手法でよさそうな入れ替えを求める。
    2. 入れ替えることによる改善の期待値をサンプルから(擬似的に)求める。
  3. 期待値が最も高い入れ替えを求める。

全てのサンプルが制約を満たすとき、サンプルの平均値も制約を満たします。これは制約に線形性があることから言えます。


線形性の補足
具体例で説明します。

 x + y \lt z

という制約があり、

 (x, y, z) = (x_1, y_1, z_1), (x_2, y_2, z_2)

が制約を満たすとします。

このとき、

 x_1 + y_1 \lt z_1

の両辺を2で割ると

 \frac{x_1}{2} + \frac{y_1}{2} \lt \frac{z_1}{2}

が得られます。よって

 (x, y, z) = (\frac{x_1}{2}, \frac{y_1}{2}, \frac{z_1}{2})

も制約を満たします。より一般的に、正の定数をかけても制約を満たすということが言えます。

さらに

 x_1 + y_1 \lt z_1 x_2 + y_2 \lt z_2 の両辺を足すと

 (x_1 + x_2) + (y_1 + y_2) \lt z_1 + z_2

が得られます。よって

 (x, y, z) = (x_1 + x_2, y_1 + y_2, z_1 + z_2)

も制約を満たします。

解に正の定数をかけたり解同士の加算を行ったりしたものも解になることから、平均や加重平均をとってもよいということがわかります。

また、この記事においてグラフを描画するときに重さの和が1になるように正規化している理由は線形性があるからです。


期待値の求め方
基本的には、全てのサンプルについて  \max{(0, 改善量)} を計算し、その平均値をとればよいです。

私の実装では、サンプルの平均もサンプルとして扱っています。サンプル数を  n として、サンプルの平均は  \sqrt{n} / 2 個分のサンプルとして期待値を計算しました。加重平均の重みが普通のサンプルの  \sqrt{n} / 2 倍になっているということです。

サンプルの平均もサンプルとして扱うことで、サンプル数が少ないときに正則化効果が生まれたと推測しています。


その他

集合の表現方法

3つの方法を紹介します。

配列

アイテムの番号を整数として配列に保持します。最も単純な表現方法だと思います。

C++だと std::vector<int> などです。

ソートしない場合、要素の追加は  O(1) ですが、要素の削除や含まれるかの判定が  O(n) で遅いです。

Index Set

tomerunさんのブログを参考にしてください。

topcoder-tomerun.hatenablog.jp

今回の問題では  N \leq 100 であるため、平衡二分木などよりも高速で便利だと思います。

bitset

C++における std::bitset<100> です。

次のような利点があります。

  • コピーが高速である。
  • 集合の演算をビット演算で簡潔に記述でき、処理が高速である。

使い分け

配列と bitset を使用しました。イテレーションを回すことが多い場所では配列を使用し、そうでないところでは bitset を使用しました。


bitset を配列に変換する方法

C++20 の std::countr_zero などを使用することで、要素数の線形時間で変換できます。定数倍が軽くはなく、 N 回ループを回したときと速度があまり変わらないかもしれません。

vector<int> bitset_to_vector(const bitset<100>& l) {
    ull low = (l & bitset<100>(ULLONG_MAX)).to_ullong();
    ull high = (l >> 64).to_ullong();

    vector<int> ret;
    ret.reserve(l.count());

    while (low) {
        int v = countr_zero(low);
        ret.push_back(v);
        low ^= (1ull << v);
    }
    while (high) {
        int v = countr_zero(high);
        ret.push_back(v + 64);
        high ^= (1ull << v);
    }
    return ret;
}


順位表の使い方

ACを一部のテストケースに制限し、自分が苦手なタイプのケースを特定していました。

assert(Q < 8 * N);

コンテスト終盤は  Q / N が小さいケースに弱いことがわかっていたため、 Q / N が小さいケースを中心に改善点を探りました。

結局のところ、最後まで  Q / N が小さいケースに弱かったため、サンプリングをうまく使えていなかったのかもしれません。

スコアの評価方法

各テストケースのコストの対数を求め、その和を最小化しました。

パラメータ調整

大雑把なパラメータ調整は実装したときに行い、細かいパラメータ調整はコンテスト最終日に行いました。

 N, D, Q によって決定されるパラメータ調整は、関数の探索になるのでかなり面倒でした。

パラメータ調整は全て手動で行いました。Optuna のような自動最適化ツールは使っていません。

高速化の寄与

高速化によってどれだけコストが下がるかは、実行時間を長くすることにより、高速化を行う前に推測できます。

私の方針だと高速化は非常に重要で、参照渡しにし忘れていた部分を直しただけで0.8%ほど順位表のスコアが上がったことがありました。

提出コード

atcoder.jp

最後に

いくらか考察も書きましたが、そのほとんどは実装後やコンテスト後に考えたことです。実際には試行錯誤の繰り返しでした。

コンテストを開いてくださったAtCoder様、writerのwataさん、コンテスト参加者の皆さん、最後まで読んでくださった読者の皆さん、ありがとうございました。

AHCで赤になりました

2023/08/23 に終了した AHC022 をもって AtCoder の Heuristic Rating が 2801 となり、レッドコーダーに昇進しました!

この記事では、今までの出来事や頑張ったことなどを述べます。レーティングを上げたい人やコンテストで勝ちたい人にとって、少しでも参考になれば幸いです。

コンテスト実績など

コンテスト実績

コンテスト成績表

短期コンテストの延長戦順位(記事執筆時点)

コンテスト 延長戦順位(本番相当) 延長戦順位(他者の延長戦を含む)
Chokudai Contest 005 1 1
AHC002 3 18
AHC005 2 10
AHC006 10 46
AHC007 1 6
AHC009 10 30
AHC010 1 1
AHC012 1 5
AHC015 1 1
HACK TO THE FUTURE 2023 本選 1 2
AHC020 1 7
AHC021 1 9
新ジャッジテストコンテスト -Heuristic- 1 1
TOYOTA Programming Contest 2023 Summer final 1 1

主な出来事

年/月 出来事
2020/03 高校を卒業する
2020/03 プログラミングを始める(言語は Python
2020/05 Algorithm で水色になる
2020/11 Algorithm で青色になる
2021/03 AHC001 で黄色相当のパフォーマンスを取る
2021/06 Algorithm で黄色になる
2021/09 RECRUIT 日本橋ハーフマラソン 2021 で初めて賞金を獲得する
2021/11 Heuristic で黄色になる
2022/01 C++ の勉強をする
2022/06 AHC011 で初めて橙色相当のパフォーマンスを取る
2022/10 AHC015 で優勝して橙色になる
2023/08 Heuristic で赤色になる

AHC015 で優勝し、その他に橙色相当のパフォーマンスを7回取って赤になりました。rated の AHC に全参加しているにもかかわらず、赤色相当のパフォーマンスを1回しか取っていないので、レッドコーダーの中ではかなり実力が低いと思います。

レッドコーダーの中には AHC001 からいい成績を残している方も多くいらっしゃいますが、私の場合は AHC の参加や復習を通じて実力をつけました。AHC が始まったころはヒューリスティックコンテストの経験が乏しく、黄色相当のパフォーマンスを出すだけで精一杯でした。なぜか AHC014 のあたりから橙色相当のパフォーマンスを高確率で出せるようになり、レーティングが大きく上がりました。突然いい結果を出せるようになるというのは面白い現象ですね。

長期コンテストの色々

スケジューリング

コンテストに集中できるように、コンテスト中は大学の課題をできるだけやらないようにしました。コンテスト前に片付けられる課題は終わらせ、コンテスト中に出た課題はできるだけコンテスト後に回しました。

一般的に社会人は忙しいはずなので、長期コンテストで上位に入る社会人は凄いとつくづく思います。

洞察

日常生活の中でもコンテストのことを考えている気がします。食事、風呂、睡眠中など様々なタイミングで洞察を得ています。

特に睡眠中は無意識状態で筋のよい考察をしているような気がしています。そのため睡眠時間はあまり削らないです。

疲労

AHC が始まったころは頑張りすぎてつらかったです。朝5時に目覚めたり、座りすぎて腰を痛めたり、キーを打ちすぎて指を痛めたりしました。コンテストが終わった翌日に疲れが押し寄せて何もやる気が起きなくなることもありました。軽度の燃え尽き症候群だったかもしれません。今振り返ると、質の悪いアイディアをひたすら実装していた気がします。

最近は以前よりもリラックスして取り組んでいます。コンテスト中だと睡眠の質が少し落ちるような気がしますが、生活に支障が出ることはほとんどないです。

実力向上に寄与した(する)と思うこと

コンテストへの参加

本番だからこそ味わえる緊張感、悔しさ、嬉しさがあると思います。特に負けたときの悔しさは成長につながる気がします。

解説や参加記を読む

コンテスト終了後、Twitter (X)、ブログ、解説放送などで他の人(特に上位陣)の解法を見るようにしています。自分の知らない手法、考察の仕方、自分の解法との差など、非常に多くのことを学べます。

延長戦

コンテスト終了後や暇なときに短期 AHC を復習することがあります。基本的にはコンテスト上位の解法を参考にして実装します。

短期 AHC は高いスコアを比較的簡単に出せるので復習にお勧めです。

アルゴリズム

Algorithm のコンテストで強い人は、典型的なアルゴリズムを知っているだけでなく、アルゴリズムの使い方が上手かったり、実装速度が速かったりする傾向があり、Heuristic コンテストでも役立つと思います。

Algorithm も強くなりたい!

コンピュータサイエンス

機械学習はときどき使えそうですね。(私は AHC で機械学習を使うほど詳しくはないです。)

コンピュータアーキテクチャ、OS、計算量理論なども少しはコンテストに役立つかもしれません。

プログラミング言語

AHC で橙以上を目指す場合、C++ や Rust などの高速な言語を使えるようにしておくことを強く勧めます。焼きなまし法、ビームサーチ、モンテカルロ法など、処理が速いほどスコアが上がる(ことが多い)手法がよく使われるからです。

私はコンテスト後に Python (PyPy) のコードを C++ に書き換えただけで大きくスコアが伸びたことがあり、とても悔しい思いをしました。それ以降は C++ を使うようにしています。

ただし、低速な言語を使用することについて否定するつもりはありません。一部のインタラクティブ問題では、実行時間をほとんど使わずに高得点が出ることもあります。

謝辞

AtCoder および writer の方へ

コンテスト環境および質の高い問題をいつも提供していただきありがとうございます。

コンテスト主催企業様へ

コンテストを主催していただきありがとうございます。おかげさまで Amazon での買い物には(今のところ)困らないです。

コンテスト参加者へ

競う相手がいるというだけでもありがたいことですが、解法を共有してくださると復習しやすくて大変ありがたいです。私が知っているヒューリスティックの手法の大部分はコンテスト参加者から学びました。

最後に

入赤は通過点で目標は世界一です。

TOYOTA Programming Contest 2023 Summer final 延長戦の解法

seed = 0, score = 997530864

TOYOTA Programming Contest 2023 Summer final における 297,762,906,271 点の解法を紹介します。延長戦の解法であり、実装に2日ほど要しました。コンテストは3時間30分しかないため、この解法をコンテスト中に実装しきるのは難しいと思います。

優勝者である terry さんの解説記事と一部内容が被ります。あらかじめご了承ください。

www.terry-u16.net

前提知識

木上のビームサーチについては rhoo さんの記事が参考になると思います。

qiita.com

問題概要

 D \times D マスの倉庫にコンテナを運び込み、その後コンテナを運び出します。

運び込み・運び出しでは、出入り口から空きマスを辿って到達できる空きマス・コンテナにしか操作が行えません。

コンテナを運び込む順番はインタラクティブに与えられます。

運び出すコンテナ番号の転倒数を最小化してください。

詳細は公式の問題文を参考にしてください。

atcoder.jp

制約

 D = 9

解法の概要

運び込みパートと運び出しパートで分けて考えます。

運び込みパートではモンテカルロ法を使用し、評価関数に従った貪欲法を改善します。

運び出しパートでは木上のビームサーチを使用します。

運び込みパート

貪欲法

各コンテナがどのあたりに配置されるとよいかを定めておき、その周辺に運び込むことを考えます。

番号が小さいコンテナから運び出したいので、番号が小さいコンテナを出入り口の近くに、番号が大きいコンテナを出入り口から遠いところに置くとよいと考えられます。そこで、幅優先探索を用いて出入り口からの各マス  i までの距離  a_i を求め、距離の配列  a をソートしたものを  b としたとき、コンテナ  t は出入り口からの距離が  b_t 程度になるマスに置くのがよさそうです。

また、空きマスの連結性を確保するために、端のマスからコンテナを詰めていくとよさそうです。

そこで、次のような貪欲法が考えられます。

  1. 運び込むコンテナ番号  t を読み込む。
  2. low link で関節点を列挙する。
  3. 各非関節点  i について評価値  f(t, i) を計算し、評価値が最も高いマスにコンテナ  t を運び込む。評価値  f(t, i) の定義は以下の通り。

 f(t, i) = - (2 \times |a_i - b_t| + (iと隣接する空きマスの数))

ここで簡単な考察をします。運び出すときの転倒数を0にするにはコンテナ0は出入り口と接している必要があります。一方で、番号が最大のコンテナは出入り口から遠いところにあるとよいというだけで、置ける場所の自由度は比較的大きいです。自由度の違いに着目すると、番号が大きいコンテナを出入り口から遠くに置くことよりも、番号が小さいコンテナを出入り口の近くに置くことを優先するとよさそうです。

そこで、コンテナ  t を運び込むときはコンテナ  t + 3 の位置に運び込むようにします。すなわち、評価関数の  b_t b_{t + 3} に修正します。ただし、コンテナ0についてはコンテナ0の位置に運び込むようにします。

出入り口から近いマスを3つほど確保しておき、必要になったときに使うというようなイメージです。

もしかしたら、出入り口から遠いところに置くことを優先し、空きマスの連結性を確保しやすくなるというメリットもあるかもしれません。

モンテカルロ法

上記の貪欲法をモンテカルロ法で改善しました。アルゴリズムは以下の通りです。

  1. 評価値が高いマスを候補点として4つ選ぶ。
  2. 一定回数または一定時間のプレイアウトを行う。プレイアウトは以下の手順。
    1. 候補点にコンテナ  t を運び込む。
    2. 残りのコンテナの順番をランダムに生成する。
    3. 上記の貪欲法で運び込みをシミュレートする。
    4. 貪欲に運び出しを行い(後述)、スコアを求める。
  3. 平均スコアが最も高いものを採用する。

補足

  • 貪欲に運び出すとは、運び出しが可能なコンテナのうち、番号が最小のものを取り続けることを表すものとする。優先度付きキューを用いて時間計算量  O(D^{2} \log{D}) で計算可能。(Dijkstra 法とよく似ている。)
  • プレイアウトの乱数列は全ての候補点で同じものを使用する。同じ条件でプレイアウトを行わないと比較しづらい。
  • 残りのコンテナ数が6個以下の場合、順番を全探索する。厳密なプレイアウトの期待値が求まる。

各ターンのプレイアウト回数を揃えた場合、毎ターン1000回ほどのプレイアウトを行うことができました。候補の数が4つの場合、各候補について250回ほどのプレイアウトが行われます。

二次元配列を平坦化して一次元配列で置き換えるなど、定数倍高速化に注意して実装しました。(low link の結果を保存し、同じ盤面が現れたときに結果を使い回すというようなことも行いましたが、あまり速くならなかったです。)

運び出しパート

ビームサーチをします。履歴のコピー回数を減らすために木上のビームサーチを使用しました。

既に運び出したコンテナの集合と、直後に運び出せるコンテナの集合を状態としてもち、差分計算を行います。

評価関数は以下のように設定しました。

 -2 \times (確定転倒数) + (直後に運び出せるコンテナの数)

確定転倒数とは、(連結性を無視して)未処理のコンテナを昇順に運び出した場合の転倒数のことです。既に運び出したコンテナの集合を bitset などで管理することにより、ビットシフトと popcount を用いて確定転倒数の差分計算を高速に行うことができます。

基本的には直後に運び出せるコンテナの数だけ、1つのノードから分岐します。ただし、番号が最小のコンテナを取り出せるときは、それが最善手であることから、それ以外の遷移を枝刈りしました。

既に運び出したコンテナの集合が一致するノードが複数出現した場合、確定転倒数が最小のものだけを残し、他のノードを削除します。既に運び出したコンテナの集合を bitset で管理することで、重複ノードの検出を厳密に行えます。(通常のビームサーチではハッシュ値を使用して重複探索を削除するのですが、今回はハッシュ値を使わなくてもよいということです。)

ビーム幅は1000にしました。ビームサーチの実行時間は0.1秒程度です。モンテカルロ法に時間を使うのが得策だと判断しました。

スコアへの寄与

プログラムを一箇所変えたものをいくつか用意し、それぞれについて100ケースの平均得点率と平均実行時間を測定しました。

変更点 平均得点率(%) 平均実行時間(秒)
なし 99.400 1.94
運び込みを貪欲法にする 98.221 0.11
モンテカルロ法のプレイアウト回数を半分にする 99.327 1.02
運び込みの評価関数において+3を行わない 98.992 1.95
運び出しを貪欲法にする 99.375 1.85

ビームサーチの寄与は比較的小さいことがわかりました。

使用した典型まとめ

今までのコンテストで身につけた典型的な手法をいくつか使用しました。典型手法とそれを使える代表的なコンテストを列挙しておきます。

考察

  • 複数の問題に分けて考える。(完全に分離するとよくない場合もある。)
    • AHC011
  • オブジェクトの自由度の違いに着目する。
    • AHC021

手法

  • 貪欲法
    • 評価関数を工夫する。
  • モンテカルロ法
    • 貪欲法でシミュレーションする。
    • 終盤は全探索して厳密な期待値を求める。
    • AHC015
  • ビームサーチ
    • 木上のビームサーチ
      • AHC021
    • 重複探索の削除
    • 明らかに無駄な探索を枝刈りする。
    • 評価関数は確定スコアに状態のよさを加えたものにする。

高速化

  • 多次元配列の平坦化
  • ビット演算

アルゴリズム

提出コード

木上のビームサーチはライブラリ化したものを使用しています。

atcoder.jp

最後に

延長戦なら満点を取れるかもしれないと思いながら復習しましたが、全ケース満点は無謀な試みだったようです。ただ、何百ケースか試すと1ケースほど満点が取れるみたいです。

(2023/8/30 追記: 解説放送を見て、思ったより満点を取れるらしいことがわかりました。評価関数あたりに改善の余地がありそうです。)

オンサイトコンテストを主催してくださったトヨタ自動車様、AtCoder の皆様、競技者の皆様、そして読者の皆様、ありがとうございました。