第8回: アニーリング基礎7. PyquboのDecode活用

これまで色々なアニーリング問題を整理してきたが、**PyquboのDecodeクラス**を有効活用できていなかったので、今回整理してみる。

**前回**: [第7回: アニーリング基礎6.補助スピン活用基礎2（シフト最適化問題）](https://qiita.com/opt_d/items/3d1eaebb419faf4222ff)
**参考**: [公式チュートリル](https://pyqubo.readthedocs.io/en/latest/reference/model.html#decodedsample)
**Keyword** ：Pyqubo, Decode

### 1. 基本準備/設定
#### 1.1 量子アニーリング用モジュールのインストール
まず、初めに、以下のモジュールをインストールする。
- ***dwave-ocean-sdk***：Dwaveアニーリングマシンのモジュール
- ***openjij***：アニーリングシュミレーター
- ***pyqubo***：Qubo行列作成サポートモジュール
- ***optuna***：ハイパーパラメタの最適化

```Python
pip install dwave-ocean-sdk
pip install openjij
pip install pyqubo
pip install optuna
```

#### 1.2 Dwaveの初期設定と基本的なライブラリのインポート
続いて、Dwaveのtoken、enpoint(アクセスURL)、使用するマシンの指定まで実施。併せて、Jijシュミレーターも設定する。
```Python
###基本設定
token = "DEV-73ac79954ef21XXXXXXXXXXXXXXXXXXXXXX"
endpoint = "https://cloud.dwavesys.com/sapi/"
solver = "Advantage2_prototype1.1"
###Dwaveの設定
from dwave.system import DWaveSampler, EmbeddingComposite
dw_sampler = DWaveSampler(solver=solver, token=token, endpoint = endpoint)
###Jijシュミレーターの設定
from openjij import SQASampler
sampler_jij = SQASampler()
###基本的なライブラリのインポート
import numpy as np
import pandas as pd
import seaborn as sns
```
### 2. 問題設定、定式化、最適化計算
#### 2.1 問題
複数の制約式があり、エネルギー式が明確な以下の問題を例にとって、PyquboのDecodeクラスを整理する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2668554/41893fb4-ddda-66d0-2eb4-7b6f02f5f8ad.png)
スコア表の値$S$、スコア表の行数$R$、列数$C$を以下のように定義する。
また、選択箇所を示すスピン$x_{rc}$も設定する。$r行c列$を選択する時、$x_{rc}=1$となる。
```Python
S = [[10,20,15],[30,40,35],[50,60,55]]
R, C = len(S),len(S)
from pyqubo import Array
x = Array.create(name="x", shape=(R,C), vartype="BINARY")
```
#### 2.2 定式化
エネルギ-式$E(x)$は選択したスコアの合計値として、制約1、2併せて以下のように定式化できる。
```math
E(x)=-\sum_{r=1}^R\sum_{c=1}^CS_{rc}x_{rc}\\
C_{1}(x) = \sum_{c=1}^C(\sum_{r=1}^Rx_{rc}-1)^2\\
C_{2}(x) = \sum_{r=1}^R(\sum_{c=1}^Cx_{rc}-1)^2\\
H(x)=\lambda_{e1}E(x)+\lambda_{c1}C_{1}(x)+\lambda_{c2}C_{2}(x)
```
```Python
### 第1項の定式化
E = 0
for r in range(R):
  for c in range(C):
    E = E - S[r][c]*x[r][c]

### 第2項の定式化
C1 = 0
for c in range(C):
  _c = 0
  for r in range(R):
    _c = _c + x[r][c]
  C1 = C1 + (_c-1)**2

### 第3項の定式化
C2 = 0
for r in range(R):
  _c = 0
  for c in range(C):
    _c = _c + x[r][c] 
  C2 = C2 + (_c-1)**2

### lambdaの設定
from pyqubo import Placeholder
lambda_e1,lambda_c1,lambda_c2 = Placeholder("lambda_e1"),Placeholder("lambda_c1"),Placeholder("lambda_c2")
```
PyquboのDecodeクラスは、Dwaveで得られた結果を見やすく整理してくれる。具体的には、`.sample`スピン毎の結果、`.energy`エネルギー値、`.constraints()`/`.subh`制約成立有無、を取得することが出来る。ただし、事前準備として、ハミルトニアンを作成する時に、エネルギー式、制約式を区別する必要がある。Pyquboの制約式をラベリングするConstraint関数により区別が可能。コーディングは以下の通り。`Constraint(制約式,"制約式名")`でラベリングされる。尚、ラベリングされなかったものがエネルギー式として認識される。
```Python
### 制約式をラベリングするConstraintクラスのインポート
from pyqubo import Constraint

### ハミルトニアンの作成
H = lambda_e1*E + lambda_c1*Constraint(C1,"C1") + lambda_c2*Constraint(C2,"C2")
```
#### 2.3 最適化計算
Qubo行列の取出しから、求解まで一挙に進める。
```Python
###Qubo行列の取り出し
model = H.compile()
feed_dict = {"lambda_e1":1,"lambda_c1":100,"lambda_c2":100}
Qubo,offset = model.to_qubo(feed_dict=feed_dict)
###Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
result = sampler_dw.sample_qubo(Qubo, num_reads=10)

result.record
----------------------------
<Output>
rec.array([([1, 0, 0, 0, 0, 1, 0, 1, 0], -705., 1, 0.),
           ([0, 1, 0, 1, 0, 0, 0, 0, 1], -705., 1, 0.),
           ([0, 0, 1, 0, 1, 0, 1, 0, 0], -705., 6, 0.),
           ([0, 1, 0, 0, 0, 1, 1, 0, 0], -705., 2, 0.)],
          dtype=[('sample', 'i1', (9,)), ('energy', '<f8'), ('num_occurrences', '<i8'), ('chain_break_fraction', '<f8')])
```
### 3. PyquboのDecodeの検証
#### 3.1 Decodeの使い方
decode_sampleset関数にDwaveでの計算結果と、$\lambda$を渡すことでdecodeされた結果を得ることができる。
```Python
decoded_result = model.decode_sampleset(sampleset=result,feed_dict=feed_dict)
decoded_result
----------------------------
<Output>
[DecodedSolution({x[0][1]:0, x[0][0]:0, x[1][0]:0, x[2][1]:0, x[0][2]:1, x[2][2]:0, x[1][2]:0, x[2][0]:1, x[1][1]:1}, energy=-105.000000),
 DecodedSolution({x[0][1]:1, x[0][0]:0, x[1][0]:0, x[2][1]:0, x[0][2]:0, x[2][2]:0, x[1][2]:1, x[2][0]:1, x[1][1]:0}, energy=-105.000000),
 DecodedSolution({x[0][1]:1, x[0][0]:0, x[1][0]:1, x[2][1]:0, x[0][2]:0, x[2][2]:1, x[1][2]:0, x[2][0]:0, x[1][1]:0}, energy=-105.000000),
 DecodedSolution({x[0][1]:0, x[0][0]:1, x[1][0]:0, x[2][1]:1, x[0][2]:0, x[2][2]:0, x[1][2]:1, x[2][0]:0, x[1][1]:0}, energy=-105.000000)]
```
試しに1回目の結果を、decode有り無しで比較すると以下の通り。decodeすることにより、スピンの座標と値の双方が示されるとともに、ハミルトニアン合計値ではなく、エネルギー値のみが返される。
```Python
print(result.record[0])
print(decoded_result[0])
----------------------------
<Output>
([0, 0, 1, 0, 1, 0, 1, 0, 0], -705., 6, 0.)
DecodedSolution({x[0][1]:0, x[0][0]:0, x[1][0]:0, x[2][1]:0, x[0][2]:1, x[2][2]:0, x[1][2]:0, x[2][0]:1, x[1][1]:1}, energy=-105.000000)
```
得られた結果から、`.sample`スピン毎の結果、`.energy`エネルギー値、`.constraints()`/`.subh`制約確認、を取得することが出来る。
```Python
decoded_result[0].sample
----------------------------
<Output>
{'x[0][1]': 0,
 'x[0][0]': 0,
 'x[1][0]': 0,
 'x[2][1]': 0,
 'x[0][2]': 1,
 'x[2][2]': 0,
 'x[1][2]': 0,
 'x[2][0]': 1,
 'x[1][1]': 1}
```

```Python
decoded_result[0].energy
----------------------------
<Output>
-105.0
```

```Python
decoded_result[0].constraints()
----------------------------
<Output>
{'C2': (True, 0.0), 'C1': (True, 0.0)}
```

```Python
decoded_result[0].subh
----------------------------
<Output>
{'C1': 0.0, 'C2': 0.0}
```
スピン毎の結果を示す`.sample`の出力は見難いので、行列で書き直す。
```Python
def ans(i, decoded_result):
  spin = decoded_result[i].sample
  ans = np.zeros(R*C).reshape(R,C) #結果を代入する3x3の空行列を作成
  for key,value in zip(spin.keys(), spin.values()):
      key = key.replace("]","").split("[") # x[0][1]→["x","0","1"]に変換
      key = [int(_v) for _v in key[1:]] # ["x","0","1"]→ [0,1]に変換
      ans[key[0]][key[1]]=value #結果の代入
  return ans
ans(1,decoded_result)
----------------------------
<Output>
array([[0., 1., 0.],
       [0., 0., 1.],
       [1., 0., 0.]])
```
制約確認も使用しやすく整理する。

```Python
def judge(i,decoded_result):
  judge_result=[]
  judge = decoded_result[i].constraints()
  for key in reversed(judge.keys()): #Constraints()は逆順のリストを返すため、reversedで逆転させて降順にする
    _temp = 0
    if judge[key][0] == True: #制約を満たす場合、つまりTrueが返される場合、「1」を返す。そうでない場合は「0」
      _temp = 1
    judge_result.append(_temp)
  return judge_result
judge(1,decoded_result)
----------------------------
<Output>
[1, 1]
```

必要な情報をひとまとめにする関数を定義する。
```Python
def summary(i, decoded_result):
  return [abs(decoded_result[i].energy), judge(i,decoded_result), ans(i,decoded_result)]
summary(1,decoded_result)
----------------------------
<Output>
[105.0, [1, 1], array([[0., 1., 0.],[0., 0., 1.],[1., 0., 0.]])]
```

#### 3.2 検証1: 制約条件を満たさない場合の挙動
$\lambda$の値を変えながら、decodeの挙動を確認してみる。

```Python
###λの調整
feed_dict = {"lambda_e1":1,"lambda_c1":100,"lambda_c2":10}

###Qubo行列の取り出し
model = H.compile()
Qubo,offset = model.to_qubo(feed_dict=feed_dict)

###Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
result = sampler_dw.sample_qubo(Qubo, num_reads=10)

###Decode
decoded_result = model.decode_sampleset(sampleset=result,feed_dict=feed_dict)

for i in range(len(result)):
  print(summary(i,decoded_result))
----------------------------
<Output>
[125.0, [1, 0], array([[0., 0., 0.],[0., 1., 0.],[1., 0., 1.]])]
[125.0, [1, 0], array([[0., 0., 0.],[1., 0., 0.],[0., 1., 1.]])]
[105.0, [1, 0], array([[0., 0., 0.],[1., 1., 0.],[0., 0., 1.]])]
[105.0, [1, 0], array([[0., 1., 0.],[0., 0., 0.],[1., 0., 1.]])]
[105.0, [1, 1], array([[0., 1., 0.],[1., 0., 0.],[0., 0., 1.]])]
[105.0, [1, 1], array([[1., 0., 0.],[0., 1., 0.],[0., 0., 1.]])]
```
制約条件が不成立時には確かにFalseを返していることが確認できた。

#### 3.3 検証2: エネルギー式が複数ある場合の挙動
エネルギー式が複数存在する場合でも、PyquboのConstraintで制約式とラベリングしたもの以外がエネルギー式として認識されるはずである。上の定式化のE(x)を単純に複製して、エネルギー値が2倍なっているか確認してみる。
```Python
### ハミルトニアンの作成
lambda_e2 = Placeholder("lambda_e2")
H = lambda_e1*E + lambda_c1*Constraint(C1,"C1") + lambda_c2*Constraint(C2,"C2") + lambda_e2*E #λのみ新たに定義してエネルギー式を複製

###λの調整
feed_dict = {"lambda_e1":1,"lambda_c1":100,"lambda_c2":100, "lambda_e2":1}

###Qubo行列の取り出し
model = H.compile()
Qubo,offset = model.to_qubo(feed_dict=feed_dict)

###Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
result = sampler_dw.sample_qubo(Qubo, num_reads=10)

###Decode
decoded_result = model.decode_sampleset(sampleset=result,feed_dict=feed_dict)

for i in range(len(result)):
  print(summary(i,decoded_result))
----------------------------
<Output>
[210.0, [1, 1], array([[0., 1., 0.],[0., 0., 1.],[1., 0., 0.]])]
[210.0, [1, 1], array([[0., 1., 0.],[1., 0., 0.],[0., 0., 1.]])]
```
確かにエネルギー式が2倍になっていることが確認できた。

#### 3.4 検証3: 制約条件が「-1」でない場合の挙動
例題では$(\sum_{}x_{rc}-1)^2$としたが、$(\sum_{}x_{rc}-2)^2$と2カ所変更する場合の`.constraints()`の挙動を確認してみる。
```Python
### 第1項の定式化
E = 0
for r in range(R):
  for c in range(C):
    E = E - S[r][c]*x[r][c]

### 第2項の定式化
C1 = 0
for c in range(C):
  _c = 0
  for r in range(R):
    _c = _c + x[r][c]
  C1 = C1 + (_c-2)**2 # 「-1」から「-2」に変更した

### 第3項の定式化
C2 = 0
for r in range(R):
  _c = 0
  for c in range(C):
    _c = _c + x[r][c] 
  C2 = C2 + (_c-2)**2 # 「-1」から「-2」に変更した

### lambdaの設定
from pyqubo import Placeholder
lambda_e1,lambda_c1,lambda_c2 = Placeholder("lambda_e1"),Placeholder("lambda_c1"),Placeholder("lambda_c2")

### ハミルトニアンの作成
H = lambda_e1*E + lambda_c1*Constraint(C1,"C1") + lambda_c2*Constraint(C2,"C2")

###λの調整
feed_dict = {"lambda_e1":1,"lambda_c1":10,"lambda_c2":10} #制約式のλに小さい値10を入れる

###Qubo行列の取り出し
model = H.compile()
Qubo,offset = model.to_qubo(feed_dict=feed_dict)

###Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
result = sampler_dw.sample_qubo(Qubo, num_reads=10)

###Decode
decoded_result = model.decode_sampleset(sampleset=result,feed_dict=feed_dict)

for i in range(len(result)):
  print(summary(i,decoded_result),decoded_result[i].constraints()) #制約式の成立可否の生データを追加した。
----------------------------
<Output>
[265.0, [0, 0], array([[0., 1., 1.],[1., 1., 1.],[1., 1., 1.]])]
{'C2': (False, 2.0), 'C1': (False, 2.0)}
```
制約条件の数値が変わっても`.constraints()`が機能していることを確認。一方で(False, 2.0)の(2.0)が何を示すかが不透明。

#### 3.5 一連の流れの関数化
```Python
###Decode
decoded_result = model.decode_sampleset(sampleset=result,feed_dict=feed_dict)

###結果の可視化
def ans(i, decoded_result):
  spin = decoded_result[i].sample
  ans = np.zeros(R*C).reshape(R,C) #結果を代入する3x3の空行列を作成
  for key,value in zip(spin.keys(), spin.values()):
      key = key.replace("]","").split("[") # x[0][1]→["x","0","1"]に変換
      key = [int(_v) for _v in key[1:]] # ["x","0","1"]→ [0,1]に変換
      ans[key[0]][key[1]]=value #結果の代入
  return ans

###制約達成の可視化
def judge(i,decoded_result):
  judge_result=[]
  judge = decoded_result[i].constraints()
  for key in reversed(judge.keys()): #Constraints()は逆順のリストを返すため、reversedで逆転させて降順にする
    _temp = 0
    if judge[key][0] == True: #制約を満たす場合、つまりTrueが返される場合、「1」を返す。そうでない場合は「0」
      _temp = 1
    judge_result.append(_temp)
  return judge_result

###結果の表示
def summary(i, decoded_result):
  return [abs(decoded_result[i].energy), judge(i,decoded_result), ans(i,decoded_result)]
```

### おわりに
これまでPyquboのDecodeを使用できていなかったが、制約達成の確認、エネルギー値の取得は極めて便利なことがよくわかった。一方、制約達成をどのように確認しているかは不透明。時間ある時にメカニズムを確認してみたい。。。