Yuulis.log

Yuulis.log

トンネルを抜けるとそこは参照エラーであった。

【AtCoder】Introduction to Heuristics Contest - AtCoder Contest Scheduling (上位2%解法) | 水コーダーが挑むAHC

AtCoder-Heuristics系 AtCoder-貪欲法 AtCoder-山登り法 AtCoder-焼きなまし法 AtCoder Heuristics系 貪欲法 山登り法 焼きなまし法 競技プログラミング C++

atcoder.jp

実行時間制限: 2 sec / メモリ制限: 1024 MiB


本コンテストは、現在定期的に開催されている AtCode Heuristics Contest (AHC) の入門向けのものである。

通常の Algorthm 系コンテストとは異なり、最適解を出すのが難しい問題に対して、制約の範囲内で出来るだけ良い解を作成するものである(過去に「マラソン型」と呼ばれていた、最適化問題を解く形態)。

本記事では、「延長戦順位表」において上位2%以上の成績を出すことを目標とした解法について考察していく。


問題概要

タイプ 1 からタイプ 26 と番号付けられた 26 種類のコンテストについて、ユーザーの満足度が出来るだけ高くなるように、 D 日分のコンテストの日程を決定することを考える。毎日ちょうど一つのコンテストを開催し、各コンテストはその日のうちに終了する。

ここで、満足度は以下のように算出される。

  • 満足度は負の値も取り得る。
  • 1日目の開始時点における満足度は 0 である。
  • d 日目にタイプ i のコンテストを開催した場合、満足度は s_{d,i} 増加する。
  • コンテストのタイプごとに満足度の下がりやすさ c_i が定まっており、各 d=1,2,\dots,D 日目の終わりに以下のように満足度が低下する。
    • d 日目以前にタイプ i のコンテストを開催した最後の日を \mathrm{last}(d,i) とする。ただし、まだ一度もタイプ i のコンテストを開催していない場合は 0 とする。
    • d 日目の終わりに、 \displaystyle \sum_{i=1}^{26} c_i \times (d-\mathrm{last}(d,i)) だけ満足度が低下する。

最終的な満足度が出来るだけ大きくなるようなコンテストの日程を求めよ。

得点

D 日目の終了時点での満足度を S として、 \max(10^6 + S, 0) の得点を得られる。

制約

  • D = 365 で固定。
  • 各 c_i は整数値で 0\leq c_i \leq 100 を満たす。
  • 各 s_{d,i} は整数値で 0\leq s_{d,i} \leq 20000 を満たす。

考察

とりあえず正の得点を得る

AHC でまずやるべきなのは、正の得点を得ること。

一見ばかばかしく思えるかもしれないが、問題内容が正しく把握できているかを確認するためにも重要なことである。


ここでは簡単に、開催するコンテストタイプを1日目から順に 1 \rightarrow 2 \rightarrow \dots \rightarrow 26 \rightarrow 1 \rightarrow \dots のように繰り返していく方針で解を出力してみる。

提出コード(クリックで展開)

#include <bits/stdc++.h>
using namespace std;

// ======================================== //

struct Input
{
    /// @brief 日数(D = 365 で固定)
    int length;
    /// @brief 各コンテストの満足度の下がりやすさ
    vector<int> satisfaction_decrease;
    /// @brief ある日の各コンテストの満足度の増加値
    vector<vector<int>> satisfaction_increase;

    /// @brief 入力データの読み込み
    /// @return 入力データ
    static Input read()
    {
        int D;
        cin >> D;
        vector<int> c(26);
        for (int i = 0; i < 26; ++i)
        {
            cin >> c[i];
        }
        vector<vector<int>> s(D, vector<int>(26));
        for (int i = 0; i < D; ++i)
        {
            for (int j = 0; j < 26; ++j)
            {
                cin >> s[i][j];
            }
        }

        return Input{D, c, s};
    }
};

struct Output
{
    /// @brief ある日に開催されるコンテストのタイプ
    vector<int> contest_types;

    /// @brief 解を出力する
    void print()
    {
        for (int i = 0; i < (int)contest_types.size(); i++)
        {
            cout << contest_types[i] + 1 << endl;
        }

        return;
    }
};

Output solve(const Input &input)
{
    /// @brief contest_types[d] := d 日目に開催されるコンテストのタイプ
    vector<int> contest_types(input.length);

    for (int i = 0; i < input.length; i++)
    {
        contest_types[i] = i % 26;
    }

    return Output{contest_types};
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    Input input = Input::read();

    Output output = solve(input);

    output.print();

    return 0;
}

これを提出すると、 13,639 点 を得ることができた。大半のseedでは全体 Score は 0 だが、運が良いと正の全体 Score が得られるようだ(seed = 7とか)。

atcoder.jp

満足度(スコア)を計算する関数の作成

具体的な解法を考える前に、もう少し準備をしたい。

この問題では最終的に D 日分のコンテスト日程を求めることになるが、決定したコンテスト日程を定量的に評価する指標が欲しいのである。

本問では、問題文にもある D 日目終了時点での満足度(以下、単に「スコア」と言うこともある)を用いればよいだろう。満足度が高いほど、より良い日程ということになる。

この満足度の計算をC++コードに落とし込むと、以下のようになる。

/// @brief 満足度(score)計算用の関数
/// @param input 入力データ
/// @param contest_types 各日のコンテストのタイプを表す配列(長さ d_l)
/// @return d_l 日目時点での満足度
ll evaluate(const Input &input, const vector<int> &contest_types)
{
    ll score = 0;
    vector<int> last(26, 0);

    for (int d = 0; d < (int)contest_types.size(); d++)
    {
        last[contest_types[d]] = d + 1;

        for (int i = 0; i < 26; i++)
        {
            score -= input.satisfaction_decrease[i] * ((d + 1) - last[i]);
        }

        score += input.satisfaction_increase[d][contest_types[d]];
    }

    return score;
}

この部分の内容は、本コンテストのB問題に当たる。特に計算量を意識することもないので、最近の ABC だと250点問題くらいだろうか。

簡単な貪欲解を作る

本問では、各 d = 1, 2, \dots, D 日目それぞれについて、開催するコンテストを逐次的に決定可能である。しかも、決定済みの部分的な日程について、その満足度も計算可能だ。

したがって、以下の単純な貪欲的な解法が成立する。

  • 各 d = 1, 2, \dots, D 日目には、 d 日目終了時点での満足度が最も高くなるようなコンテストタイプを開催する。

C++コードに落とし込むと、solve関数内は以下のようになる。

solve関数のコード(クリックで展開)

Output solve(const Input &input)
{
    /// @brief contest_types[d] := d 日目に開催されるコンテストのタイプ
    vector<int> contest_types;

    for (int d = 0; d < input.length; d++)
    {
        /// @brief d 日目に開催するコンテストとして最良なもの
        int best_type = -1;
        /// @brief 現時点での最良な満足度
        ll best_score = -INFL;

        for (int i = 0; i < 26; i++)
        {
            contest_types.push_back(i);

            ll current_score = evaluate(input, contest_types);
            if (current_score > best_score)
            {
                best_score = current_score;
                best_type = i;
            }

            contest_types.pop_back();
        }

        // d 日目に開催するコンテストのタイプを最良のものに確定
        contest_types.push_back(best_type);
    }

    return Output{contest_types};
}

これを提出すると、 62,634,806 点(延長戦 866 位)を取ることができた。ちなみに、seed = 0の全体 Score は 1,254,667 である。

atcoder.jp

順位表はこの得点で一つの団子状態(200人超！)となっており、とりあえず最初の一歩は踏み出せた形である。

山登り法

とりあえず貪欲解を構築できたので、この解をさらに改良していきたい。

ここで重要なのが、先述した「開催するコンテストを逐次的に決定可能」であるという点。

つまり、既存の解の一部を少し変化させることを制限時間いっぱいまで繰り返して、より良いスコアを目指すことができるということだ。もちろん、スコアが悪くなれば元に戻す。

このような解の探索手法を山登り法という。


さて、問題は「どのように解を変化させるか」(近傍の取り方という)だが、ここでは以下の2種類を試してみたい。

一点更新

  • D 日間の中からランダムに1日を選び、その日に開催するコンテストタイプをランダムに変更する。

これは比較的簡単に実装できる。先ほどの貪欲解をgreedy_solve()で構築し、それを初期解として一点更新を時間いっぱいまで繰り返す。

C++での実装コードは以下の通り。

solve関数のコード(クリックで展開)

Output solve(const Input &input)
{
    /// @brief contest_types[d] := d 日目に開催されるコンテストのタイプ
    /// @note 貪欲法で求めた解を初期解とする
    vector<int> contest_types = greedy_solve(input).contest_types;
    /// @brief 現時点での満足度
    ll current_score = evaluate(input, contest_types);

    /// @brief 乱数生成器
    mt19937 rand{42};
    /// @brief 山登り法の開始時刻
    auto start_time = chrono::system_clock::now();
    /// @brief 制限時間(ミリ秒単位)
    const int LIMIT = 1900;
    /// @brief 試行回数
    int iteration = 0;

    while (true)
    {
        auto current_time = chrono::system_clock::now();

        // 制限時間を超えたら終了
        if (chrono::duration_cast<chrono::milliseconds>(current_time - start_time).count() >= LIMIT)
        {
            break;
        }

        // 変更する日数とコンテストのタイプをランダムに選ぶ
        int d = rand() % input.length;
        int type = rand() % 26;

        // 元のコンテストのタイプを保存
        int old_type = contest_types[d];

        // d 日目のコンテストのタイプを変更
        contest_types[d] = type;

        // 満足度を計算
        ll new_score = evaluate(input, contest_types);

        if (current_score < new_score)
        {
            cerr << "iteration: " << iteration << ", new score = " << new_score << endl;
            current_score = new_score;
        }
        // 満足度が改善されなかった場合, 元に戻す
        else
        {
            contest_types[d] = old_type;
        }

        iteration++;
    }

    cerr << "total iterations: " << iteration << endl;

    return Output{contest_types};
}

これを提出すると、 107,379,543 点(延長戦 251 位)を取ることができた。seed = 0の全体 Score は 1,994,060 である。

atcoder.jp

先ほどの貪欲解から飛躍的に全体 Score を改善することができた。


なお、このコードでは、スコアが更新されたときのイテレーションと新しいスコア、そして山登り法全体で回ったイテレーションを標準エラー出力に出力するようにしている。

こうすることで、山登り法においてどのくらいの更新が行われているのかが分かるようになる。

ちなみに、seed = 0での結果はこんな感じ。

iteration: 7, new score = 257914
iteration: 29, new score = 263844
iteration: 56, new score = 275447
...
iteration: 57489, new score = 1026615
iteration: 73499, new score = 1027139
iteration: 73725, new score = 1036725
iteration: 78792, new score = 1036744
total iterations: 327759

全体では約32.7万回ものイテレーションが回っているのに対し、実際にスコアが更新されているのは7.8万回までで、それ以降は更新されていないようだ。

どうやらこの近傍の取り方では、探索できる解空間は限られてしまっているらしい...

一点更新 + 二点 Swap

そでは、一点更新に加えて以下のような「二点 Swap」による近傍の取り方も試してみよう。

  • D 日間の中で比較的近い2つの日をランダムに選び、それらコンテストのタイプを入れ替える。

ここで、「比較的近い」とはどこまでを指すのかについてだが、これは実験して探すしかないのだ。

色々な値を入れてスコア計算を繰り返すと、私の場合は「ある日から14日以内の別の日を選んで入れ替える」ことが最善だと分かったので、実装コードはその方針のものを載せている。


また、二点 Swap を行うだけでは部分的にしか探索できないので、確率で先ほどの一点更新も行うようにする。

以上をC++で実装すると以下のようになる。

solve関数のコード(クリックで展開)

Output solve(const Input &input)
{
    /// @brief contest_types[d] := d 日目に開催されるコンテストのタイプ
    /// @note 貪欲法で求めた解を初期解とする
    vector<int> contest_types = greedy_solve(input).contest_types;
    /// @brief 現時点での満足度
    ll current_score = evaluate(input, contest_types);

    /// @brief 乱数生成器
    mt19937 rand{42};
    /// @brief 山登り法の開始時刻
    auto start_time = chrono::system_clock::now();
    /// @brief 制限時間(ミリ秒単位)
    const int LIMIT = 1900;
    /// @brief 試行回数
    int iteration = 0;

    while (true)
    {
        auto current_time = chrono::system_clock::now();

        // 制限時間を超えたら終了
        if (chrono::duration_cast<chrono::milliseconds>(current_time - start_time).count() >= LIMIT)
        {
            break;
        }

        // 確率 1/2 で, 一点更新か二点 Swap かのいずれかを選ぶ
        // 一点更新 : ある日のコンテストのタイプを変更
        if (rand() % 2 == 0)
        {
            // 変更する日数とコンテストのタイプをランダムに選ぶ
            int d = rand() % input.length;
            int type = rand() % 26;

            // 元のコンテストのタイプを保存
            int old_type = contest_types[d];

            // d 日目のコンテストのタイプを変更
            contest_types[d] = type;

            // 満足度を計算
            ll new_score = evaluate(input, contest_types);

            if (current_score < new_score)
            {
                cerr << "iteration: " << iteration << ", new score = " << new_score << endl;
                current_score = new_score;
            }
            // 満足度が改善されなかった場合, 元に戻す
            else
            {
                contest_types[d] = old_type;
            }
        }
        // 二点 Swap: 比較的近い2つの日のコンテストのタイプを入れ替える
        else
        {
            // 入れ替える日をランダムに選ぶ
            int d1 = rand() % input.length;
            int d2 = d1 + (rand() % (min(input.length - 1, d1 + 14) - d1 + 1));

            // 実際に Swap
            swap(contest_types[d1], contest_types[d2]);

            ll new_score = evaluate(input, contest_types);
            if (current_score < new_score)
            {
                cerr << "iteration: " << iteration << ", new score = " << new_score << endl;
                current_score = new_score;
            }
            // 満足度が改善されなかった場合, 元に戻す
            else
            {
                swap(contest_types[d1], contest_types[d2]);
            }
        }

        iteration++;
    }

    cerr << "total iterations: " << iteration << endl;

    return Output{contest_types};
}

これを提出すると、 117,071,989 点(延長戦 76 位)を取ることができた。seed = 0の全体 Score は 2,050,521 である。

atcoder.jp

seed = 0の標準エラー出力はこんな感じ。

iteration: 0, new score = 255614
iteration: 19, new score = 261544
iteration: 37, new score = 273147
...
iteration: 86837, new score = 1148424
iteration: 116130, new score = 1148738
iteration: 134033, new score = 1149010
total iterations: 326919

結局13万回程度で更新は止まってしまっているが、それでも得点は結構伸びている。

焼きなまし法

さて、ここまでは「解を少しずつ変化させ、スコアが改善したら採用する」という方針でやってきた。

しかし、目的は最終的なスコアを最大化することであって、途中段階のスコアを改善することではない。

つまり、スコアが良いものだけを採用していては、解空間の局所最適解に嵌まって抜け出せなくなる可能性があるということだ。

したがって、多少スコアが悪くなっても確率でその変化を採用し、局所最適解から脱出するような方針を取りたい。

そのための便利なやり方が、「焼きなまし法」である。焼きなまし法についての理論的な話は、公式解説か AHF の資料、あるいは以下の記事を参照するとよい。

qiita.com


実装では、一点更新と二点 Swap それぞれで、変化させた解を採用するif文に確率の要素を加えることになる。

この確率 p(\Delta, T) は、

p(\Delta, T) = \left\{ \begin{array}{ll} 1 & (\Delta \geq 0) \\ e^{\tfrac{\Delta}{T}} & (\Delta \lt 0) \end{array} \right.

と書ける。ただし、 \Delta は解を変化させたことによるスコアの変化量、 T は温度(焼きなまし法の探索進度パラメータ)である。

この式は、スコアの変化が良い方向に向かったときは必ず採用し、悪い方向に向かったとしても探索進度が早いときほどより大きな確率で採用することを言っている。

温度 T は、焼きなまし法開始からの時間を区間 [0, 1] の範囲で正規化した値 t を用いて、

  • t \leftarrow T_0^{1-t} \times T_1^t

と指数的に更新する。ここで、 T_0 は開始時の温度、 T_1 は終了時の温度である。

焼きなまし法では、この2つの温度をうまい具合に決めるのが重要である。今回は T_0 = 2000, \, T_1 = 600 とした。


以上をC++で実装すると以下のようになる。

solve関数のコード(クリックで展開)

Output solve(const Input &input)
{
    /// @brief contest_types[d] := d 日目に開催されるコンテストのタイプ
    /// @note 貪欲法で求めた解を初期解とする
    vector<int> contest_types = greedy_solve(input).contest_types;
    /// @brief 現時点での満足度
    ll current_score = evaluate(input, contest_types);

    /// @brief 乱数生成器
    mt19937 rand{42};
    /// @brief 0.0 以上 1.0 未満の乱数を生成する
    uniform_real_distribution<double> rand_01(0.0, 1.0);

    /// @brief 焼きなまし法の開始時刻
    auto start_time = chrono::system_clock::now();
    /// @brief 制限時間(ミリ秒単位)
    const int LIMIT = 1900;

    /// @brief 焼きなましの開始温度
    const double START_TEMP = 2e3;
    /// @brief 焼きなましの終了温度
    const double END_TEMP = 6e2;
    /// @brief 現在の温度
    double current_temperature = START_TEMP;

    /// @brief 試行回数
    int iteration = 0;

    while (true)
    {
        auto current_time = chrono::system_clock::now();

        // 制限時間を超えたら終了
        if (chrono::duration_cast<chrono::milliseconds>(current_time - start_time).count() >= LIMIT)
        {
            break;
        }

        // 確率 1/2 で, 一点更新か二点 Swap かのいずれかを選ぶ
        // 一点更新 : ある日のコンテストのタイプを変更
        if (rand() % 2 == 0)
        {
            // 変更する日数とコンテストのタイプをランダムに選ぶ
            int d = rand() % input.length;
            int type = rand() % 26;

            // 元のコンテストのタイプを保存
            int old_type = contest_types[d];

            // d 日目のコンテストのタイプを変更
            contest_types[d] = type;

            // 満足度を計算
            ll new_score = evaluate(input, contest_types);

            if (current_score < new_score || rand_01(rand) < exp((new_score - current_score) / current_temperature))
            {
                cerr << "iteration: " << iteration << ", new score = " << new_score << endl;
                current_score = new_score;
            }
            // 満足度が改善されなかった場合, 元に戻す
            else
            {
                contest_types[d] = old_type;
            }
        }
        // 二点 Swap: 比較的近い2つの日のコンテストのタイプを入れ替える
        else
        {
            // 入れ替える日をランダムに選ぶ
            int d1 = rand() % input.length;
            int d2 = d1 + (rand() % (min(input.length - 1, d1 + 14) - d1 + 1));

            // 実際に Swap
            swap(contest_types[d1], contest_types[d2]);

            ll new_score = evaluate(input, contest_types);
            if (current_score < new_score || rand_01(rand) < exp((new_score - current_score) / current_temperature))
            {
                cerr << "iteration: " << iteration << ", new score = " << new_score << endl;
                current_score = new_score;
            }
            // 満足度が改善されなかった場合, 元に戻す
            else
            {
                swap(contest_types[d1], contest_types[d2]);
            }
        }

        iteration++;

        double progress = (double)chrono::duration_cast<chrono::milliseconds>(current_time - start_time).count() / (double)LIMIT;
        current_temperature = pow(START_TEMP, 1.0 - progress) * pow(END_TEMP, progress);
    }

    cerr << "total iterations: " << iteration << endl;

    return Output{contest_types};
}

これを提出すると、 120,603,130 点(延長戦 34 位)を取ることができた！延長戦順位表では上位約2%の得点である。

ちなみに、seed = 0の全体 Score は 1,978,491 であった。あれ、さっきよりも下がってる...？

atcoder.jp

実は、今まで全体 Score が 0 であった seed = 11が 97,729 となっている。このような局所最適解に嵌まりやすいケースで大きな改善が見られたと言えるだろう。


とりあえず延長戦順位表の2ページ目に入るところまでは行けたので、今回はここまでにしようと思う。

コメント

当初は上位10%を目指してやっていたのだが、ここまで伸びるとは思わなかった。焼きなまし法恐るべし。

今後も AHC の過去問をこんな感じでやっていきたい。