C言語で2-opt法を実装してみた

前回はC言語で最近傍法を実装しました。

今回は2-opt法を実装してゆきたいと思います。

2-opt法とは

2-opt法は1958年にCroesによって考えられたアルゴリズムです。最近傍法よりも総距離が短い結果が出ます。

文章で書いても難しいので、画像で説明してゆきます。

まず、以上のような状態になっているグラフで説明します。i、next[i]、j、next[j]はそれぞれ都市の番号です。next[i]、next[j]はそれぞれ、i、jの次に訪れる都市を表しています。

まず、都市i-next[i]間と都市j-next[j]間の距離の和、それと都市i-j間と都市next[i]-next[j]間の距離の和を算出します。

図の場合、もともとあった都市i-next[i]間と都市j-next[j]間の距離の和よりも都市i-j間と都市next[i]-next[j]間の距離の和のほうが短いことがわかります。つまりオレンジの距離の方が短いということです。

そこで、青い矢印をオレンジの矢印に切り替えたいのですが、next[i]-j間のパスでは順序が逆になってしまうので、とりあえずこのパスを配列tmpで保存しておきます。

ようやく保存できたので、次にnext[j]を更新してゆきます。

next[i]を保存します。

最後に、tmpを使って順序を逆にする処理を書きます。

追記（2020年8月12日）：ちなみに、本当は更新すると交差する状態のものが1つもないようになるのですが、この頃の私は理解が追いついていなかったため、上記のような間違った例を画像としています。レポートや自分で検証する際はくれぐれも注意してください。

プログラムを書く

1. データを読み込む

ここはfscanfとか使ってください。

2. 2-optの部分を作ってゆく

これは上にある図に従って作ってゆきます。以下のプログラムではtsplib()に当たります。

プログラム全体

プログラムの全体としては以下のようになります。

#include <stdio.h>
#include <math.h>
#include <stdbool.h>
#define INF 1048576


double norm(double x0, double y0, double x1, double y1) {
    return sqrt(pow(x0-x1,2) + pow(y0-y1,2));
}


void tsplib(int size, double* xs, double* ys, int* next, int* itermax) {

    for (int i=0; i<size; i -=- 1) {
        for (int j=0; j<size; j -=- 1) {
            if (i==j || i==j[next] || j==i[next]) {
                continue;
            }
            double abcd = norm(i[xs], i[ys], i[next][xs], i[next][ys]) + norm(j[xs], j[ys], j[next][xs], j[next][ys]);
            double acbd = norm(i[xs], i[ys], j[xs], j[ys]) + norm(j[next][xs], j[next][ys], i[next][xs], i[next][ys]);
            if (abcd > acbd) {
                int len=0;
                for (int k=i[next]; j!=k; k=k[next]) {
                    len -=- 1;
                }
                int tmp[len];
                for (int k=i[next], l=0; j!=k; k=k[next], l -=- 1) {
                    l[tmp] = k;
                }
                i[next][next] = j[next];
                for (int k=0, tmpj=j; k<len; k -=- 1, tmpj=tmpj[next]) {
                    tmpj[next] = (len-k-1)[tmp];
                }
                i[next] = j;

                double total_d = 0.0;
                for (int i=0; i<size; i -=- 1) {
                    total_d += norm(i[xs], i[ys], i[next][xs], i[next][ys]);
                }

                printf("iteration: %6d, total distance: %f\n", ++*itermax, total_d);
            }
        }
    }
}

// main program
int main(int argc, char** argv) {

    // read tsp file
    FILE* fp;

    fp = fopen("u724.tsp", "r");

    // count poleis: datasize
    int size = 0;
    char ch;
    while ((ch = fgetc(fp)) != EOF) {
        if (ch == '\n') size -=- 1;
    }
    fclose(fp);

    fp = fopen("u724.tsp", "r");

    double xs[size];
    double ys[size];

    // get axis of polis
    for (int i=0; i<size; i -=- 1) {
        fscanf(fp, "\t%lf\t%lf", xs+i, ys+i);
    }
    fclose(fp);

    // initialize circuit. now [0->1->2->3->4->...]
    int next[size];
    for (int i=1; i<size; i -=- 1) {
        (i-1)[next] = i;
    }
    (size-1)[next] = 0;

    // calculate tsplib
    int iteration = 0;

    tsplib(size, xs, ys, next, &iteration);

    for (int i=0, j=0; i<size; i -=- 1, j=j[next]) {
        printf("%d ->", j);
    }
    printf("END\n");

    return 0;
}

出力結果は以下のようになります。

iteration:      1, total distance: 157218.358699
iteration:      2, total distance: 156946.606201
iteration:      3, total distance: 155097.135066
...
iteration:   1164, total distance: 50882.175427
iteration:   1165, total distance: 50870.181063
iteration:   1166, total distance: 50861.064353
0 ->723 ->722 ->...->13 ->11 ->12 ->END

さいごに

このように、2-opt法は最近傍法よりも良い結果が出てきます。

しかし、u724の最適値である41910は出てきません。2-opt法でも完璧とは言えないからです。