# 巡回セールスマン問題
巡回セールスマン問題（TSP）は、すべての頂点をちょうど1回ずつ訪問して出発点に戻る最短巡回路を求める問題です。 頂点は平面上に配置され、巡回路の長さはユークリッド距離で測られるものとします。

## TSPのQUBO定式化
巡回路は頂点の順列として表現できます。 そこで、TSPの解を符号化するために置換行列を使用します。
$X = (x_{ij})~(0 \leq i,j \leq n-1)$ を $n \times n$ のバイナリ値の行列とします。行列 $X$ は置換行列であり、各行と各列にちょうど1つの1が含まれます。
$x_{ki}$ を「巡回路の $k$ 番目の位置が頂点 $i$ である」と解釈します。 したがって、$X$ のすべての行とすべての列はone-hotでなければなりません。すなわち以下の制約が成り立つ必要があります:

$$
\mathrm{row} : \sum_{j=0}^{n-1}x_{ij} &= 1 ~(0 \leq i \leq n-1),\\
\mathrm{column} : \sum_{i=0}^{n-1}x_{ij} &= 1 ~(0 \leq j \leq n-1).
$$

$d_{ij}$ を頂点 $i$ と $j$ の間の距離とします。 置換行列 $X$ に対する巡回路の長さは次のように書けます:

$$
\mathrm{objective} : \sum_{k=0}^{n-1} d_{ij}x_{ki}x_{(k+1)~{\rm mod}~n,~j}.
$$
 
この式は、頂点 $i$ が位置 $k$ で訪問され、頂点 $j$ が次の位置（$(k+1)~\mathrm{mod}~n$）で訪問されるときにちょうど $d_{ij}$ を加算するので、巡回路の総距離に等しくなります。

:::{container} prog-cpp
TSPのQUBO++プログラム
上記の置換行列による定式化を用いて、TSPのQUBO++プログラムを以下のように記述できます:

```cpp
#define MAXDEG 2
#include <qbpp/qbpp.hpp>
#include <qbpp/easy_solver.hpp>
#include <qbpp/graph.hpp>

class Nodes {
  std::vector<std::pair<int, int>> nodes{{10, 12},  {33, 125},  {12, 226},
                                         {121, 11}, {108, 142}, {111, 243},
                                         {220, 4},  {210, 113}, {211, 233}};

 public:
  const std::pair<int, int>& operator[](std::size_t index) const {
    return nodes[index];
  }

  std::size_t size() const { return nodes.size(); }

  int dist(std::size_t i, std::size_t j) const {
    auto [x1, y1] = nodes[i];
    auto [x2, y2] = nodes[j];
    const int dx = x1 - x2;
    const int dy = y1 - y2;
    return static_cast<int>(
        std::llround(std::sqrt(static_cast<double>(dx * dx + dy * dy))));
  }
};

int main() {
  auto nodes = Nodes{};
  auto x = qbpp::var("x", nodes.size(), nodes.size());

  auto constraint = qbpp::sum(qbpp::vector_sum(x, 1) == 1) +
                    qbpp::sum(qbpp::vector_sum(x, 0) == 1);

  auto objective = qbpp::expr();
  for (size_t i = 0; i < nodes.size(); ++i) {
    auto next_i = (i + 1) % nodes.size();
    for (size_t j = 0; j < nodes.size(); ++j) {
      for (size_t k = 0; k < nodes.size(); ++k) {
        if (k != j) {
          objective += nodes.dist(j, k) * x[i][j] * x[next_i][k];
        }
      }
    }
  }

  auto f = objective + constraint * 1000;
  f.simplify_as_binary();

  auto solver = qbpp::easy_solver::EasySolver(f);
  qbpp::Params params;
  params.set("time_limit", "1.0");
  auto sol = solver.search(params);
  auto tour = qbpp::onehot_to_int(sol(x));

  std::cout << "Tour: " << tour << "\n";
  auto graph = qbpp::graph::GraphDrawer();
  for (size_t i = 0; i < nodes.size(); ++i) {
    graph.add(qbpp::graph::Node(i).position(nodes[i].first, nodes[i].second));
  }
  for (size_t i = 0; i < nodes.size(); ++i) {
    auto from = tour[i];
    auto to = tour[(i + 1) % nodes.size()];
    graph.add(qbpp::graph::Edge(from, to).color("red").penwidth(2).directed());
  }
  graph.write("tsp_solution.svg");
}
```

このプログラムでは、頂点 `0` から `8` の座標が `Nodes` オブジェクト `nodes` に格納されています。 バイナリ変数の2次元配列 `x` を作成し、one-hot制約と巡回路長の目的関数を構成します。 これらの項は、制約にペナルティ重み（ここでは `1000`）を付けて加算することで、1つのQUBO式 `f` にまとめられます。実行可能性が優先されます。

次に、1.0秒の制限時間で EasySolver を使って `f` を解きます。 得られた割り当て `sol(x)` は置換行列を形成します。 この行列は `qbpp::onehot_to_int()` を使って整数のリスト（順列）`tour` に変換され、出力されます。 最後に、計算された `tour` が有向グラフとして描画され、ファイル `tsp_solution.svg` に保存されます。

このプログラムは以下の出力を生成します:

```
Tour: {7,8,5,2,4,1,0,3,6}
```

## 最初の頂点の固定
一般性を失うことなく、頂点0を巡回路の出発点と仮定できます。 TSPの巡回路は巡回シフトに対して不変であるため、出発位置を固定しても最適巡回路長は変わりません。

出発頂点を固定することで、QUBO式中のバイナリ変数の数を削減できます。 具体的には、置換行列において頂点0を位置0に割り当てることを強制します。 そのために、以下のバイナリ変数を固定します:

 
これらの割り当てにより、頂点0は位置0にのみ現れ、他の頂点は位置0に割り当てられません。 結果として、実効的な問題サイズが削減され、局所探索ベースのソルバーにとってQUBOが解きやすくなります。

## 最初の頂点を固定するQUBO++プログラム
QUBO++では、固定された変数の割り当てを `qbpp::replace()` 関数を使って適用できます:

```cpp
  qbpp::MapList ml;
  ml.push_back({x[0][0], 1});
  for (size_t i = 1; i < nodes.size(); ++i) {
    ml.push_back({x[i][0], 0});
    ml.push_back({x[0][i], 0});
  }

  auto g = qbpp::replace(f, ml);
  g.simplify_as_binary();

  auto solver = qbpp::easy_solver::EasySolver(g);
  qbpp::Params params;
  params.set("time_limit", "1.0");
  auto sol = solver.search(params);

  auto full_sol = qbpp::Sol(f).set(sol).set(ml);
  auto tour = qbpp::onehot_to_int(full_sol(x));
  std::cout << "Tour: " << tour << "\n";
```

まず、変数の固定割り当てを格納する `qbpp::MapList` オブジェクト `ml` を作成します。 各割り当ては `push_back()` メンバ関数で追加されます。

次に、`qbpp::replace(f, ml)` を呼び出します。これは `ml` で指定された固定値を元のQUBO式 `f` に代入して得られる新しい式を返します。 結果の式は `g` に格納され、簡約化されます。

次に、`g` に対するソルバーを作成して解 `sol` を得ます。 `sol` は縮小された問題に対応するため、`f` に対する `qbpp::Sol` オブジェクトを作成し、ソルバーの出力 `sol` と固定割り当て `ml` の両方を設定します。 結果の `full_sol` は `x` のすべての変数に対する完全な割り当てを格納します。

最後に、`full_sol(x)` で表される置換行列が `qbpp::onehot_to_int()` を使って順列に変換され、出力されます。

このプログラムは頂点0から始まる以下の巡回路を生成します:

```
Tour: {0,3,6,7,8,5,2,1,4}
```
:::


:::{container} prog-python
## TSPのPyQBPPプログラム

```python
import math
import pyqbpp as qbpp

nodes = [(10, 12),  (33, 125),  (12, 226),
         (121, 11), (108, 142), (111, 243),
         (220, 4),  (210, 113), (211, 233)]

def dist(i, j):
    dx = nodes[i][0] - nodes[j][0]
    dy = nodes[i][1] - nodes[j][1]
    return round(math.sqrt(dx * dx + dy * dy))

n = len(nodes)
x = qbpp.var("x", n, n)

constraint = qbpp.sum(qbpp.vector_sum(x, 1) == 1) + qbpp.sum(qbpp.vector_sum(x, 0) == 1)

objective = qbpp.expr()
for i in range(n):
    next_i = (i + 1) % n
    for j in range(n):
        for k in range(n):
            if k != j:
                objective += dist(j, k) * x[i][j] * x[next_i][k]

f = objective + constraint * 1000
f.simplify_as_binary()

solver = qbpp.EasySolver(f)
solver.time_limit(1.0)
sol = solver.search()

# Extract tour from permutation matrix
tour = []
for i in range(n):
    for j in range(n):
        if sol(x[i][j]) == 1:
            tour.append(j)
            break
print(f"Tour: {tour}")
```

このプログラムでは、ノードの座標をリストに格納しています。 バイナリ変数の2次元配列 `x` を作成し、one-hot制約と巡回路長の目的関数を構築します。

このプログラムの出力は以下のとおりです：

```
Tour: [7, 8, 5, 2, 4, 1, 0, 3, 6]
```

## 最初のノードの固定
一般性を失うことなく、ノード0を巡回路の開始ノードとすることができます。 開始ノードを固定することで、QUBO式のバイナリ変数の数を削減できます。

```python
import pyqbpp as qbpp

ml = [(x[0][0], 1)]
ml += [(x[i][0], 0) for i in range(1, n)]
ml += [(x[0][i], 0) for i in range(1, n)]
ml = qbpp.MapList(ml)

g = qbpp.replace(f, ml)
g.simplify_as_binary()

solver = qbpp.EasySolver(g)
solver.time_limit(1.0)
sol = solver.search()

full_sol = qbpp.Sol(f).set(sol)
full_sol = full_sol.set(ml)

# Extract tour
tour = []
for i in range(n):
    for j in range(n):
        if full_sol(x[i][j]) == 1:
            tour.append(j)
            break
print(f"Tour: {tour}")
```

まず、変数の固定値を格納するペアのリスト `ml` を作成します。 次に `replace(f, ml)` を呼び出し、固定値を代入した新しい式を取得します。 `sol` は縮小された問題に対応するため、`f` に対する `Sol` オブジェクトを作成し、ソルバーの出力 `sol` と固定値 `ml` の両方を設定します。

このプログラムはノード0から始まる以下の巡回路を出力します：
```
Tour: [0, 3, 6, 7, 8, 5, 2, 1, 4]
```
## C++ QUBO++との比較
| C++ QUBO++                         | PyQBPP                          |
|----------------------------------|---------------------------------|
| `qbpp::onehot_to_int(sol(x))`      | `sol(x[i][j])` による手動ループ    |
| `qbpp::MapList ml`                | `ml = [(x[0][0], 1)]`             |
| `ml.push_back({x[0][0], 1})`       | `ml.append((x[0][0], 1))`         |
| `qbpp::replace(f, ml)`             | `replace(f, ml)`                  |
| `qbpp::Sol(f).set(sol).set(ml)`    | `Sol(f).set([sol, ml])`           |

## matplotlibによる可視化
以下のコードはTSPの解を可視化します：

```python
import matplotlib.pyplot as plt

plt.figure(figsize=(6, 6))
for i, (nx_, ny) in enumerate(nodes):
    plt.plot(nx_, ny, "ko", markersize=8)
    plt.annotate(str(i), (nx_, ny), textcoords="offset points", xytext=(5, 5))

for i in range(n):
    fr = tour[i]
    to = tour[(i + 1) % n]
    plt.annotate("", xy=(nodes[to][0], nodes[to][1]),
                 xytext=(nodes[fr][0], nodes[fr][1]),
                 arrowprops=dict(arrowstyle="->", color="#e74c3c", lw=2))
plt.title("TSP Tour")
plt.savefig("tsp.png", dpi=150, bbox_inches="tight")
plt.show()
```

巡回路はノードを結ぶ赤い有向矢印で表示されます。
:::