wildcatで巡回セールスマン問題を試してみた

#はじめに
wildcatの使い方は https://github.com/mdrft/wildcat を見てください。
こちらも参考になるかもしれません。[wildcatで4x4の数独を解いてみた](https://qiita.com/gyu-don/items/0c59423b12256821c38c)
また、巡回セールスマン問題の詳細はこちらの記事が参考になります。[物理のいらない量子アニーリング入門](http://blog.brainpad.co.jp/entry/2017/04/20/160000)
それでは標題の通り試してみます。
#ハミルトニアンの作成
参考[[2](http://blog.brainpad.co.jp/entry/2017/04/20/160000)]から、巡回セールスマン問題のハミルトニアンは下記の様になります。

```math 
H = \sum_{t,a,b}D_{a,b}x_{t,a}x_{t+1,b} + A\sum_{t}\left(\sum_{a}x_{t,a}-1\right)^2 + A\sum_{a}\left(\sum_{t}x_{t,a}-1\right)^2
```
また巡回する町の数を $n$ とします。
#QUBOに変形
ハミルトニアンを下記の形式に変形します。

```math 
H = \sum_{i,j}Q_{i,j}x_{i}x_{j}
```
ここで $H$ を見ると $x$ の一次の項が出そうですが、$x$ が0か1しか値が取らないという性質を利用して、

```math 
x_{t,a} = x_{t,a}^2
```
とすることで、目的の形式に持っていきます。
##ステップ１
いきなりは難しいので、まず下記の形式に持っていきます。

```math 
H = \sum_{t_1,t_2,a,b}J_{t_1,t_2,a,b}x_{t_1,a}x_{t_2,b}
```

###第一項

```math 
\begin{eqnarray}
\sum_{t,a,b}D_{a,b}x_{t,a}x_{t+1,b} = \sum_{t_1,t_2,a,b}\delta_{t_1+1,t_2}D_{a,b}x_{t_1,a}x_{t_2,b}
\end{eqnarray}
```
ここで $\delta_{a,b}$ はクロネッカーデルタです。
###第二項

```math 
\begin{eqnarray}
A\sum_{t}\left(\sum_{a}x_{t,a}-1\right)^2 = A\sum_{t_1,t_2,a,b}\left(\delta_{t_1,t_2}-2\delta_{a,b}\delta_{t_1,t_2}\right)x_{t_1,a}x_{t_2,b} + An
\end{eqnarray}
```
ここで $x_{t,a} = x_{t,a}^2$ を使いました。

###第三項

```math 
\begin{eqnarray}
A\sum_{a}\left(\sum_{t}x_{t,a}-1\right)^2 = A\sum_{t_1,t_2,a,b}\left(\delta_{a,b}-2\delta_{a,b}\delta_{t_1,t_2}\right)x_{t_1,a}x_{t_2,b} + An
\end{eqnarray}
```
ここで $x_{t,a} = x_{t,a}^2$ を使いました。
したがって、$H$ は下記の様になりました。

```math 
H=\sum_{t_1,t_2,a,b}\left(\delta_{t_1+1,t_2}D_{a,b}+A\delta_{t_1,t_2}+A\delta_{a,b}-4A\delta_{a,b}\delta_{t_1,t_2}\right)x_{t_1,a}x_{t_2,b} + 2An
```

ここで、定数項はエネルギーの基準をずらしてゼロにして考えます。あとは $J$ を計算して、それを変形して $Q$ を計算します。これで数学の計算は終わりにします。

#コード
それではコードで書いていきます。
まず、ライブラリを読み込みます。そして断熱計算の設定を行います。今回はローカルで計算をします。また、計算が成功する確率を上げるため、初期の温度を1000度にしました。

```python
from wildcat.solver.qubo_solver import QuboSolver
from wildcat.network.local_endpoint import LocalEndpoint
from wildcat.annealer.simulated.simulated_annealer import SimulatedAnnealer
from wildcat.annealer.simulated.single_spin_flip_strategy import SingleSpinFlipStrategy
from wildcat.annealer.simulated.temperature_schedule import TemperatureSchedule
from wildcat.util.matrix import hamiltonian_energy
    
schedule = TemperatureSchedule(initial_temperature=1000, last_temperature=0.1, scale=0.8)
strategy = SingleSpinFlipStrategy(repetition=10)
annealer = SimulatedAnnealer(schedule=schedule, strategy=strategy)
local_endpoint = LocalEndpoint(annealer=annealer)
```
次に参考[[2](http://blog.brainpad.co.jp/entry/2017/04/20/160000)]を参考に町の位置と距離の計算を定義し、また定数は適当に400としました。

```python
#座標
positions = np.array((
    (24050.0000, 123783),
    (24216.6667, 123933),
    (24233.3333, 123950),
    (24233.3333, 124016),
    (24250.0000, 123866),
    (24300.0000, 123683),
    (24316.6667, 123900),
    (24316.6667, 124083),
    (24333.3333, 123733),
))
#距離の計算
def dist(a, b):
    a = np.array(a)
    b = np.array(b)
    return np.sqrt(((a - b)**2).sum())
#町の数
N = len(positions)
#定数A
A = 400
```
それでは一気に $J$ を計算していきます。


```python
#行列Dを計算
D = np.empty((N, N), dtype = np.float64)
for i in range(N):
    for j in range(N):
        D[i, j] = dist(positions[i], positions[j])
#クロネッカーデルタを定義
def delta(i, j):
    if(i == j):
        return 1
    else:
        return 0
#行列Jを計算
J = np.empty((N, N, N, N), dtype = np.float64)
for a in range(N):
    for t1 in range(N):
        for b in range(N):
            for t2 in range(N):
                J[a, t1, b, t2] = delta(t1+1, t2) * D[a, b] + A * delta(t1, t2) + A * delta(a, b)\
                                  - 4 * A * delta(a, b) * delta(t1, t2)
#定数項の計算
C = 2 * A * N
```

それでは式を変形して $Q$ を$J$から計算していきます。

```python
#Qの計算
Q = np.empty((N*N,N*N), dtype = np.float64)
x = 0
for a in range(N):
    for t1 in range(N):
        y = 0
        for b in range(N):
            for t2 in range(N):
                Q[x,y] = J[a,t1,b,t2]
                y = y + 1
        x = x + 1
```
$Q$ が出来たので、実際に最適化の計算を行います。

```python
#最適化計算
solver = QuboSolver(Q)

def callback(arrangement):
    e = solver.hamiltonian_energy(arrangement)
    print("Energy: ", solver.hamiltonian_energy(arrangement))
    print("Spins: ", arrangement)

arrangement = solver.solve(callback, endpoint=local_endpoint).result()

[Out]

Energy:  -179479.61895467353
Spins:  [-1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1
 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1 -1
 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1 -1
 -1 -1 -1  1 -1 -1 -1 -1 -1]
```
ここでEnergyは低いほうが正解近くなります。またSpinsは量子ビットのスピンの向きを表しています。
結果をわかりやすくすると下記のようになります。

```python
result = np.empty((N, N), dtype = np.int)
i = 0
for x in range(N):
    for y in range(N):
        if(arrangement[i] == 1):
            result[x, y] = 1
        else:
            result[x, y] = 0
        i += 1
print("Result: \n", result)

[Out]

Result: 
 [[0 1 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 1]
 [0 0 0 0 0 1 0 0 0]
 [0 0 1 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [1 0 0 0 0 0 0 0 0]
 [0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 0 1 0 0]
 [0 0 0 1 0 0 0 0 0]]
```
見てみると、巡回の条件を満たしていない部分があることがわかります。
もう一度やってみると、今度は成功しました。

```python
Result: 
 [[1 0 0 0 0 0 0 0 0]
 [0 1 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 1 0]
 [0 0 0 0 0 1 0 0 0]
 [0 0 0 1 0 0 0 0 0]
 [0 0 1 0 0 0 0 0 0]
 [0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 0 0 0 1]
 [0 0 0 0 0 0 1 0 0]]
```
成功したときの距離を計算すると下記のようになります。

```python
#距離を計算する
length = 0
for t in range(N):
    for a in range(N):
        for b in range(N):
            if(t < N-1):
                if(result[a, t] == 1 and result[b, t+1] == 1):
                    length += D[b, a]
            else:
                for c in range(N):
                    if(result[a, t] == 1 and result[c, 0] == 1):
                        length += D[a, c]
print("Length:", length)

[Out]

Length: 5203.634602384536
```

興味がある方はぜひ試してみてください。コードは上げておきます。https://github.com/Yokomakura/wildcat-TSP/blob/master/wildcat%E5%B7%A1%E5%9B%9E%E3%82%BB%E3%83%BC%E3%83%AB%E3%82%B9%E3%83%9E%E3%83%B3.ipynb
急ぎ足でやったので、ミスってたらごめんなさい。
#参考
1. [wildcatで4x4の数独を解いてみた](http://blog.brainpad.co.jp/entry/2017/04/20/160000)
2. [物理のいらない量子アニーリング入門](http://blog.brainpad.co.jp/entry/2017/04/20/160000)