import numpy as np


# 確率 p でエラーを発生させたサイズ d のデータ量子ビット格子を生成
# 確率 p で 1（エラーあり）、確率 1-p で 0（エラーなし）が格納される
def generate_qubit_lattice(d: int, p: float, seed: int | None) -> np.ndarray:
    rng = np.random.default_rng(seed)
    qubits = np.where(rng.random((d * 2 - 1, d)) < p, 1, 0)
    qubits[1::2, -1] = (
        0  # 1,3,5,... 行目はデータ量子ビットが 1 つ少ないので、右端を 0 で埋める
    )
    return qubits


# データ量子ビットのエラー情報から補助量子ビットでの測定結果をシミュレート
def measure_stabilizers(qubits: np.ndarray) -> np.ndarray:
    d = len(qubits[0])
    stabilizers = np.zeros((d - 1, d), dtype=int)
    for i in range(d - 1):
        for j in range(d):
            stabilizers[i, j] = (
                qubits[i * 2, j]
                ^ qubits[i * 2 + 1, j - 1]
                ^ qubits[i * 2 + 1, j]
                ^ qubits[i * 2 + 2, j]
            )
    return stabilizers


d = 11  # 格子サイズ
p = 0.08  # データ量子ビットのエラー率
qubits = generate_qubit_lattice(
    d, p, 0
)  # 確率 p でエラーを発生させたサイズ d のデータ量子ビット格子
stabilizers = measure_stabilizers(qubits)  # 補助量子ビットでの測定結果

import cv2
from IPython.display import display, Image


# 量子ビット格子の様子を図示
def show_sc_img(
    qubits: np.ndarray,  # データ量子ビットのエラー情報
    modified: np.ndarray | None,  # 修正したデータ量子ビットの情報
    stabilizers: np.ndarray,  # 補助量子ビットの情報
    radius: int = 10,
    thickness: int = 2,
) -> None:
    d = len(qubits[0])
    size = 6 * radius * d
    img = np.full((size, size, 3), 255, dtype=np.uint8)

    # データ量子ビット（丸、辺）を描く
    for i in range(d * 2 - 1):
        for j in range(d - i % 2):
            center = ((i % 2 + j * 2 + 1) * 3 * radius, (i + 1) * 3 * radius)

            # データ量子ビットが反転している場合、丸と辺を赤で塗りつぶす
            if qubits[i, j]:
                # 辺の伸び方（0, 2, ... 行目は縦、1, 3, ... 行目は横）
                dx = 3 * radius if i % 2 == 1 else 0
                dy = 3 * radius if i % 2 == 0 else 0
                cv2.line(
                    img=img,
                    pt1=(center[0] - dx, center[1] - dy),
                    pt2=(center[0] + dx, center[1] + dy),
                    color=(0, 0, 255),
                    thickness=thickness * 2,
                )
                cv2.circle(
                    img,
                    center=center,
                    radius=radius,
                    color=(0, 0, 255),
                    thickness=-1,
                )

            # 修正された箇所は緑色、そうでなければ黒色で丸の輪郭を描く
            outline_color = (
                (0, 255, 0) if modified is not None and modified[i, j] else (0, 0, 0)
            )
            cv2.circle(
                img=img,
                center=center,
                radius=radius,
                color=outline_color,
                thickness=thickness,
            )

    # 補助量子ビット（四角）を描く
    for i in range(d - 1):
        for j in range(d):
            pt1 = ((j * 6 + 2) * radius, (i * 6 + 5) * radius)
            pt2 = ((j * 6 + 4) * radius, (i * 6 + 7) * radius)

            # 補助量子ビットまわりで反転している量子ビットの数が奇数の場合、黒で塗りつぶす
            if stabilizers[i, j]:
                cv2.rectangle(img, pt1=pt1, pt2=pt2, color=(0, 0, 0), thickness=-1)

            # 黒色で四角の輪郭を描く
            cv2.rectangle(
                img=img, pt1=pt1, pt2=pt2, color=(0, 0, 0), thickness=thickness
            )

    # 描画したものを表示
    img = cv2.resize(img, (512, 512))
    _, buf = cv2.imencode(".png", img)
    display(Image(data=buf.tobytes()))

show_sc_img(qubits, None, stabilizers)

from amplify import VariableGenerator, equal_to, PolyArray, Model
from amplify import sum as amplify_sum


# 補助量子ビットの測定結果から組み合わせ最適化モデルを構築
def construct_model(d: int, stabilizers: np.ndarray) -> tuple[PolyArray, Model]:
    m = stabilizers  # 補助量子ビットでの測定結果

    gen = VariableGenerator()

    # データ量子ビットのエラー有無を表すバイナリ変数
    x = gen.array("Binary", d * 2 - 1, d)
    # 1,3,5,... 行目はデータ量子ビットが 1 つ少ないので、右端を定数 0 で埋める
    x[1::2, -1] = 0

    # 補助バイナリ変数
    z = gen.array("Binary", d - 1, d)

    # sum メソッドで目的関数を構築
    objective = x.sum()

    # 隣接するデータ量子ビットのエラー個数
    S = [
        [
            x[i * 2, j] + x[i * 2 + 1, j - 1] + x[i * 2 + 1, j] + x[i * 2 + 2, j]
            for j in range(d)
        ]
        for i in range(d - 1)
    ]

    # 制約条件の左辺
    left_values = [
        S[i][j] - 2 * z[i, j] - m[i, j] for i in range(d - 1) for j in range(d)
    ]

    # equal_to 関数と sum 関数で制約条件を構築
    constraints = amplify_sum([equal_to(left_value, 0) for left_value in left_values])

    # 目的関数と制約条件からモデルを構築
    model = Model(objective, constraints)

    return (x, model)


x, model = construct_model(d, stabilizers)

from amplify import FixstarsClient, solve, Result
from datetime import timedelta


def solve_optimization(model: Model, timeout_ms: int = 2000) -> Result:
    # クライアントを設定
    client = FixstarsClient()
    # client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # ローカル環境等で使用する場合は、Amplify AEのアクセストークンを入力してください。
    client.parameters.timeout = timedelta(
        milliseconds=timeout_ms
    )  # 実行時間を設定 (デフォルトは 2 秒)

    # 問題を解く
    result = solve(model, client)

    # result が空の場合、制約条件を満たす解が得られなかったことを意味する
    if len(result) == 0:
        raise RuntimeError("制約条件を満たす解が見つかりませんでした")

    return result


result = solve_optimization(model)

solution = x.evaluate(result.best.values).astype(int)  # 求めた解を元の変数に代入

show_sc_img(qubits, solution, stabilizers)

modified_qubits = qubits ^ solution  # エラー訂正後のエラー状況
modified_stabilizers = measure_stabilizers(
    modified_qubits
)  # エラー訂正後の補助量子ビットの測定結果

show_sc_img(modified_qubits, None, modified_stabilizers)

def check_success_or_failure(
    modified_qubits: np.ndarray, modified_stabilizers: np.ndarray
) -> None:
    if modified_stabilizers.sum() != 0:
        print(
            "復号失敗: 隣接するデータ量子ビットのエラー個数が奇数の補助量子ビットがあります"
        )
    elif modified_qubits[0].sum() % 2 == 1:
        print("復号失敗: 一番上の行にあるデータ量子ビットのエラー個数が奇数です")
    else:
        print("復号成功です！")


check_success_or_failure(modified_qubits, modified_stabilizers)

# 一連の処理を実行する
def simulate(d: int, p: float, seed: int | None, timeout_ms: int) -> None:
    # データ量子ビット格子を生成
    qubits = generate_qubit_lattice(d, p, seed)

    # 補助量子ビットでの測定結果
    stabilizers = measure_stabilizers(qubits)

    # 生成されたエラー状況を表示
    show_sc_img(qubits, None, stabilizers)

    # 組み合わせ最適化モデルを構築
    x, model = construct_model(d, stabilizers)

    # 最適化の実行
    result = solve_optimization(model, timeout_ms)

    # 推定エラー位置を表示
    solution = x.evaluate(result.best.values).astype(int)
    show_sc_img(qubits, solution, stabilizers)

    # 訂正後のエラー状況を表示
    modified_qubits = qubits ^ solution
    modified_stabilizers = measure_stabilizers(modified_qubits)
    show_sc_img(modified_qubits, None, modified_stabilizers)

    # 復号の成否を判定
    check_success_or_failure(modified_qubits, modified_stabilizers)

d = 15  # 格子サイズ
p = 0.08  # エラー発生確率
seed = None  # 乱数のシード値 (整数を指定: 問題固定、None: ランダム生成)
timeout_ms = 3000  # 計算時間 (単位はミリ秒)
simulate(d, p, seed, timeout_ms)

イジングマシンによる量子誤り訂正の解法¶

Surface Code とは¶

エラーの訂正方法¶

Surface Code 復号問題の定式化¶

最適化問題の入力および決定変数の定義¶

目的関数の定式化¶

制約条件の定式化¶

Note¶

Surface Code 復号問題の求解¶

データ量子ビット・補助量子ビットの添え字の割り当て¶

問題の生成¶

Note¶

組合せ最適化モデルの構築¶

最適化の実行¶

結果の確認¶

〈補足〉「どの補助量子ビット $j$ についても $S_j \le \underline{3}$ を満たす最適解が必ず存在する」ことの説明¶