21-20. 量子アニーリング（QUBO）で天体ショー（Galaxies）を解く

import numpy as np
from pyqubo import Binary
from dwave_qbsolv import QBSolv

#パラメータ
h = 5
w = 5
dot_num = 6
dots = [(0, 1.5), (1, 2), (2, 2), (3, 2), (4, 0), (4, 3)]


#量子ビットを用意する
for i in range(h):
    for j in range(w):
        for k in range(dot_num):
            command = f'q{i}{j}{k} = Binary(\'q{i}{j}{k}\')'
            print(command)
            exec(command)

#各マスは一つのレイヤーに所属する（＝kの次元がワンホット）
H = 0
for i in range(h):
    for j in range(w):
        command = f'H += (q{i}{j}0 + q{i}{j}1 + q{i}{j}2 + q{i}{j}3 + q{i}{j}4 + q{i}{j}5 - 1)**2'
        print(command)
        exec(command)

#各レイヤー繰り返し
for n in range(dot_num):
    #点の座標
    d = dots[n]
    #各マスをチェックして
    for i in range(h):
        for j in range(w):
            target = (int(d[0] + (d[0] - i)), int(d[1] + (d[1] - j)))
            #点と一致していればそのマスは「1にする」（半分重なりに注意）
            if -0.5 <= (i - d[0]) <= 0.5 and -0.5 <= (j - d[1]) <= 0.5:
                command = f'H += (q{i}{j}{n} - 1)**2'
                print(command)
                exec(command)
            #点を挟んで対称のマスが存在すればその2マスは「2個の量子ビットを揃える」
            elif 0 <= target[0] < h and 0 <= target[1] < w:
                """
                ここ注意、2回設定されるので重みは各0.5にすること
                """
                command = f'H += 0.5 * (q{i}{j}{n} - q{target[0]}{target[1]}{n})**2'
                print(command)
                exec(command)
            else:
                #相手がいなければそのマスは無関係なので「0にする」
                command = f'H += (q{i}{j}{n} - 0)**2'
                print(command)
                exec(command)

#コンパイル
model = H.compile()
qubo, offset = model.to_qubo()

#分割サンプリング
response = QBSolv().sample_qubo(qubo, num_repeats=500)

#レスポンスの確認
print('response')
for (sample, energy, occurrence) in response.data():
    print(f'Sample:{list(sample.values())} Energy:{energy} Occurrence:{occurrence}')
    
    #レイヤー番号に戻して確認
    box = np.array(list(sample.values())).reshape(h, w, dot_num)
    ans = np.zeros((h, w), int)
    for i in range(h):
        for j in range(w):
            ans[i, j] = np.argmax(box[i, j, :])
    print(ans)

response
Sample:[0, 1, 0,...,0, 0, 1] Energy:-32.0 Occurrence:7
[[1 0 0 2 1]
 [1 1 1 1 1]
 [1 2 2 2 1]
 [3 3 3 3 3]
 [4 2 5 5 5]]
Sample:[0, 1, 0,...,0, 0, 1] Energy:-32.0 Occurrence:89
[[1 0 0 1 1]
 [1 1 1 1 1]
 [1 1 2 3 1]
 [3 3 3 3 3]
 [4 3 5 5 5]]
Sample:[0, 1, 0,...,0, 0, 1] Energy:-32.0 Occurrence:3
[[1 0 0 1 1]
 [2 2 1 2 2]
 [1 1 2 3 1]
 [2 2 3 2 2]
 [4 3 5 5 5]]

#隣り合うマスは同時に１になりやすい（一塊のブロックにするため）
for k in range(dot_num):
    for i in range(h - 1):
        for j in range(w):
            command = f'H += -0.1 * (q{i}{j}{k} * q{i+1}{j}{k})'
            print(command)
            exec(command)
    for j in range(w - 1):
        for i in range(h):
            command = f'H += -0.1 * (q{i}{j}{k} * q{i}{j+1}{k})'
            print(command)
            exec(command)

response
Sample:[0, 1, 0,...,0, 0, 1] Energy:-34.10000000000001 Occurrence:234
[[1 0 0 1 1]
 [1 1 1 1 1]
 [1 1 2 3 1]
 [3 3 3 3 3]
 [4 3 5 5 5]]