反復コード

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

sns.set()
plt.rc("xtick", labelsize=20)
plt.rc("ytick", labelsize=20)
plt.rc("lines", linewidth=3)
plt.rc("font", size=20)
plt.rc("legend", fontsize="large")
plt.rc("axes", labelsize=20)
plt.rcParams["figure.figsize"] = 15, 6
plt.rcParams["legend.title_fontsize"] = 25

# Qiskit imports
from qiskit import (
    QuantumCircuit,
    QuantumRegister,
    ClassicalRegister,
)

# qiskit-ibm-runtime
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler
from qiskit_ibm_runtime.circuit import MidCircuitMeasure

service = QiskitRuntimeService()

# You can see all backends that support mid-circuit measurements by running the following code.
service.backends(filters=lambda b: "measure_2" in b.supported_instructions)

# Choose the least busy backend that supports mid-circuit measurements (`measure_2`).

backend = service.least_busy(
    filters=lambda b: "measure_2" in b.supported_instructions,
    operational=True,
    simulator=False,
    dynamic_circuits=True,
)
# backend = service.backend(backend_name) # alternatively, you could choose a specific backend
print(backend.name)

def build_qc(
    qreg_data,
    qreg_syndrome,
    creg_data,
    creg_syndrome,
    qreg_ref=None,
    creg_ref=None,
):
    """Build a typical error correction circuit"""
    if qreg_ref:
        return QuantumCircuit(
            qreg_data,
            qreg_syndrome,
            creg_data,
            creg_syndrome,
            qreg_ref,
            creg_ref,
        )
    else:
        return QuantumCircuit(
            qreg_data, qreg_syndrome, creg_data, creg_syndrome
        )


def encode_bit_flip(circuit, qreg_data, qreg_ref=None) -> QuantumCircuit:
    """Encode bit-flip. This is done by simply adding a cx"""

    for q in qreg_data:
        circuit.x(q)

    if qreg_ref:
        for q in qreg_ref:
            circuit.x(q)

    circuit.barrier()

    return circuit


def measure_syndrome_bit(
    circuit, qreg_data, qreg_syndrome, creg_syndrome, qreg_ref=None
):
    """
    Measure the syndrome by measuring the parity.
    We reset our ancilla qubits after measuring the stabilizer
    so we can reuse them for repeated stabilizer measurements.
    Because we have already observed the state of the qubit,
    we can write the conditional reset protocol directly to
    avoid another round of qubit measurement if we used
    the `reset` instruction.
    """
    circuit.cx(qreg_data[0], qreg_syndrome[0])
    circuit.cx(qreg_data[1], qreg_syndrome[0])
    circuit.cx(qreg_data[0], qreg_syndrome[1])
    circuit.cx(qreg_data[2], qreg_syndrome[1])
    circuit.barrier()

    for q_measure, c_measure in zip(qreg_syndrome, creg_syndrome):
        circuit.append(MidCircuitMeasure(), [q_measure], [c_measure])

    with circuit.if_test((creg_syndrome[0], 1)):
        circuit.x(qreg_syndrome[0])
    with circuit.if_test((creg_syndrome[1], 1)):
        circuit.x(qreg_syndrome[1])

    circuit.barrier()
    return circuit


def apply_correction_bit(circuit, qreg_data, creg_syndrome):
    """We can detect where an error occurred and correct our state"""
    with circuit.if_test((creg_syndrome, 3)):
        circuit.x(qreg_data[0])
    with circuit.if_test((creg_syndrome, 1)):
        circuit.x(qreg_data[1])
    with circuit.if_test((creg_syndrome, 2)):
        circuit.x(qreg_data[2])
    circuit.barrier()
    return circuit


def apply_final_readout(
    circuit, qreg_data, creg_data, qreg_ref=None, creg_ref=None
):
    """Read out the final measurements"""

    circuit.barrier()
    if qreg_ref:
        circuit.measure(qreg_ref, creg_ref)
    circuit.measure(qreg_data, creg_data)

    return circuit

def build_error_correction_sequence(
    num_cycles, cycles_per_circuit, nq_ref=3, delay=None
) -> QuantumCircuit:
    circuits = []
    reps = []

    qreg_data = QuantumRegister(3, name="qdata")
    qreg_syndrome = QuantumRegister(2, name="qsyndrome")
    creg_data = ClassicalRegister(3, name="cdata")
    creg_syndrome = ClassicalRegister(2, name="csyndrome")
    qreg_ref = QuantumRegister(nq_ref, name="qreference")
    creg_ref = ClassicalRegister(nq_ref, name="creference")

    circuit = build_qc(
        qreg_data,
        qreg_syndrome,
        creg_data,
        creg_syndrome,
        qreg_ref=qreg_ref,
        creg_ref=creg_ref,
    )
    circuit = encode_bit_flip(circuit, qreg_data, qreg_ref=qreg_ref)

    circuit_n = circuit.copy()
    circuit_n = apply_final_readout(
        circuit_n, qreg_data, creg_data, qreg_ref=qreg_ref, creg_ref=creg_ref
    )
    circuits.append(circuit_n)
    reps.append(0)

    for i in range(1, num_cycles + 1):
        if delay:
            circuit.delay(delay, unit="us")
            circuit.barrier()
        circuit = measure_syndrome_bit(
            circuit,
            qreg_data,
            qreg_syndrome,
            creg_syndrome,
            qreg_ref=qreg_ref,
        )
        circuit = apply_correction_bit(circuit, qreg_data, creg_syndrome)
        circuit_n = circuit.copy()
        if i % cycles_per_circuit == 0:
            circuit_n = apply_final_readout(
                circuit_n,
                qreg_data,
                creg_data,
                qreg_ref=qreg_ref,
                creg_ref=creg_ref,
            )
            circuits.append(circuit_n)
            reps.append(i)

    return circuits, np.array(reps)


num_cycles = 40
cycles_per_circuit = 4
nq_ref = 3
circuits, rep_array = build_error_correction_sequence(
    num_cycles, cycles_per_circuit, nq_ref=3, delay=5
)
circuits[1].draw(output="mpl", cregbundle=False, fold=50)

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_circuit_ref = pm.run(circuits[-1])

init_layout = isa_circuit_ref.layout.initial_index_layout(
    filter_ancillas=True
)
print(init_layout)

# get all qubits ordered by T1
t1_data = []
for i in range(backend.num_qubits):
    try:
        t1_us = backend.properties().t1(i) * 1e6
    except Exception:
        t1_us = 0.0
    t1_data.append((i, t1_us))

t1_data_sorted = sorted(t1_data, key=lambda x: x[1], reverse=True)

# exclude the qubits we have already mapped the error correcting code to
t1_data_sorted = [
    t1_data for t1_data in t1_data_sorted if t1_data[0] not in init_layout[:5]
]

# use the best qubits in terms of T1 for the reference qubits
init_layout = init_layout[:5] + [t1_data[0] for t1_data in t1_data_sorted[:3]]
print(init_layout)

# These are the resulting T1 times
properties = backend.properties()
print("Amplitude damping decoherence times for code data qubits:")
for q in init_layout[:3]:
    t1 = properties.t1(q)
    print(f"qubit {q}: T1 = {t1 * 1e6:.0f} mus")

print("\nAmplitude damping decoherence times for reference qubits:")
for q in init_layout[-3:]:
    t1 = properties.t1(q)
    print(f"qubit {q}: T1 = {t1 * 1e6:.0f} mus")

# now we transpile all circuits to this initial layout; this way each circuit is run on the same qubits and we can make a fair comparison
pm = generate_preset_pass_manager(
    backend=backend,
    optimization_level=3,
    initial_layout=init_layout,
)
isa_circuits = pm.run(circuits)

isa_circuits[1].draw("mpl", cregbundle=False, fold=50)

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-REPCODE"]
sampler.options.max_execution_time = 600  # this workload is expected to be under 10s, but it is generally a good habit to set a max execution time (here 600s = 10m)

job = sampler.run(isa_circuits, shots=1000)
print(job.job_id())

job.status()

results = job.result()

def correct_counts(counts_dict):
    """
    Corrects the measured logical qubit encoded in the repetition code using majority vote
    """

    result = {"000": 0, "111": 0}
    for bitstring, count in counts_dict.items():
        key = "111" if bitstring.count("1") > 1 else "000"
        result[key] += count
    return result


accuracy = []  # logical qubit
accuracies_ref = np.zeros(
    (len(results), nq_ref)
)  # accuracies on individual reference qubits

for n, pub_result in enumerate(results):
    # logical accuracy (one minus error probability) for active error correction with repetition code
    counts = pub_result.data.cdata.get_counts()
    shots = sum(counts.values())
    counts_corrected = correct_counts(counts)
    accuracy.append(counts_corrected.get("111", 0) / shots)

    # accuracy for individual physical reference qubits without any error correction
    for i in range(nq_ref):
        counts = pub_result.data.creference.slice_bits(i).get_counts()
        accuracies_ref[n, i] = counts.get("1", 0) / shots

accuracy = np.array(accuracy)

def error_proba_t1(N, t1):
    """
    Exponential fitting function for amplitude damping vs. number of cycles
    """

    t_cycle = 7.3e-6  # approximate time per cycle = 5 mus delay + 2.3 mus for error correction
    return 1 - np.exp(-t_cycle * N / t1)


fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(
    rep_array, (1.0 - accuracy) * 100.0, "ko-", linewidth=3, label="rep code"
)

for i in range(nq_ref):
    accuracy_1q = accuracies_ref[:, i]
    if i == 0:
        ax.plot(
            rep_array,
            (1.0 - accuracy_1q) * 100.0,
            "go-",
            linewidth=1,
            label="1q reference",
        )
    else:
        ax.plot(rep_array, (1.0 - accuracy_1q) * 100.0, "go-", linewidth=1)

params_bf, pcov = curve_fit(
    error_proba_t1, rep_array, 1.0 - accuracy, bounds=([0, 5e-3])
)
t1_bf = params_bf[0]
print(f"Best-fit effective T1 = {t1_bf * 1e6:.0f} us")
error_prob_bf = error_proba_t1(
    rep_array, t1_bf
)  # np.array([1 - np.exp(-t_cycle * nt/t1) for nt in rep_array])
ax.plot(rep_array, error_prob_bf * 100.0, "k--", linewidth=1)

ax.set_xlabel("error correction rounds")
ax.set_ylabel("error [%]")
ax.set_ylim(bottom=0)
ax.legend(fontsize=15);