절단된 벨 쌍을 이용한 동적 회로 벤치마킹

from qiskit import QuantumCircuit

from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler
from qiskit.transpiler import generate_preset_pass_manager

import numpy as np
import matplotlib.pyplot as plt


def create_bell_stab(initial_layouts):
    """
    Create a circuit for a 1D chain of qubits (number
    of qubits must be a multiple of 4), where a middle
    Bell pair is consumed to create a Bell at the edge.
    Takes as input a list of lists, where each element
    of the list is a 1D chain of physical qubits that is
    used as the initial_layout for the transpiled circuit.
    Returns a list of length-2 tuples, each tuple
    contains a circuit to measure the ZZ stabilizer and
    a circuit to measure the XX stabilizer of the edge
    Bell state.
    """
    bell_circuits = []
    for (
        initial_layout
    ) in initial_layouts:  # Iterate over chains of physical qubits
        assert (
            len(initial_layout) % 4 == 0
        ), "The length of the chain must be a multiple of 4, "
        f"len(inital_layout)={len(initial_layout)}"
        num_pairs = len(initial_layout) // 4

        bell_parallel = QuantumCircuit(4 * num_pairs, 4 * num_pairs)

        for pair_idx in range(num_pairs):
            (q0, q1, q2, q3) = (
                pair_idx * 4,
                pair_idx * 4 + 1,
                pair_idx * 4 + 2,
                pair_idx * 4 + 3,
            )
            (c0, c1) = pair_idx * 4, pair_idx * 4 + 3  # edge qubits
            (ca0, ca1) = pair_idx * 4 + 1, pair_idx * 4 + 2  # middle qubits

            bell_parallel.h(q0)
            bell_parallel.h(q1)
            bell_parallel.cx(q1, q2)
            bell_parallel.cx(q0, q1)
            bell_parallel.cx(q2, q3)
            bell_parallel.h(q2)

        # add barrier BEFORE measurements and add id in conditional
        bell_parallel.barrier()
        for pair_idx in range(num_pairs):
            (q0, q1, q2, q3) = (
                pair_idx * 4,
                pair_idx * 4 + 1,
                pair_idx * 4 + 2,
                pair_idx * 4 + 3,
            )
            (ca0, ca1) = pair_idx * 4 + 1, pair_idx * 4 + 2  # middle qubits

            bell_parallel.measure(q1, ca0)
            bell_parallel.measure(q2, ca1)
        # bell_parallel.barrier() #remove barrier after measurement

        for pair_idx in range(num_pairs):
            (q0, q1, q2, q3) = (
                pair_idx * 4,
                pair_idx * 4 + 1,
                pair_idx * 4 + 2,
                pair_idx * 4 + 3,
            )
            (ca0, ca1) = pair_idx * 4 + 1, pair_idx * 4 + 2  # middle qubits
            with bell_parallel.if_test((ca0, 1)):
                bell_parallel.x(q3)
            with bell_parallel.if_test((ca1, 1)):
                bell_parallel.z(q0)
                bell_parallel.id(q0)  # add id here for correct alignment

        bell_zz = bell_parallel.copy()
        bell_zz.barrier()
        bell_xx = bell_parallel.copy()
        bell_xx.barrier()
        for pair_idx in range(num_pairs):
            (q0, q1, q2, q3) = (
                pair_idx * 4,
                pair_idx * 4 + 1,
                pair_idx * 4 + 2,
                pair_idx * 4 + 3,
            )
            bell_xx.h(q0)
            bell_xx.h(q3)
        bell_xx.barrier()
        for pair_idx in range(num_pairs):
            (q0, q1, q2, q3) = (
                pair_idx * 4,
                pair_idx * 4 + 1,
                pair_idx * 4 + 2,
                pair_idx * 4 + 3,
            )
            (c0, c1) = pair_idx * 4, pair_idx * 4 + 3  # edge qubits

            bell_zz.measure(q0, c0)
            bell_zz.measure(q3, c1)

            bell_xx.measure(q0, c0)
            bell_xx.measure(q3, c1)

        bell_circuits.append(bell_zz)
        bell_circuits.append(bell_xx)

    return bell_circuits


def get_mse(result, initial_layouts):
    """
    given a result object and the initial layouts,
    returns a dict of layouts and their mse
    """
    layout_mse = {}
    for layout_idx, initial_layout in enumerate(initial_layouts):
        layout_mse[tuple(initial_layout)] = {}

        num_pairs = len(initial_layout) // 4

        counts_zz = result[2 * layout_idx].data.c.get_counts()
        total_shots = sum(counts_zz.values())

        # Get ZZ expectation value
        exp_zz_list = []
        for pair_idx in range(num_pairs):
            exp_zz = 0
            for bitstr, shots in counts_zz.items():
                bitstr = bitstr[::-1]  # reverse order to big endian
                b1, b0 = (
                    bitstr[pair_idx * 4],
                    bitstr[pair_idx * 4 + 3],
                )  # parse bitstring to get edge measurements for each 4-q chain
                z_val0 = 1 if b0 == "0" else -1
                z_val1 = 1 if b1 == "0" else -1
                exp_zz += z_val0 * z_val1 * shots
            exp_zz /= total_shots
            exp_zz_list.append(exp_zz)

        counts_xx = result[2 * layout_idx + 1].data.c.get_counts()
        total_shots = sum(counts_xx.values())

        # Get XX expectation value
        exp_xx_list = []
        for pair_idx in range(num_pairs):
            exp_xx = 0
            for bitstr, shots in counts_xx.items():
                bitstr = bitstr[::-1]  # reverse order to big endian
                b1, b0 = (
                    bitstr[pair_idx * 4],
                    bitstr[pair_idx * 4 + 3],
                )  # parse bitstring to get edge measurements for each 4-q chain
                x_val0 = 1 if b0 == "0" else -1
                x_val1 = 1 if b1 == "0" else -1
                exp_xx += x_val0 * x_val1 * shots
            exp_xx /= total_shots
            exp_xx_list.append(exp_xx)

        mse_list = [
            ((exp_zz - 1) ** 2 + (exp_xx - 1) ** 2) / 2
            for exp_zz, exp_xx in zip(exp_zz_list, exp_xx_list)
        ]

        print(f"layout {initial_layout}")
        for idx in range(num_pairs):
            layout_mse[tuple(initial_layout)][
                tuple(initial_layout[4 * idx : 4 * idx + 4])
            ] = mse_list[idx]
            print(
                f"qubits: {initial_layout[4*idx:4*idx+4]}, mse:, "
                f"{round(mse_list[idx],4)}"
            )
            # print(f'exp_zz: {round(exp_zz_list[idx],4)},
            # exp_xx: {round(exp_xx_list[idx],4)}')
        print(" ")
    return layout_mse


def plot_mse_ecdfs(layouts_mse, combine_layouts=False):
    """
    Plot CDF of MSE data for multiple layouts.
    Optionally combine all data in a single CDF
    """

    if not combine_layouts:
        for initial_layout, layouts in layouts_mse.items():
            sorted_layouts = dict(
                sorted(layouts.items(), key=lambda item: item[1])
            )  # sort layouts by mse

            # get layouts and mses
            layout_list = list(sorted_layouts.keys())
            mse_list = np.asarray(list(sorted_layouts.values()))

            # convert to numpy
            x = np.array(mse_list)
            y = np.arange(1, len(x) + 1) / len(x)

            # Prepend (x[0], 0) to start CDF at zero
            x = np.insert(x, 0, x[0])
            y = np.insert(y, 0, 0)

            # Create the plot
            plt.plot(
                x,
                y,
                marker="x",
                linestyle="-",
                label=f"qubits: {initial_layout}",
            )

            # add qubits labels for the edge pairs
            for xi, yi, q in zip(x[1:], y[1:], layout_list):
                plt.annotate(
                    [q[0], q[3]],
                    (xi, yi),
                    textcoords="offset points",
                    xytext=(5, -10),
                    ha="left",
                    fontsize=8,
                )

    elif combine_layouts:
        all_layouts = {}
        all_initial_layout = []
        for (
            initial_layout,
            layouts,
        ) in layouts_mse.items():  # puts together all layout information
            all_layouts.update(layouts)
            all_initial_layout += initial_layout

        sorted_layouts = dict(
            sorted(all_layouts.items(), key=lambda item: item[1])
        )  # sort layouts by mse

        # get layouts and mses
        layout_list = list(sorted_layouts.keys())
        mse_list = np.asarray(list(sorted_layouts.values()))

        # convert to numpy
        x = np.array(mse_list)
        y = np.arange(1, len(x) + 1) / len(x)

        # Prepend (x[0], 0) to start CDF at zero
        x = np.insert(x, 0, x[0])
        y = np.insert(y, 0, 0)

        # Create the plot
        plt.plot(
            x,
            y,
            marker="x",
            linestyle="-",
            label=f"qubits: {sorted(list(set(all_initial_layout)))}",
        )

        # add qubit labels for the edge pairs
        for xi, yi, q in zip(x[1:], y[1:], layout_list):
            plt.annotate(
                [q[0], q[3]],
                (xi, yi),
                textcoords="offset points",
                xytext=(5, -10),
                ha="left",
                fontsize=8,
            )

    plt.xscale("log")
    plt.xlabel("Mean squared error of ⟨ZZ⟩ and ⟨XX⟩")
    plt.ylabel("Cumulative distribution function")
    plt.title("CDF for different initial layouts")
    plt.grid(alpha=0.3)
    plt.show()

from qiskit_aer import AerSimulator

# 4-qubit chain for simulation
sim_layout = [[0, 1, 2, 3]]

aer_backend = AerSimulator()
sim_circuits = create_bell_stab(sim_layout)
sim_circuits[1].draw("mpl", fold=-1, idle_wires=False)

pm_sim = generate_preset_pass_manager(
    optimization_level=0, backend=aer_backend, initial_layout=sim_layout[0]
)
isa_sim_circuits = pm_sim.run(sim_circuits)
isa_sim_circuits[1].draw("mpl", fold=-1, idle_wires=False)

# Run on noiseless simulator
sampler_sim = Sampler(mode=aer_backend)
sim_job = sampler_sim.run(isa_sim_circuits)
sim_mse = get_mse(sim_job.result(), sim_layout)

service = QiskitRuntimeService()
backend = service.least_busy(operational=True)

from itertools import chain
from collections import defaultdict


def stripes16_from_backend(backend):
    """
    Creates stripes of 16 qubits, four non-overlapping
    four-qubit chains, that cover as much of the coupling
    map as possible. Returns any unused qubits as leftovers.
    """
    # get the undirected adjacency list
    edges = backend.coupling_map.get_edges()
    graph = defaultdict(set)
    for u, v in edges:
        graph[u].add(v)
        graph[v].add(u)

    qubits = sorted(graph)  # all qubit indices that appear

    # greedy search for 4-long linear chains (blocks) ────────────
    used = set()  # qubits already placed in a block
    blocks = []  # each block is a four-qubit list

    for q in qubits:  # deterministic order for reproducibility
        if q in used:
            continue  # already consumed by earlier block

        # depth-first "straight" walk of length 3 without revisiting nodes
        def extend(path):
            if len(path) == 4:
                return path
            tip = path[-1]
            for nbr in sorted(graph[tip]):  # deterministic
                if nbr not in path and nbr not in used:
                    maybe = extend(path + [nbr])
                    if maybe:
                        return maybe
            return None

        block = extend([q])
        if block:  # found a 4-node path
            blocks.append(block)
            used.update(block)

    # bundle four four-qubit blocks into one 16-qubit
    # stripe (max number of measurement compatible with if-else)
    stripes = [
        list(chain.from_iterable(blocks[i : i + 4]))
        for i in range(0, len(blocks) // 4 * 4, 4)  # full groups of four
    ]

    leftovers = set(qubits) - set(chain.from_iterable(stripes))
    return stripes, leftovers

initial_layouts, leftover = stripes16_from_backend(backend)

# -------------------------Step 1-------------------------
circuits = create_bell_stab(initial_layouts)

# -------------------------Step 2-------------------------
isa_circuits = []
for ind, init_layout in enumerate(initial_layouts):
    pm = generate_preset_pass_manager(
        optimization_level=0, backend=backend, initial_layout=init_layout
    )
    isa_circ = pm.run(circuits[ind * 2 : ind * 2 + 2])
    isa_circuits.extend(isa_circ)
isa_circuits[1].draw("mpl", fold=-1, idle_wires=False)

# -------------------------Step 3-------------------------
sampler = Sampler(mode=backend)
sampler.options.environment.job_tags = ["TUT_BDC"]
job = sampler.run(isa_circuits)

# -------------------------Step 4-------------------------
layouts_mse = get_mse(job.result(), initial_layouts)

plot_mse_ecdfs(layouts_mse, combine_layouts=True)