初めての～焼きなまし法～　抽象的なコードと具体的な実装までのメモ

この記事について
焼きなまし法って知ってますか???
1. とある確率ってなんやねん
焼きなまし法を理解するにあたって参考にした記事
～”抽象的な焼きなまし法”の完成形のコードはこちら～
~抽象～　細かく説明していきます
1. 抽象状態クラス
2. 抽象焼きなましクラス
~実際にコンテストで使ってみる～　具体例編
さいごに

この記事について

この記事は焼きなまし法について、筆者が学んだことをまとめた記事です
具体的には、以下の事を書きました

抽象クラス(状態・焼きなまし法)のコード
↑を用いた具体的な使用例のコード

焼きなまし法って知ってますか???

皆さん、焼きなまし法って知ってますか？(私は最近まで聞いたことがある。程度の認識でした)
知らない人のために簡単に説明すると

初期状態を作る
現在の状態から少し変更を加えた状態を作る
現在の状態と新しい状態を比較して、
新しい状態の方がスコアが良い、または状態が悪化しても”とある確率”で状態を新しい状態で更新する
2.に戻る

を繰り返して、より良い状態(解)を求める方法のことです

とある確率ってなんやねん

上の説明で「とある確率」と書いたものは、温度と呼ばれます。
始めてすぐのうちはポカポカ🌡️です。徐々に冷めていきます🧊
ポカポカの時のほうが、状態の大きい悪化も受け入れやすいです

どうして悪化を受け入れるの？？🤔
良い質問です。👍
悪化を受け入れないと、局所最適解に陥りやすいのです。

はい。ザックリとした焼きなまし法の説明終わりです。
具体的なコードで見ないと分かりづらいと思いますので、
これより、抽象的な焼きなまし法のコードと具体的な焼きなまし法のコードを書いていきます
(使用言語はC++です)

焼きなまし法を理解するにあたって参考にした記事

以下に書いたクラスのコードは↑こちらの記事のコードを基本として、再利用しやすいようにクラスにまとめたものです。
(注意: 基本的な処理の流れはだいたい同じですが、少し自分好みに作り替えています…)

こちらの記事を一度読んだうえで当記事を読むと、よりわかりやすいかも知れません

～”抽象的な焼きなまし法”の完成形のコードはこちら～

#include <bits/stdc++.h>

/// @brief 焼きなまし法の対象となる状態
template <typename score_type>
class state {
public:
    /// @brief 近傍の状態を取得
    virtual void get_neighbor_state() = 0;

    /// @brief スコア計算 
    virtual score_type score_caluculation() = 0;

    /// @brief 近傍の状態へ遷移
    virtual void transition_to_neighbor_state() = 0;
};

/// @brief 焼きなまし法クラス
template <typename state_type, typename score_type>
class simulated_annealing {
public:
    /// @brief コンストラクタ
    /// @param time_limit 焼きなます制限時間(秒)
    /// @param start_temp 初期温度
    /// @param end_temp 終了温度
    /// @param initial_state 初期状態(解)
    simulated_annealing(double time_limit, double start_temp, double end_temp, state_type initial_state)
        : time_limit(time_limit), start_temp(start_temp), end_temp(end_temp), current_temp(start_temp),
        current_state(initial_state), engine(seed_gen()), distribution(random_min_number, random_max_number) {
    }

    /// @brief 焼きなまし法の実行
    virtual state_type execute_simulated_annealing() {
        start_time = std::chrono::high_resolution_clock::now();
        current_temp = start_temp;
        while (true) {
            /*
            経過時間のチェック
            */
            current_time = std::chrono::high_resolution_clock::now();
            elapsed_time = std::chrono::duration_cast<std::chrono::duration<double> >(current_time - start_time).count();
            if (elapsed_time > time_limit) break;
            
            /*
            遷移部分
            */
            state_type new_state = current_state;
            new_state.get_neighbor_state();

            score_type current_score = current_state.score_caluculation();
            score_type new_score = new_state.score_caluculation();

            current_temp = temperature_calculation();

            double transition_probability = transition_probability_calculation(current_score, new_score, current_temp);

            /*
            遷移確率が遷移条件を満たせば更新
            */
            if (is_transition_allowed(transition_probability)) {
                current_state = new_state;
                current_state.transition_to_neighbor_state();
            }
        }
        return current_state;
    }


protected:
    /*
    時間管理
    */
    std::chrono::high_resolution_clock::time_point start_time;
    std::chrono::high_resolution_clock::time_point current_time;
    double elapsed_time;
    double time_limit;

    /*
    温度
    */
    double start_temp;
    double end_temp;
    double current_temp;

    /*
    状態
    */
    state_type current_state;

    /*
    乱数生成
    */
    std::random_device seed_gen;
    std::mt19937 engine;
    const double random_min_number = 0.0;
    const double random_max_number = 1.0;
    std::uniform_real_distribution<> distribution;


    /*
    各計算式メソッド
    */
    /// @brief 温度計算
    virtual double temperature_calculation() = 0;

    /// @brief 遷移確率計算
    virtual double transition_probability_calculation(score_type current_score, score_type new_score, double current_temp) = 0;

    /// @brief 遷移のための条件式 
    virtual bool is_transition_allowed(const double transition_probability) = 0;
};

~抽象～　細かく説明していきます

両方ともテンプレートを使った抽象クラスで作られているので、実際に使うときにはこれらを継承して実装することになります。

抽象状態クラス

状態は以下3つのことができればよいです

現在の状態から少し変更した状態(近傍)を取得(作る)
現在の状態のスコア計算
状態を実際に更新する

まだ抽象クラスなのでこれらは純粋仮想メソッドとして定義するだけです
使うときは、問題に応じてこれらをoverrideして使います

抽象焼きなましクラス

メンバ

変数

焼きなまし法に必要なメンバ変数を用意します

時間管理のための変数(単位は秒)
温度管理のための変数
保持する状態変数
乱数生成
これは基本的に遷移の条件式メソッド内で使う予定があるので予め用意しています。

メソッド

現在の温度計算
{スコア・現在温度}を使った遷移確率の計算
遷移の条件式

これらも純粋仮想メソッドとして用意します

コンストラクタ

変数	値
time_limit	制限時間秒で指定します(○.○秒とかもOK)
start_temp	参考にさせていただいた記事によると初期温度 = “一回の遷移で動きうるスコア幅の最大値程度”
end_temp	終了温度 = “一回の遷移で動きうるスコア幅の最小値程度” が目安となるそうです。
initial_state	初期状態です。てきと~に作ったものを入れてあげてくださいランダムで生成してもよいですし貪欲で決めてもよいです (「ここら辺は問題に応じて」考えてくださいね♬)

焼きなまし法のメインとなる処理をするメソッド: execute_simulated_annealing()

このメソッドは仮想メソッド(純粋でない)なので、ねじ込みたい処理があればoverrideしてねじ込んでもOKです

コードの内容を日本語でかいてみた…

経過時間のチェック

経過時間 = 現在の時刻 – 開始時刻
より、経過時間を取得します
もし制限時間を超えていればループを打ち切ります

遷移部分

まずは新しい状態(new_state)を現在の状態をコピーして作ります
次に、new_stateで近傍状態を得ます

{現在の状態, 近傍の状態}の二つで、それぞれのスコアを計算します

現在の温度を計算します
(また後ほど、計算方法を書きます)

遷移確率を計算します
(こちらも後ほど…)

更新部分

遷移確率を受け取った遷移の条件式メソッド の判定結果で状態を更新するか判断します
更新する場合(戻り値:true)は近傍状態を現在の状態として実際に更新します

~実際にコンテストで使ってみる～　具体例編

今回は練習問題として、こちらの問題をやってみました(参考記事で紹介されていたので…)

A - Just Write Numbers!

AtCoder is a programming contest site for anyone from beginners to experts. We hold weekly programming contests online.

提出コード(長いので一部(上で書いた抽象クラスの部分)は省略)

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

/*
省略した部分
抽象クラス(状態・焼きなまし)の定義
がここに書かれています
*/


/*
これ以下で具象を定義する!!
*/

/*
chokudai contest 004用の"状態"
*/
template <typename score_type>
class state_for_chokudai_contest_004 : public state<score_type> {
public:
    state_for_chokudai_contest_004(int n, int b_1, int b_2, int b_3)
    : n(n), b_1(b_1), b_2(b_2), b_3(b_3), score(-1), change_row(-1), change_column(-1), changed_number(-1) {
        table.assign(n, vector<int>(n, -1));
        table_min.assign(n, vector<int>(n, -1));
        table_max.assign(n, vector<int>(n, -1));
    }

    /*
    表の初期化と入出力メソッド
    */
    /// @brief 表[i][j]に入れられる値の範囲を指定する
    /// @param is_min true 最小値か  false 最大値
    void set_table_i_j_value_range(int row, int column, int number, bool is_min) {
        if (is_min) table_min[row][column] = number;
        else table_max[row][column] = number;
    }

    /// @brief 初期状態のテーブルを構築する
    void set_init_table() {
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                //まずは各テーブルの値がmin <= maxの範囲内のランダムな値になるようにする
                table[i][j] = table_min[i][j] + rand() % (table_max[i][j] - table_min[i][j] + 1);
            }
        }

        //構築と同時にスコアを求める
        /*和がb_iとなる区間の個数*/
        int b_1_count = 0;
        int b_2_count = 0;
        int b_3_count = 0;
        /*横向きの区間*/
        //行
        for (int i = 0; i < n; i++) {
            //列(区間の左端)
            for (int j = 0; j < n; j++) {
                int now_sum = 0;
                //列(区間の右端)
                for (int k = j; k < n; k++) {
                    now_sum += table[i][k];
                    if (now_sum == b_1) b_1_count++;
                    if (now_sum == b_2) b_2_count++;
                    if (now_sum == b_3) b_3_count++;
                }
            }
        }
        /*縦向きの区間*/
        //列
        for (int i = 0; i < n; i++) {
            //行(区間の上端)
            for (int j = 0; j < n; j++) {
                int now_sum = 0;
                //行(区間の下端)
                for (int k = j; k < n; k++) {
                    now_sum += table[k][i];
                    if (now_sum == b_1) b_1_count++;
                    if (now_sum == b_2) b_2_count++;
                    if (now_sum == b_3) b_3_count++;
                }
            }
        }

        score = b_1 * b_1_count + b_2 * b_2_count + b_3 * b_3_count;
    }

    /// @brief 表の出力
    void output_table() {
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                cout << table[i][j] << " ";
            }
            cout << endl;
        }
    }

    /*
    焼きなまし法を適用するためのメソッド
    */
    void get_neighbor_state() override {
        change_row = rand() % n;
        change_column = rand() % n;
        changed_number = table_min[change_row][change_column] + rand() % (table_max[change_row][change_column] - table_min[change_row][change_column] + 1);
    }

    score_type score_caluculation() override {
        //遷移する方でない場合
        if (changed_number == -1) return score;

        //"仮の変更"前の値を記録
        int before_change_number = table[change_row][change_column];

        /*和がb_iとなる区間の個数*/
        int before_b_1_count = 0;
        int before_b_2_count = 0;
        int before_b_3_count = 0;
        //横向き区間
        for (int i = 0; i < n; i++) {
            int now_sum = 0;
            for (int j = i; j < n; j++) {
                now_sum += table[change_row][j];
                if (now_sum == b_1) before_b_1_count++;
                if (now_sum == b_2) before_b_2_count++;
                if (now_sum == b_3) before_b_3_count++;
            }
        }
        //縦向き区間
        for (int i = 0; i < n; i++) {
            int now_sum = 0;
            for (int j = i; j < n; j++) {
                now_sum += table[j][change_column];
                if (now_sum == b_1) before_b_1_count++;
                if (now_sum == b_2) before_b_2_count++;
                if (now_sum == b_3) before_b_3_count++;
            }
        }
        //変更を行う前に、関係する部分のスコアを排除する
        score -= b_1 * before_b_1_count + b_2 * before_b_2_count + b_3 * before_b_3_count;

        //テーブルを一旦変更する
        table[change_row][change_column] = changed_number;

        /*和がb_iとなる区間の個数*/
        int after_b_1_count = 0;
        int after_b_2_count = 0;
        int after_b_3_count = 0;
        //横向き区間
        for (int i = 0; i < n; i++) {
            int now_sum = 0;
            for (int j = i; j < n; j++) {
                now_sum += table[change_row][j];
                if (now_sum == b_1) after_b_1_count++;
                if (now_sum == b_2) after_b_2_count++;
                if (now_sum == b_3) after_b_3_count++;
            }
        }
        //縦向き区間
        for (int i = 0; i < n; i++) {
            int now_sum = 0;
            for (int j = i; j < n; j++) {
                now_sum += table[j][change_column];
                if (now_sum == b_1) after_b_1_count++;
                if (now_sum == b_2) after_b_2_count++;
                if (now_sum == b_3) after_b_3_count++;
            }
        }
        //変更後のスコアを求める
        score += b_1 * after_b_1_count + b_2 * after_b_2_count + b_3 * after_b_3_count;
        //テーブルを元に戻す
        table[change_row][change_column] = before_change_number;
        return score;
    }

    void transition_to_neighbor_state() override {
        table[change_row][change_column] = changed_number;
        change_row = -1;
        change_column = -1;
        changed_number = -1;
    }

private:
    /*
    表関連
    */
    int n;
    static vector<vector<int> > table;
    /*
    table[i][j]に代入できる値の 最小・最大値
    */
    static vector<vector<int> > table_min;
    static vector<vector<int> > table_max;

    /*
    スコア計算のための値
    */
    int b_1;
    int b_2;
    int b_3;
    score_type score;

    /*
    表に加わる変更
    */
    int change_row;
    int change_column;
    int changed_number;
};

//外部でも定義する
template <typename score_type>
vector<vector<int> > state_for_chokudai_contest_004<score_type>::table;
template <typename score_type>
vector<vector<int> > state_for_chokudai_contest_004<score_type>::table_min;
template <typename score_type>
vector<vector<int> > state_for_chokudai_contest_004<score_type>::table_max;


/*
chokudai contest 004用の焼きなまし法
*/
template <typename state_type, typename score_type>
class simulated_annealing_for_chokudai_contest_004 : public simulated_annealing<state_type, score_type> {
public:
    simulated_annealing_for_chokudai_contest_004(double time_limit, double start_temp, double end_temp, state_type init_state)
    : simulated_annealing<state_type, score_type>(time_limit, start_temp, end_temp, init_state) {}


protected:
    double temperature_calculation() override {
        return this->start_temp + (this->end_temp - this->start_temp) * this->elapsed_time / this->time_limit;
    }

    double transition_probability_calculation(int current_score, int new_score, double current_temp) override {
        /*テンプレ最大化 今回は最大化する問題だからこっち!!*/
        return exp((new_score - current_score) / current_temp);

        /*
        テンプレ最小化
        return exp((current_score - new_score) / current_temp);
        */
    }

    bool is_transition_allowed(const double transition_probability) override {
        return transition_probability > this->distribution(this->seed_gen);
    }
};


int main() {
    int n, b_1, b_2, b_3;
    cin >> n >> b_1 >> b_2 >> b_3;
    
    //初期状態
    state_for_chokudai_contest_004<int> init_state(n, b_1, b_2, b_3);
    //各最小値の設定
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            int min_num;
            cin >> min_num;
            init_state.set_table_i_j_value_range(i, j, min_num, true);
        }
    }
    //各最大値の設定
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            int max_num;
            cin >> max_num;
            init_state.set_table_i_j_value_range(i, j, max_num, false);
        }
    }
    //初期テーブル(表)の構築
    init_state.set_init_table();

    //焼きなまし法をインスタンス化
    simulated_annealing_for_chokudai_contest_004<state_for_chokudai_contest_004<int>, int> sa(2.9, 50, 0, init_state);

    state_for_chokudai_contest_004<int> answer = sa.execute_simulated_annealing();

    //出力
    answer.output_table();
}

具象状態クラス

具象状態クラスを宣言する上での注意点

問題にもよりますが、状態を表すのに配列などを用いる場合(盤面のような状態を持つ場合等)は状態のコピーが重い処理となるので

「状態クラスに静的変数として宣言する」
「同じ情報を参照するようにする」

などをした上で、差分のみを持って近傍を考えるようにすると余計なコピーをせずに済むので良きかもです!

差分で考えること自体は正解だと思うのですが、それを実現するための方法として上記2つの方法が良いのかは分かりません
(筆者の経験不足です…強強erたちはどのようにしてコピーをなくしているのか気になる…)

メンバ変数　(とくに重要なもの)

メンバ変数	役割
table	表、実際に各マスに入れる値を記録する静的メンバなので、インスタンスが違くても共有される
table_min table_max	各マスに入れられる最小値、最大値を記録する同じく静的メンバ。
change_row	変更する行番号を記録する
change_column	変更する列番号を記録する
changed_number	表[change_row][change_column]に入れる値
change_~~系についてのまとめ	change_~~で近傍状態を表すことで、表のコピーをなくすことができる。要は差分だけを状態として持つわけである…

overrideしたメソッドの説明

メソッド名	実装
get_neighbor_state	表の範囲内で、ランダムな位置を選び、change_{row, column}に格納するその位置に入れられる値 {min_{r, c} <= changed_number <= max_{r, c}} をランダムに選んで格納する
score_caluculation	1.「状態を変更した場合」についてのスコア計算をする必要がある 2.毎回全ての範囲について計算するのは、重たい処理すぎるので… ↓ 変更される位置と値 {change_row, change_column, changed_number} を使って関係する行と列のみ再計算してスコアを高速に求める。 table_変更位置の値を仮変更し、スコア計算後に元の値に戻すことで求めている。
transition_to_neighbor_state	共通している表(table)に「行った変更」を実際に適用する。「これを呼びだした状態」が次の「現在の状態」になるわけなので、change_~~系の変数に無効な値を入れなおす。

具象焼きなましクラス

焼きなまし法に関するメソッドをoverrideする

ここで計算式にしているもの*1*2は、先ほどの参考記事を見て書いております
(あの記事には、まじで感謝です🥰🥰🥰)

メソッド名	実装
temperature_calculation	“線形でstart_tempからend_tempに減少する”*1
transition_probability_calculation	スコア最大化2 “ exp((new_score – current_score) / current_temp) “ スコア最小化* exp((current_score – new_score) / current_temp) こちらもこれで良さそう今回は最大化する方を使う
is_transition_allowed	遷移確率 > 0.0から1.0の乱数で遷移する (遷移確率計算式よりスコア改善時は遷移確率 >= 1となる)