平面 traveling salesman problem の厳密解を求める．

平面tspの分岐限定法のコードをなんとなく書いてみたけれども，全然早くならなかった．

Traveling Salesman Problem (TSP)

• N個の都市と都市間の距離が与えられる．全ての都市をちょうど一度ずつ巡って始点に戻る順回路の中で移動距離が最小のものを出力する．

平面TSP

• 都市がユークリッド座標上に存在して，自由に頂点間を行き来出来る．他はTSPと同じ．

実装

要件

• 入力: 都市の座標の集合．vector<Point>．
• 出力: 都市の辿り方．vector<int>．i番目の要素は，i番目にどの都市へ向かえば良いか記述されている．

動的計画法

• 本記事ではbitDPはやらない．
• 蟻本を買いましょう．
• このままの実装だとbitDPの方が速いかも

分岐限定法

• 深さ優先探索(dfs・再帰)で次の頂点を選ぼうとしているとする．
• どの頂点を選んでも暫定解を更新できないならば，それ以上は探索しなくても良い．
• 『暫定解を更新できない』ことを確かめるにはどうしたら良いか？

暫定解が更新できなさそう

• 通ってきたパスの長さ+残りの頂点を巡回するパスの最短の長さが暫定解を超えたらその先は探索しない．
• 残りの頂点を巡回するパスの最短の長さはそもそも部分問題なので，厳密に計算する訳にはいかない．
  • 代わりに最小全域木を使う．最小全域木の辺の長さの総和は，巡回するパスの最短の長さ以下．
  • 最小全域木は比較的高速に求めることができる．
• つまり，通ってきたパスの長さ+残りの頂点から構成される最小全域木が暫定解を超えたらその先は探索しない．

実装

• 以下に実装を示す．
• 定数倍の高速化はあまり意識しないようにした．
  • 例えばpriority_queueを外に出せばちょっとは早くなるんじゃないかな(適当)
• 厳密には，下の実装は残りの頂点を巡回するパスの最短の長さ = 残りの頂点から構成される最小全域木の長さではない．

// union-find
class unionfind {
public:
    vector<int> data;
    unionfind(int size) : data(size, -1) { }
    bool connect(int x, int y) {
        x = root(x); y = root(y);
        if (x != y) {
            if (data[y] < data[x]) swap(x, y);
            data[x] += data[y]; data[y] = x;
        }
        return x != y;
    }
    inline bool same(int x, int y) {
        return root(x) == root(y);
    }
    inline int root(int x) {
        return data[x] < 0 ? x : data[x] = root(data[x]);
    }
    inline int size(int x) {
        return -data[root(x)];
    }
};


typedef double numeric;
// 座標
struct Point {
    numeric x;
    numeric y;
    Point(double _x = 0, double _y = 0) :x(_x), y(_y) { }
};

// 平面TSP
vector<int> tsproblem(const vector<Point>& points) {
    // 頂点数
    const int n = static_cast<int>(points.size());
    const numeric inf = 1e100;

    // 暫定解の距離
    numeric best = inf;
    // 暫定解(idx番目の頂点はval番目に訪れる．val=0は未だ訪れていない)
    vector<int> best_perm(n);

    // 頂点u,v間の距離
    function<numeric(int, int)> calcdist = [&points](int u, int v) {
        return hypot(points[u].x - points[v].x, points[u].y - points[v].y);
    };

    // dist[i][j] iとjの距離
    vector<vector<numeric>> dist(n);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            dist[i].push_back(calcdist(i, j));

    // distrank[i][j] : iにj番目に近い頂点番号
    vector<vector<int>> distrank(n);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; j += ++j == i)
            distrank[i].push_back(j);
        sort(distrank[i].begin(), distrank[i].end(), [&dist, &i](int u, int v) {return dist[u][i] < dist[v][i]; });
    }

    // visited[v]=0な頂点について，最小全域木の長さを求める．クラスカル法
    function<numeric(const vector<int>&)> spantreedist = [&](const vector<int>& visited) {
        unionfind uf(n);
        int nc = 0;
        vector<pair<numeric, pair<int, int>>> pq;

        for (int i = 0; i < n; ++i) {
            if (visited[i]) continue;
            for (int j = i + 1; j < n; ++j) {
                if (visited[j]) continue;
                pq.emplace_back(dist[i][j], make_pair(i, j));
            }
            ++nc;
        }
        sort(pq.begin(), pq.end());
        numeric ans = 0;
        for (int i = 0; 1 < nc; ++i) {
            auto p = pq[i].second;// pq.pop();
            if (!uf.connect(p.first, p.second)) continue;
            --nc;
            ans += dist[p.first][p.second];
        }
        return ans;
    };

    function<void(vector<int>, int, int, numeric)> dfs = [&](vector<int> visited, int pos, int cnt, numeric len) {
        // 最後の頂点に到達
        if (cnt == n) {
            // 暫定解と比較
            len += dist[pos][0];
            if (len < best) {
                best_perm = visited;
                best = len;
            }
            return;
        }

        // 残りの頂点を巡回するパスの最短の長さ(近似)
        numeric leftdistmin = spantreedist(visited);
        {
            numeric l = inf;
            for (int i = 1; i < n; ++i)
                if (!visited[i])
                    l = min(l, dist[0][i]);
            leftdistmin += l;
        }

        ++cnt;
        // 近い頂点から順にチェック
        for (int i : distrank[pos]) {
            if (visited[i]) continue;
            visited[i] = cnt;
            numeric l = len + dist[pos][i];
            if (len + leftdistmin < best)
                dfs(visited, i, cnt, l); // 暫定解を更新できる可能性がある場合のみ，進める
            visited[i] = 0;
        }
    };

    best_perm[0] = 1;
    dfs(best_perm, 0, 1, 0.0);

    // 変換
    vector<int> tourist(n);
    for (int i = 0; i < n; ++i)
        tourist[best_perm[i] - 1] = i;
    return tourist;
}

検証用

• 出力は https://csacademy.com/app/geometry_widget/ の形になってます．

mt19937_64 randdev(8901016);
inline double rand_range(double l, double h) {
    return uniform_real_distribution<double>(l, h)(randdev);
}

int main() {

    {
        const int n = 19;
        vector<Point> points;
        for(int i = 0; i < n; ++i)
            points.emplace_back(rand_range(-10.0, 10.0), rand_range(-10.0, 10.0));

        for (auto p : points)
            printf("%.4f %.4f\n", p.x, p.y);

        auto ans = tsproblem(points);

        int p = ans.back();
        for (auto to : ans) {
            printf("Segment %.4f %.4f %.4f %.4f\n",
                points[p].x, points[p].y, points[to].x, points[to].y);
            p = to;
        }

        printf("Line 0 1 0\nLine 1 0 0\n");
    }


    return 0;
}