焼きなまし法によるTSP (paizaランク S 相当)

巡回セールスマン問題とは、都市の集合と各都市間の距離が与えられ、全都市をちょうど1回ずつ訪れたのち出発した都市に戻ってくるような経路 (巡回路) のうち最も短いものを求める問題です。

ここでは、焼きなまし法と呼ばれるアルゴリズムを学習しましょう。このアルゴリズムは発見的解法 (ヒューリスティクス) と呼ばれるもので、ある程度良い解が出力されることが期待されるものの、解の精度は全く保証されません。

一般的な焼きなまし法の概要は、以下の通りです。

・ 解を適当に作り、暫定解 X とする
・ 仮想的な温度 T を設定する
・ 温度 T を受け取り 0 以上 1 以下の実数を返す関数 P を用意する。P は T が大きければ 1 に近い値を、小さければ 0 に近い値を返すように設計する
・ 一定の回数か、解があまり改善されなくなるか、温度 T がある程度下がるまで以下を繰り返す
    ・ 暫定解 X を少しだけ変更したもの (X の近傍) を 1 つ選択して X' とする。X' が X より真に良い場合は確率 1 で、そうでない場合は確率 P(T) で X を X' で更新する。最後に温度 T を少し下げる
・ X を解として出力する

アルゴリズムの最初の方では温度 T が高いため、悪くなってしまう場合であっても解を更新する確率が高くなります。つまり、アルゴリズムの初期段階では解が大胆に変わることが多くなります。アルゴリズムの実行が進み温度 T が低くなっていくと、良くなる場合でないと滅多に更新されなくなっていきます。このようにすることで、解が局所的最適解に留まってしまうことが防がれ、大域的最適解へと近づいていくことが期待できます (期待ができるというだけで、必ず大域的最適解に辿り着く保証があるというわけではない)。

このヒューリスティクスを巡回セールスマン問題に適用しましょう。本問では、以下のようにします。

・ 巡回路を適当に作り、初期解とする
・ 仮想的な温度 T を設定する
・ 温度 T を受け取り 0 以上 1 以下の実数を返す関数 P を用意する。P は T が大きければ 1 に近い値を、小さければ 0 に近い値を返すように設計する
・ 温度 T がある程度下がるまで以下を繰り返す
    ・ (*) 巡回路を成す辺を 2 本選び、その辺を繋ぎ変えることを考える。繋ぎ変えることによって巡回路長が短くなるなら繋ぎ変える。短くならないなら、確率 P(T) で繋ぎ変える。最後に温度 T を少し下げる
・ 巡回路を出力する

辺の繋ぎ変え方は 2-opt 法と同様です。例えば、辺 (a, d) と (b, c) を繋ぎ変えると 辺 (a, b) と (c, d) になります。



n 個の都市 (都市 0、都市 1、...、都市 n-1) のデータと、n 個の都市を巡る初期巡回路が与えられます。焼きなまし法を実装し、巡回路を求めてください。ただし、ジャッジの都合上で解を一意にするために、以下の指示に従ってください。

・ 温度 T の初期値は 1000
・ 温度 T の更新は T ← T x 0.99 でおこなう
・ (*) の処理は T >= 0.1 の間繰り返す
・ 関数 P は exp((d_before-d_after)/T) とする (ただし、exp は底がネイピア数 e である指数関数を、d_before と d_after はそれぞれ辺を繋ぎ変える前と後の巡回路長を表す)

(※ 関数 P は d_before ≦ d_after のときに定義されるため、(d_before-d_after)/T は常に 0 以下となります。従って、関数 P は 0 以上 1 以下の値を返します)

(※ なお、これらのパラメータはジャッジの都合により探索回数が少なくなるように設定されています。温度 T の下げ方を緩やかにしたり、関数 P を変えるなどすると、最終的に得られる解が良くなったり悪くなったりします。興味のある人はパラメータを調整しながら実験をやってみてください)

辺を 2 本選ぶ際と確率を計算する際には、以下で指定する疑似乱数生成プロシージャ (xorshift32) を用いてください。(unsigned int は符号なし 32 bit 整数型)

unsigned int state = 813;

unsigned int xorshift32(){
    unsigned int x = state;
    x ^= x << 13;
    x ^= x >> 17;
    x ^= x << 5;
    state = x;
    return state;
}

void pick_two(int &a, int &b){
    while(true){
        a = xorshift32()%n;
        b = xorshift32()%n;
        if(a > b){
            swap(a, b);
        }
        if(!(a+1 < b) || (a == 0 && b == n-1)){
            continue;
        }
        return;
    }
}

double rnd(){
    return (xorshift32()-1)/4294967294;
}

焼きなまし法の疑似コードを以下に示します。なお、SWAP(a, d, b, c)は巡回路 tour の辺 (a, d) と (b, c) を繋ぎ変えて辺 (a, b) と (c, d) にする処理を表します。

void simulated_annealing(int tour[]){
    double t = 1.0;
    while(t >= 0.1){
        int a, b;
        pick_two(a, b);
        
        d_before = (現在の巡回路 tour の巡回路長)
        d_after = (現在の巡回路 tour の都市 tour[a] と tour[(a+1)%n] を結ぶ辺と、都市 tour[b] と tour[(b+1)%n] を結ぶ辺を繋ぎ変えた巡回路の巡回路長)

        if(d_after < d_before || rnd() <= exp((d_before-d_after)/t)){
            // 辺を繋ぎ変え巡回路を更新する
            SWAP(tour[a], tour[a+1], tour[b], tour[(b + 1)%n]);
        } 

        // 温度を更新する
        t *= 0.99;
    }
}

なお、各都市は二次元平面上の点として与えられ、都市間の距離にはユークリッド距離を用いるものとします。
巡回路の始点は任意です。巡回路が正しければ正解と判定されます。