量子アニーリングで今晩のおかずを最適化


食生活がガタガタなそこのあなた，たまには栄養を気にしてご飯を決めてみませんか．

## この記事の仕様
**量子アニーリングの実装を行ったことがない著者が，実装までのイメージを掴むために，解けると自分が嬉しい問題を選んで，試行錯誤した記録です**

- 私の基本的な興味は量子ゲート方式に傾いていますが，量子イジングマシン方式にも少し触れておいた方が可能性が広がるかも，という心意気です
    - ひとまずやってみて，なにかと痛い目を見て，基礎に帰結しちゃんと勉強をしはじめる，というサイクルが，経験的に自分の学習意欲を維持できるような気がしているので(今の所)，今回は「ひとまずやってみる」の段階です．
    - 自分の理解の整理のために，そしてツッコミを受け取れるように，記事にして残します．
- 荷造りをする予定もないし，セールスマンでもないので，ひとまず手軽な問題として，今日の晩ごはんを決めてみます
- 手探り状態なので，先輩方にツッコミをいただけると，とてもとても嬉しいです
- 量子アニーリングの理論に関しては，[「量子アニーリングの基礎」](https://www.kyoritsu-pub.co.jp/bookdetail/9784320035386)の5章の途中くらいまで読んだ程度です
- 著者は栄養学に関する知識は一切ありません

#### 参考にさせていただいたありがたい記事
- [以前作った4x4の数独を新しいWildqatに移植してみた](https://qiita.com/gyu-don/items/dbe005fd5099bfd8e946)
- [QUBOを使って日本地図の彩色問題やってみた](https://qiita.com/gyu-don/items/dc83392cb890db6c8d61)
- [プログラミング初心者がQUBOを使って個別指導塾の時間割を最適化してみる](https://qiita.com/papepupi/items/b7b84ff4ec69490fe70e)
- [整数ナップサック問題のIsing表現ハミルトニアン：Knapsack with Integer Weights](https://qard.is.tohoku.ac.jp/T-Wave/?p=69)

# 概要

**実現したいこと：健康的な今日の晩ごはんの組み合わせを探索する**

- ”健康的”の指標は，厚生労働省が出してる[「日本人の食事摂取基準(2015年版)」](https://www.glico.co.jp/navi/e07.html)(リンクは抜粋記事)に従う
    - 18〜29歳女性の項目を参照
-  1日の理想的な栄養素の比率がどうやら3：3：4らしい(ソースは[こちら](https://style.nikkei.com/article/DGXKZO16220180R10C17A5NZBP01?channel=DF130120166090))ので，これに倣い，健康的な晩ごはんには，1日に推奨されている栄養素の40％をふくむことにする(朝，昼は理想的な食事をしたということにして評価しない)
- 単位はkcalに限定する
    - 単位の変換は[こちらの資料](https://www.glico.co.jp/navi/dic/dic_01.html)に倣う
- 栄養素に関してはまったく詳しくないので，3大栄養素とよばれる「たんぱく質」，「脂質」，「炭水化物」のみを評価する
- [ちょうど良い資料](http://www.sej.co.jp/i/products/allergy/)があったので，組み合わせの対象はセブンイレブンのホットスナックに限定((健康的とは))
    - エントリーは以下の6つ(チョイスは私の好みです)

||料理名|たんぱく質[kcal]|脂質[kcal]|炭水化物[kcal]|
|:--|:-------|:-------|:-------|:-------|
|0|BIGポークフランク|11.5*4=46|27.4*9=246.6|7.4*4=29.6|
|1|ビッグアメリカンドッグ|6.4*4=25.6|16.0*9=144|37.3*4=149.2|
|2|からあげ棒(竜田揚げ)|7.1*4=28.4|13.4*9=120.6|13.8*4=55.2|
|3|チキンナゲット ５個入り|11.8*4=47.2|18.8*9=169.2|14.9*4=59.6|
|4|フライドポテト|2.4*4=9.6|8.2*9=73.8|26.9*4=107.6|
|5|牛肉コロッケ|3.7*4=14.8|16.5*9=148.5|18.7*4=74.8|
         
## 評価方法
- 選択した料理の，各栄養素の合計が，推奨値に近ければ近いほど健康的であるとする
- 評価するパラメータは，たんぱく質，脂質，炭水化物，そしてそ合計エネルギーの4つ
    - つまり，コスト関数が4つできる
    - 制約条件は特に設けない(栄養素の値に上限はなく，いくつ料理を選んでも構わない)
- 合計エネルギーはたんぱく質[kcal]+脂質[kcal]+炭水化物[kcal]で表す
- 同じ料理を2つ以上選ばない
    - バランスよくいろんなもの食べないとですよね(どれも揚げm)
   
### 各パラメータの推奨値

|パラメータ名|推奨値[kcal]|
|:-------|:-------|
|たんぱく質|20[g]*4＝80[kcal]|
|脂質|20~30[g]*9=180~270[kcal] → 240[kcal]|
|炭水化物|390~507[kcal] → 450[kcal]|
|エネルギー|780[kcal]|

# 計算
では，計算にかかります．量子アニーリングには，MDR社の[blueqat](https://github.com/Blueqat/Blueqat)を使用します．
今回作成したコードのJupyter Notebookは[こちら](https://nbviewer.jupyter.org/gist/nisuikiba/529edd0a5dd67ac0d4ea734924aae55b)からとべます．

## 目的関数(ハミルトニアン)を作成
各パラメータの目的関数を作成します．今回は，制約条件を定めないので，作成するのはコスト関数のみです．

### 計算の準備
まず，なにかと準備をします．

```python
from blueqat import opt
import numpy as np
from sympy import *

N = 6 #料理の数

P = 80 #たんぱく質の推奨値

L = 240 #脂質の推奨値

C = 450 #炭水化物の推奨値

E = 780 #エネルギーの推奨値

protain = [46, 25.6, 28.4, 47.2, 9.6, 14.8]
lipid = [246.6, 144, 120.6, 169.2, 73.8, 148.5]
carbon = [29.6, 149.2, 55.2, 59.6, 107.6, 74.8]

energy = []
for i in range (N):
    energy.append(round(protain[i]+lipid[i]+carbon[i], 1))

```
各料理のパラメータは配列に入れました．たとえば，`protain[0]`は，0番目の料理(BIGポークフランク)のたんぱく質の値を示しています．

$q0, q1, q2, q3, q4, q5 ∈ {0,1}$の6つの変数を用意し，$qi=0$のとき，$i$番目の料理を食べない，$qi=1$のとき，$i$番目の料理を食べる，とします．

```python
q0, q1, q2, q3, q4, q5 = symbols("q0 q1 q2 q3 q4 q5")
q = [q0, q1, q2, q3, q4, q5]
```
### コスト関数の作成
選択した料理の，各栄養素の合計が，推奨値に近ければ近いほど，コスト関数が小さい値になるような関数を考えます．以下のような式を考えました．

```math
(P-\sum_{i=0}^N p_{i}q_{i})^2+(L-\sum_{i=0}^N l_{i}q_{i})^2+(C-\sum_{i=0}^N c_{i}q_{i})^2+(E-\sum_{i=0}^N e_{i}q_{i})^2
```

これを計算してくれるコードを書きます．このへんのコードに関しては，本筋ではないのでここにはのせません(気になる方は[Jupyter Notebook](https://nbviewer.jupyter.org/gist/nisuikiba/529edd0a5dd67ac0d4ea734924aae55b)をご参照ください)．結果のみを示します．
`cost_p`はパラメータpのコスト関数(前述のコスト関数の第一項)を表しています．
尚，$qi$が$0$でも$1$でも$qi^2=qi$なので，$qi^2=qi$を代入しています．

```python
cost_p = 2355.2*q0*q1 + 2612.8*q0*q2 + 4342.4*q0*q3 + 883.2*q0*q4 + 1361.6*q0*q5 - 5244*q0 + 1454.08*q1*q2 + 2416.64*q1*q3 + 491.52*q1*q4 + 757.76*q1*q5 - 3440.64*q1 + 2680.96*q2*q3 + 545.28*q2*q4 + 840.64*q2*q5 - 3737.44*q2 + 906.24*q3*q4 + 1397.12*q3*q5 - 5324.16*q3 + 284.16*q4*q5 - 1443.84*q4 - 2148.96*q5 + 6400
cost_l = 99249.92*q0*q1 + 87194.8*q0*q2 + 126637.28*q0*q3 + 48805.68*q0*q4 + 95563.16*q0*q5 - 54833.24*q0 + 50540.8*q1*q2 + 73402.88*q1*q3 + 28289.28*q1*q4 + 55391.36*q1*q5 - 52643.84*q1 + 64487.2*q2*q3 + 24853.2*q2*q4 + 48663.4*q2*q5 - 49319.0*q2 + 36095.52*q3*q4 + 70676.24*q3*q5 - 57043.04*q3 + 27238.44*q4*q5 - 33076.44*q4 - 51717.11*q5 + 57600
cost_c = 205434.72*q0*q1 + 131586.48*q0*q2 + 177854.4*q0*q3 + 123080.4*q0*q4 + 153431.64*q0*q5 - 186167.16*q0 + 130197.92*q1*q2 + 175977.6*q1*q3 + 121781.6*q1*q4 + 151812.56*q1*q5 - 185286.56*q1 + 112718.4*q2*q3 + 78004.4*q2*q4 + 97240.04*q2*q5 - 142082.36*q2 + 105432.0*q3*q4 + 131431.2*q3*q5 - 172224.0*q3 + 90954.2*q4*q5 - 135419.0*q4 - 157598.39*q5 + 202500
cost_e = 205434.72*q0*q1 + 131586.48*q0*q2 + 177854.4*q0*q3 + 123080.4*q0*q4 + 153431.64*q0*q5 - 398819.16*q0 + 130197.92*q1*q2 + 175977.6*q1*q3 + 121781.6*q1*q4 + 151812.56*q1*q5 - 395694.56*q1 + 112718.4*q2*q3 + 78004.4*q2*q4 + 97240.04*q2*q5 - 276854.36*q2 + 105432.0*q3*q4 + 131431.2*q3*q5 - 354384.0*q3 + 90954.2*q4*q5 - 261479.0*q4 - 314744.39*q5 + 608400
```

今回は，制約条件を定めていないので，このコスト関数がそのまま目的関数になります．

## QUBOを作成
では，先ほど作成した目的関数をQUBOに変換します．

```python
#6*6の行列を作成
qubo_p = np.zeros((6,6))
qubo_l = np.zeros((6,6))
qubo_c = np.zeros((6,6))
qubo_e = np.zeros((6,6))
```
$6\times6$の行列の上三角部分に，先ほど計算した目的関数の，該当する係数部分を代入します．例によって詳しい計算過程は長くなるので掲載しません．
最終的には以下のような行列が得られました．

```python
[[-645063.56  512474.56  352980.56  486688.48  295849.68  403788.04]
 [      0.   -637065.6   312390.72  427774.72  272344.    359774.24]
 [      0.         0.   -471993.16  292604.96  181407.28  243984.12]
 [      0.         0.         0.   -588975.2   247865.76  334935.76]
 [      0.         0.         0.         0.   -431418.28  209431.  ]
 [      0.         0.         0.         0.         0.   -526208.85]]
```

## 量子アニーリングする
[blueqat](https://github.com/Blueqat/Blueqat)を使用します．

```python
a = opt.opt() 
a.qubo = qubo
result = a.sa(shots = 100, sampler = "fast")
opt.counter(result)
```
結果がこちら(どきどき)．

```python
Counter({'000110': 19,
         '010010': 4,
         '010001': 9,
         '110000': 6,
         '100001': 7,
         '010100': 8,
         '100100': 10,
         '000101': 11,
         '001011': 10,
         '011000': 9,
         '101000': 6,
         '100010': 1})
```

ばらつきがすごい…！
でも答えっぽいものが出てきて感動です．

数回実行してみたところ，かろうじて頻度が高かった`'000110'`と`'001011'`を最終的な答えとします．
精度をあげるにはどうしたらいいんだろう…．

# 晩ごはんが決まりました！
さあ答え合わせです．
`'000110'`が示しているのは，**「チキンナゲット ５個入り」**と**「フライドポテト」**の組み合わせです．
この組み合わせの栄養素は以下の通りでした．

|パラメータ名|推奨値|`'000110'`の値|
|:-------|:-------|:-------|
|たんぱく質|80|56.8|
|脂質|240|248|
|炭水化物|450|167.2|
|エネルギー|780|472|

微妙…．

もう一方の答えに期待しましょう．
`'001011'`が示しているのは，**「からあげ棒(竜田揚げ)」**と**「フライドポテト」**と**「牛肉コロッケ」**の組み合わせです．
この組み合わせの栄養素は以下の通りでした．

|パラメータ名|推奨値|`'001011'`の値|
|:-------|:-------|:-------|
|たんぱく質|80|52.8|
|脂質|240|342.9|
|炭水化物|450|237.6|
|エネルギー|780|633.3|

なんだか少し見えてきました．

**脂質が高い…．**

どちらの組み合わせも，脂質のみが推奨値を超えており，それ以外のパラメータがほぼ推奨値を大きく下回っています．つまり，用意した料理の多くが，脂質が非常に高く(ホットスナックだから当たり前)，脂質のみが早い段階で推奨値以上になってしまったため，ほかの栄養素が推奨値よりかなり小さい値になってしまっているようです．

# まとめ
用意したデータセットに偏りがあったため，期待していたような綺麗な結果はあまり得られませんでした…．しかし，試行錯誤しながら手を動かすことで，少しだけ量子アニーリングという不思議な技術が怖くなくなりました．
実際のところ，ちゃんと計算できているのか，目的関数の設定の方法は適切であったのか，初心者の自分には定かではない点がまだいくつかあります．是非とも，先輩方にツッコミを入れていただければと思います．

あと，ホットスナックはとっても脂質が高いということを，身をもって知ることができました．
ホットスナックだけで理想的な食事を設計することはイジングマシンにも不可能なので，みなさんもホットスナックでごはんを済ませないようにしましょう！