21-4. 量子アニーリング（D-wave）でテント・アンド・ツリーパズルを解く（大規模編）

import numpy as np
from scipy import signal
from dwave.system.composites import EmbeddingComposite
from dwave.system.samplers import DWaveSampler
from pyqubo import Binary
from dwave_qbsolv import QBSolv

#============================
# パズル
#============================
pazzle = ['-----',
          '-t-t-',
          '-----',
          'tt---',
          '----t']
horizontal = [2,0,1,1,1]
vertical = [1,1,0,2,1]


#============================
# パズルの変換とフィルタの用意
#============================
tree = []
for line in pazzle:
    line = line.replace('-', '0')
    line = line.replace('t', '1')
    tree.append(list(line))
tree = np.array(tree, dtype=int)
print('tree\n{}'.format(tree))

tent = np.zeros(tree.shape, dtype=int)
print('tent\n{}'.format(tent))

horizontal = np.array(horizontal)
print('horizontal\n{}'.format(horizontal))
vertical = np.array(vertical)
print('vertical\n{}'.format(vertical))

#畳み込み用のフィルタ
filter_4 = np.array([[0, 1, 0],
                     [1, 0, 1],
                     [0, 1, 0]], dtype=int)
filter_8 = np.array([[1, 1, 1],
                     [1, 0, 1],
                     [1, 1, 1]], dtype=int)

#============================
# 量子ビットの用意
#============================
#木の４近傍のマス（＝テント候補地）
mask_candidate = (signal.convolve2d(tree, filter_4, mode='same') > 0) * (tree == 0)
print('mask_candidate\n{}'.format(mask_candidate))

#量子ビットの数
num_qbits = np.sum(mask_candidate)
print('num_qbits\n{}'.format(num_qbits))

#テント候補地にインデックスを振る（テント候補地以外は－１）
candidate_index = np.ones(tree.shape, dtype=int) * -1
candidate_index[mask_candidate==True] = range(num_qbits)
print('candidate_index\n{}'.format(candidate_index))

#量子ビットを用意する
for i in range(num_qbits):
    command = 'q{:03} = Binary(\'q{:03}\')'.format(i, i)
    print(command)
    exec(command)

#============================
# 条件を記述する
#============================
#条件を記述する関数
def H_add(qbits, num):
    command = '(' + '+'.join(['q{:03}'.format(i) for i in qbits]) + '-{})**2'.format(num)
    print(command)
    return eval(command)

#条件を記述していく
H = 0

#条件１、各木の４近傍のマスについて、１つか２つだけ１になる
#各木について
for i, j in np.array(np.where(tree==1)).T:
    #注目する木の座標
    print(i, j)
    
    #一つの木
    mask_single_tree = np.zeros(tree.shape, dtype=int)
    mask_single_tree[i, j] = 1
    
    #その周囲のテント候補地
    mask_around_single_tree = signal.convolve2d(mask_single_tree, filter_4, mode='same') * mask_candidate
    #print('mask_around_single_tree\n{}'.format(mask_around_single_tree))
    
    #そのインデックス
    qbits_index = candidate_index[mask_around_single_tree>0]
    #print('qbits_index\n{}'.format(qbits_index))
    
    #これらは１つか２つだけ１になる
    H += H_add(qbits_index, 1.5)

#条件２、８近傍で隣接する２つのマスのうち、０つまたは１つだけ１になる
#各テント候補地について
for i, j in np.array(np.where(mask_candidate==1)).T:
    #他のテント候補地が
    for k, l in np.array(np.where(mask_candidate==1)).T:
        #８近傍の関係であれば
        if (k, l) == (i, j+1) or (k, l) == (i+1, j-1) or (k, l) == (i+1, j) or (k, l) == (i+1, j+1):
            #注目する２つのテント候補地
            print((i, j), (k, l))
            
            #これらのインデックス
            qbits_index = [candidate_index[i, j], candidate_index[k, l]]
            #print(qbits_index)
            
            #これらは０つまたは１つだけ１になる
            H += H_add(qbits_index, 0.5)

#条件３、各行、各列の "テントの数" と "数字" が一致する
#vertical
for i in range(len(vertical)):
    #この列のインデックス
    qbits_index = candidate_index[i, :][candidate_index[i, :] > -1]
    print(qbits_index)
    
    #これらは指定された数だけ１になる
    H += H_add(qbits_index, vertical[i])

#horizontal
for i in range(len(horizontal)):
    #この列のインデックス
    qbits_index = candidate_index[:, i][candidate_index[:, i] > -1]
    print(qbits_index)
    
    #これらは指定された数だけ１になる
    H += H_add(qbits_index, horizontal[i])

#コンパイル
model = H.compile()
qubo, offset = model.to_qubo()
print('qubo\n{}'.format(qubo))

#D-waveの設定
sampler = EmbeddingComposite(DWaveSampler(endpoint='https://cloud.dwavesys.com/sapi',
                                          token='your_token',
                                          solver='DW_2000Q_5'))

#問題を分割してサンプリングしまくる
response = QBSolv().sample_qubo(qubo, solver=sampler, num_repeats=50)

#レスポンスの確認
count = 0
print('response')
for (sample, energy, occurrences) in response.data():
    print('Sample:{} Energy:{} Occurrences:{}'.format(list(sample.values()), energy, occurrences))
    
    #１番目を保存しておく
    if count == 0:
        best_sample = np.array(list(sample.values()))
    
    count += 1
    if count > 20:
        break

#最頻解の確認
print('best_sample\n{}'.format(best_sample))


#============================
# テントの位置を確認
#============================
tent[mask_candidate] = best_sample
print('tent\n{}'.format(tent))

#============================
# パズル
#============================
pazzle = ['---tt-----t-------',
          't-------tt---t--t-',
          't-----------t-t---',
          '-----t---t--------',
          '--t-t--------tt---',
          '-----t----------t-',
          '--t--t----t--t---t',
          '------t-t-t-----t-',
          '-tt---t-----t-----',
          '---t----t--t--t---',
          '---------t------t-',
          '------t--t------t-',
          '-t------------t---',
          '--t-t---t-----t---',
          '-----------t------',
          '---tt--t---------t',
          't-----------t-----',
          '----t-tt-t---tt-t-']
horizontal = [3,4,2,5,2,6,2,2,6,2,4,3,3,4,1,5,2,5]
vertical = [5,3,4,2,4,3,4,3,4,2,5,0,6,1,2,5,0,8]