以前作った4x4の数独を新しいWildqatに移植してみた

# 昨日くらいから流行ってます。
以前書いた[wildcatで4x4の数独を解いてみた](https://qiita.com/gyu-don/items/0c59423b12256821c38c)が、昨日くらいから「いいね」が増えはじめました。
Twitter検索してみると、有名な方に取り上げていただいたようです。大変有り難いことです。

ですが、その記事はwildcatが出たばかりの頃に書かれており、今は仕様が大きく変わっていて、最新のwildqat (wildcatから名称変更されました)では動きません。

とはいえ、アニーリングに使うQUBOは変わらないので、そこの部分は流用できます。
(wildqatに限らず、D-Waveや富士通のデジタルアニーラなどの実機にもQUBOは流用できます。ただし、特にD-Waveの場合、グラフ埋め込みなどの別の作業が必要になることがあります)

# 最新版で書き直してみた
ということで、2018年11月現在での最新版で、4x4の数独を書き直してみました。

## 忙しい人向け:

- APIが変わり、呼び出し方はシンプルになった
- 以前は何度も回さないと答えが合わなかったが、一発で合うようになった
- 今回のコードのJupyter Notebookは[ここ](https://gist.github.com/gyu-don/835755f2d2ed272651752afbe6f8fc77)にある

## 問題の定義

```py
problem_str = '''
1? ??
?? 2?

?3 ??
?? ?4'''

# Indices of array are:
#  0  1 |  2  3
#  4  5 |  6  7
# ------+-------
#  8  9 | 10 11
# 12 13 | 14 15

# Remove whitespace and newline
import re
problem_str = re.sub(r'[^0-9?]', '', problem_str)
```

前回同様、このように数独の問題を定義します。

## QUBOの作成

前回のを流用しています。新Wildqatでは、quboはnumpyのarrayではなく、listでとるので、あえてnumpy使う意義はないかもしれません。
が、面倒なので、numpyのままにしています。

### QUBOについて
QUBOと呼ばれる奇妙な二次元配列の読み方についても簡単に説明します。
QUBOを使ったアニーリングでは、最適化する対象はビットです。
QUBOの数式自体は、ここでは省略しますが、QUBOには、ビット同士の相互作用が埋め込まれており、あるビットと別のビットが両方とも1だった場合に得る、または失う「エネルギー」が記述されます。

そして、アニーリングをすると、できるだけ「エネルギー」を最小化するようなビットが求められます。
つまり、望ましいビットの組み合わせでは、QUBOの値をマイナスにして、望ましくないビットの組み合わせではQUBOの値をプラスにするように、QUBOを作ることで、望ましいビットの組み合わせをアニーリングで求められます。

今回の数独では、数字を{1, 2, 3, 4}から選ばないといけないのですが、1〜4を、4つのビットを使って表しています。

|ビット列|数字|
|:---|:---|
|1000| 1|
|0100| 2|
|0010| 3|
|0001| 4|

1100とか1111のように、複数のビットが1になっていたら、おかしなことになるのですが、そうなるのを防ぐように、そうなったときにペナルティを与える(最適解から遠ざかる)ようにします。

例えば、`qubo[1, 2] = 1.0`としてやると、
エネルギー = `qubo[1, 2] * 1ビット目 * 2ビット目` は次のようになります。

|1ビット目|2ビット目|エネルギー|
|:---|:---|:---|
|0|0|0.0|
|0|1|0.0|
|1|0|0.0|
|1|1|1.0|

このようにして、複数のビットが1になった場合のペナルティを定義することができます。

より詳細については、[前回の記事](https://qiita.com/gyu-don/items/0c59423b12256821c38c#qubo%E3%82%92%E4%BD%9C%E3%82%8B)をご覧ください。

## 今回使ったQUBO

4x4数独では、数字は16個あるので、今回は

```
# Bits arrangement:
# i = 0..15: i-th value is 1
# i = 16..31: (i-16)-th value is 2
# i = 32..47: (i-32)-th value is 3
# i = 48..63: (i-48)-th value is 4
```

のように、0〜15番目のビットは、i番目が1であることを表すビット、16〜31番めのビットは、(i-16)番目が1であることを表すビット、としています。

また、本来ならば、既に答えが埋まってるところはアニーリングする必要はないんですが、4x4なら問題自体が簡単なので、あえてできるだけ愚直な書き方をしています。

```py
import numpy as np

n = 4**3
grouping = [
    # rows
    [ 0,  1,  2,  3],
    [ 4,  5,  6,  7],
    [ 8,  9, 10, 11],
    [12, 13, 14, 15],
    # columns
    [ 0,  4,  8, 12],
    [ 1,  5,  9, 13],
    [ 2,  6, 10, 14],
    [ 3,  7, 11, 15],
    # blocks
    [ 0,  1,  4,  5],
    [ 2,  3,  6,  7],
    [ 8,  9, 12, 13],
    [10, 11, 14, 15],
]

# Bits arrangement:
# i = 0..15: i-th value is 1
# i = 16..31: (i-16)-th value is 2
# i = 32..47: (i-32)-th value is 3
# i = 48..63: (i-48)-th value is 4

# なにか数字を埋めると報酬(Hが小さくなる)がもらえるようにする
qubo = np.eye(n) * -1

# 行、列、ブロックに同じ数字があったら、ペナルティ(Hが大きくなる)を与えるようにする
for g in grouping:
    for i in range(4):
        for j in range(3):
            for k in range(j+1, 4):
                qubo[i*16 + g[j], i*16 + g[k]] += 2

# 同じマスには別の数字が入らないようにペナルティを与える
for i in range(16):
    for j in range(3):
        for k in range(j+1, k):
            qubo[i+j*16, i+k*16] += 2.5

# 既に問題文に数字が入っている場合、そこには別の数字を入れないように報酬やペナルティを調整する
for i, v in enumerate(problem_str):
    if v == '?':
        continue
    v = int(v) - 1
    qubo[v*16 + i, v*16 + i] -= 20
    for j in range(4):
        qubo[j*16 + i, j*16 + i] += 4
    for g in grouping:
        if i not in g:
            continue
        for j in g:
            qubo[v*16 + j, v*16 + j] += 3 + v*0.3
```

## アニーリングの実行

当初はWeb APIベースでやっていたのですが、今はWildqat自体でアニーリングできるようになっています。
以前よりもシンプルになりました。

```py
import wildqat as wq
a = wq.opt()
a.qubo = [...]  # QUBOをfloatのリストのリストで渡す
result = a.sa()
```

## 結果
```py
board = ['?'] * 16
for i,b in enumerate(result):
    if b == 1:
        if board[i%16] != '?':
            print('!', end='')
        board[i%16] = (i//16) + 1

print('''{}{} {}{}
{}{} {}{}

{}{} {}{}
{}{} {}{}'''.format(*board))
```

で、以下のように表示されました。

```
12 43
34 21

43 12
21 34
```

前回は、何度か回さないと正しい答えが得られなかったのですが、今回は一発で求まりました。
偶然正しかった可能性があるので、何度か回しましたが、いずれも間違った答えは出ませんでした。

また、[こちら](https://gist.github.com/gyu-don/835755f2d2ed272651752afbe6f8fc77)に今回のコードのJupyter Notebookを置いています。