Long-range entanglement with dynamic circuits

import random
from typing import List, Optional
 
import matplotlib.pyplot as plt
import numpy as np
 
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit import IfElseOp
from qiskit.circuit.classical import expr
from qiskit.primitives import BitArray
from qiskit.transpiler import CouplingMap
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime.transpiler.passes.scheduling import (
    DynamicCircuitInstructionDurations,
)

def get_1d_qubit_line(
    coupling_map: CouplingMap, start_qubit: int = 0
) -> List[int]:
    """
    Use random search to find a 1D line, repeating until we find a line of sufficient length.
    This is very inefficient, but it's fine for this demo.
    """
 
    def get_random_line(
        coupling_map: CouplingMap, start_qubit: int
    ) -> List[int]:
        """
        Do a random, self-avoiding walk on the coupling map to get a 1D line.
        """
        edge_list = list(coupling_map.get_edges())
        path = [start_qubit]
        while True:
            # Get edges connected to the current qubit that don't connect to a qubit we've already visited
            current_qubit = path[-1]
            previously_visited_qubits = path[:-1]
            candidate_edges = [
                edge
                for edge in edge_list
                if current_qubit in edge
                and not any(
                    qubit in previously_visited_qubits for qubit in edge
                )
            ]
            if candidate_edges == []:
                return path
            next_edge = random.choice(candidate_edges)
            edge_list.remove(next_edge)
            next_qubit = next(
                qubit for qubit in next_edge if qubit != path[-1]
            )
            path.append(next_qubit)
 
    # Now repeat the random walk many times to find a long line
    return max(
        [get_random_line(coupling_map, start_qubit) for _ in range(1000)],
        key=lambda x: len(x),
    )

def coupling_map_from_qubit_line(
    coupling_map: CouplingMap, qubit_line: List[List[int]]
) -> CouplingMap:
    """
    Modify the full coupling map to force linearity in the qubit layout
    """
    new_coupling_map = []
    # Iterate the line pair-wise and append edges that contain both qubits
    for qubit, next_qubit in zip(qubit_line, qubit_line[1:]):
        edge = next(
            edge
            for edge in coupling_map.get_edges()
            if qubit in edge and next_qubit in edge
        )
        new_coupling_map.append(edge)
    return CouplingMap(new_coupling_map)
 
 
# Set up access to IBM Quantum devices
service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, dynamic_circuits=True
)
if "if_else" not in backend.target.operation_names:
    backend.target.add_instruction(IfElseOp, name="if_else")
 
# Set qubit line and coupling map
QUBIT_LINE = get_1d_qubit_line(backend.coupling_map)
COUPLING_MAP_1D = coupling_map_from_qubit_line(
    backend.coupling_map, QUBIT_LINE
)
MAX_POSSIBLE_QUBITS_BTW_CNOT = len(QUBIT_LINE) - 2
 
# Use this duration class to get appropriate durations for dynamic
# circuit backend scheduling
DURATIONS = DynamicCircuitInstructionDurations.from_backend(backend)
 
print(f"Backend is: {backend.name}")
print(f"Found line of length {len(QUBIT_LINE)}.")
print("Using coupling map:\n", COUPLING_MAP_1D)
print(
    f"Maximum number of qubits between CNOT for {backend.name} is {MAX_POSSIBLE_QUBITS_BTW_CNOT} with the given qubit line."
)

# Level of optimizations the transpiler uses: There are 4 optimization levels ranging from 0 to 3,
# where higher optimization levels take more time and computational effort but may yield a more
# optimal circuit. level 0 does no explicit optimization, it will just try to make a circuit
# runnable by transforming it to match a topology and basis gate set, if necessary.
# Level 1, 2 and 3 do light, medium and heavy optimization, using a combination of passes, and by
# configuring the passes to search for better solutions.
OPTIMIZATION_LEVEL = 0
 
# Set the number of repetitions of each circuit, for sampling.
# The number of qubits between the control and target are grouped into blocks
# of length 4. The provided min and max number of qubits will be modified to
# align with these block sizes.
SHOTS = 1000
 
# Set the number of times the experimental counts are resampled.
# This is used to calculate mean and standard deviations of results.
BOOTSTRAP_SAMPLES = 10
 
# The min number of qubits between the control and target qubits on line
MIN_NUMBER_QUBITS = 0
 
# The max number of qubits between the control and target qubits on line
# The max for MIN_NUMBER_QUBITS is len(QUBIT_LINE) - 2
# max_number_qubits must satisfy MAX_NUMBER_QUBITS - MIN_NUMBER_QUBITS = 3 (mod 4)
# at this point. This is just to make things easier to break jobs up. Not a real limitation.
MAX_NUMBER_QUBITS = 15
assert (
    (MAX_NUMBER_QUBITS - MIN_NUMBER_QUBITS) % 4 == 3
), "MAX_NUMBER_QUBITS must satisfy MAX_NUMBER_QUBITS - MIN_NUMBER_QUBITS = 3 (mod 4)"
assert (MAX_NUMBER_QUBITS + 2) <= len(
    QUBIT_LINE
), "MAX_NUMBER_QUBITS must satisfy MAX_NUMBER_QUBITS + 2 <= len(QUBIT_LINE)"

def CNOT_unitary(
    qc: QuantumCircuit, control_qubit: int, target_qubit: int
) -> QuantumCircuit:
    """Generate a CNOT gate between data qubit control_qubit and data qubit target_qubit using local CNOTs
 
    Assumes that the long-range CNOT gate will be spanning a 1D chain of n-qubits subject to nearest-neighbor
    connections only with the chain starting at the control qubit and finishing at the target qubit.
 
    Assumes that control_qubit < target_qubit (as integers) and that the provided circuit qc has |0> set
    qubits control_qubit+1, ..., target_qubit-1
 
    Args:
        qc (QuantumCicruit) : A Quantum Circuit to add the long range localized unitary CNOT
        control_qubit (int) : The qubit used as the control.
        target_qubi (int) : The qubit targeted by the gate.
 
    Example:
 
        qc = QuantumCircuit(8,2)
        qc = CNOT_unitary(qc, 0, 7)
 
    Returns:
        QuantumCircuit
 
    """
    assert target_qubit > control_qubit
    n = target_qubit - control_qubit - 1
    k = int(n / 2)
    qc.barrier()
    for i in range(control_qubit, control_qubit + k):
        qc.cx(i, i + 1)
        qc.cx(i + 1, i)
        qc.cx(-i - 1, -i - 2)
        qc.cx(-i - 2, -i - 1)
    if n % 2 == 1:
        qc.cx(k + 2, k + 1)
        qc.cx(k + 1, k + 2)
    qc.barrier()
    qc.cx(k, k + 1)
    for i in range(control_qubit, control_qubit + k):
        qc.cx(k - i, k - 1 - i)
        qc.cx(k - 1 - i, k - i)
        qc.cx(k + i + 1, k + i + 2)
        qc.cx(k + i + 2, k + i + 1)
    if n % 2 == 1:
        qc.cx(-2, -1)
        qc.cx(-1, -2)
    qc.barrier()
    return qc

def build_circuits_uni(n: int) -> List[QuantumCircuit]:
    """Builds the unitary circuits needed to estimate the average gate fidelity
 
    Args:
        n (int): Number of qubits between the control and target of the CNOT
    """
    assert n >= 0, "Error: n needs to be a non-negative integer"
    circuits = []
    for i in range(n):
        circuit = QuantumCircuit(i + 2, 2)
        circuit.h(0)
        circuit = CNOT_unitary(circuit, 0, i + 1)  # Add long range CNOT
        circuit.measure([0, i + 1], [0, 1])  # Prepare circuits to measure
        circuits.append(circuit)
    return circuits

circuits_uni = build_circuits_uni(MAX_NUMBER_QUBITS)

circuits_uni[-1].draw("mpl", fold=-1)

def CNOT_dyn(
    qc: QuantumCircuit,
    control_qubit: int,
    target_qubit: int,
    c1: Optional[ClassicalRegister] = None,
    c2: Optional[ClassicalRegister] = None,
    add_barriers: Optional[bool] = True,
) -> QuantumCircuit:
    """Generate a CNOT gate between data qubit control_qubit and data qubit target_qubit using Bell Pairs.
 
    Post processing is used to enable the CNOT gate via the provided classical registers c1 and c2
 
    Assumes that the long-range CNOT gate will be spanning a 1D chain of n-qubits subject to nearest-neighbor
    connections only with the chain starting at the control qubit and finishing at the target qubit.
 
    Assumes that control_qubit < target_qubit (as integers) and that the provided circuit qc has |0> set
    qubits control_qubit+1, ..., target_qubit-1
 
    n = target_qubit - control_qubit - 1 : Number of qubits between the target and control qubits
    k = int(n/2) : Number of Bell pairs created
 
    Args:
        qc (QuantumCicruit) : A Quantum Circuit to add the long range localized unitary CNOT
        control_qubit (int) : The qubit used as the control.
        target_qubi (int) : The qubit targeted by the gate.
 
    Optional Args:
        c1 (ClassicalRegister) : Default = None. Required if n > 1. Register requires k bits
        c2 (ClassicalRegister) : Default = None. Required if n > 0. Register requires n - k bits
        add_barriers (bool) : Default = True. Include barriers before and after long range CNOT
 
    Note: This approach uses two if_test statements. A better (more performant) approach is
    to have the parity values combined into a single classical register and then use a switch
    statement. This was done in the associated paper my modifying the qasm file directly. The ability
    to use a switch statement via Qiskit in this way is a future release capability.
 
    Returns:
        QuantumCircuit
    """
    assert target_qubit > control_qubit
    n = target_qubit - control_qubit - 1
    t = int(n / 2)
 
    if add_barriers is True:
        qc.barrier()
 
    # Determine where to start the bell pairs and
    # add an extra CNOT when n is odd
    if n % 2 == 0:
        x0 = 1
    else:
        x0 = 2
        qc.cx(0, 1)
 
    # Create t Bell pairs
    for i in range(t):
        qc.h(x0 + 2 * i)
        qc.cx(x0 + 2 * i, x0 + 2 * i + 1)
 
    # Entangle Bell pairs and data qubits and measure
    for i in range(t + 1):
        qc.cx(x0 - 1 + 2 * i, x0 + 2 * i)
 
    for i in range(1, t + x0):
        if i == 1:
            qc.h(2 * i + 1 - x0)
            qc.measure(2 * i + 1 - x0, c2[i - 1])
            parity_control = expr.lift(c2[i - 1])
        else:
            qc.h(2 * i + 1 - x0)
            qc.measure(2 * i + 1 - x0, c2[i - 1])
            parity_control = expr.bit_xor(c2[i - 1], parity_control)
 
    for i in range(t):
        if i == 0:
            qc.measure(2 * i + x0, c1[i])
            parity_target = expr.lift(c1[i])
        else:
            qc.measure(2 * i + x0, c1[i])
            parity_target = expr.bit_xor(c1[i], parity_target)
 
    if n > 0:
        with qc.if_test(parity_control):
            qc.z(control_qubit)
 
    if n > 1:
        with qc.if_test(parity_target):
            qc.x(target_qubit)
 
    if add_barriers is True:
        qc.barrier()
    return qc

def build_circuits_dyn(n: int) -> List[QuantumCircuit]:
    """ """
    assert n >= 0, "Error: n needs to be a non-negative integer"
    circuits = []
 
    qr = QuantumRegister(
        n + 2, name="q"
    )  # Circuit with n qubits between control and target
    cr = ClassicalRegister(
        2, name="cr"
    )  # Classical register for measuring long range CNOT
 
    k = int(n / 2)  # Number of Bell States to be used
    c1 = ClassicalRegister(
        k, name="c1"
    )  # Classical register needed for post processing
    c2 = ClassicalRegister(
        n - k, name="c2"
    )  # Classical register needed for post processing
 
    for i in range(n):
        if n > 1:
            circuit = QuantumCircuit(qr, cr, c1, c2, name="CNOT")
        elif n == 1:
            circuit = QuantumCircuit(qr, cr, c2, name="CNOT")
        elif n == 0:
            circuit = QuantumCircuit(qr, cr, name="CNOT")
 
        circuit.h(0)
        circuit = CNOT_dyn(
            qc=circuit, control_qubit=0, target_qubit=i + 1, c1=c1, c2=c2
        )
 
        circuit.measure(0, cr[0])
        circuit.measure(i + 1, cr[1])
 
        circuits.append(circuit)
    return circuits

circuits_dyn = build_circuits_dyn(MAX_NUMBER_QUBITS)
circuits_dyn[-1].draw("mpl", fold=-1)

# Generate the ordered list of physical qubits in the chain
full_layout = set()
for edge in list(COUPLING_MAP_1D.get_edges()):
    full_layout.update(set(edge))
full_layout = list(full_layout)
 
isa_circuits_uni, isa_circuits_postproc, isa_circuits_dyn = [], [], []
for circuit_uni, circuit_dyn in zip(circuits_uni, circuits_dyn):
    # Generate the main Qiskit transpile passes.
    pm_uni = generate_preset_pass_manager(
        coupling_map=COUPLING_MAP_1D,
        initial_layout=full_layout[: circuit_uni.num_qubits],
        optimization_level=OPTIMIZATION_LEVEL,
        backend=backend,
    )
 
    pm_dyn = generate_preset_pass_manager(
        coupling_map=COUPLING_MAP_1D,
        initial_layout=full_layout[: circuit_dyn.num_qubits],
        optimization_level=OPTIMIZATION_LEVEL,
        backend=backend,
    )
 
    isa_circuits_uni.append(pm_uni.run(circuit_uni))
    isa_circuits_dyn.append(pm_dyn.run(circuit_dyn))

isa_circuits_dyn[4].draw("mpl", fold=-1, idle_wires=False)

sampler = Sampler(mode=backend)
 
sampler.options.experimental = {
    "execution_path": "gen3-experimental",
}
 
job_uni = sampler.run(isa_circuits_uni, shots=SHOTS)
job_dyn = sampler.run(isa_circuits_dyn, shots=SHOTS)

from qiskit.quantum_info import hellinger_fidelity
 
BELL_DIST = {"00": 0.5, "11": 0.5}
 
 
def get_counts_from_bitarray(instance):
    """
    Extract counts from result data
    """
    for field, value in instance.__dict__.items():
        if isinstance(value, BitArray):
            return value.get_counts()
    return None
 
 
def resample_single_dictionary(d):
    """Resample a single dictionary based on its weights."""
    keys = list(d.keys())
    weights = list(d.values())
    sum_weights = sum(weights)
 
    if sum_weights == 0:
        return d
 
    resampled_keys = random.choices(keys, weights=weights, k=sum_weights)
 
    # Count the occurrences of each key in the resampled keys
    resampled_counts = {}
    for key in resampled_keys:
        resampled_counts[key] = resampled_counts.get(key, 0) + 1
 
    return resampled_counts
 
 
def resample_dict_list(dict_list, n_samples):
    """Resample the entire list of dictionaries n_samples times."""
    resampled_lists = []
 
    for _ in range(n_samples):
        new_version = [resample_single_dictionary(d) for d in dict_list]
        resampled_lists.append(new_version)
 
    return resampled_lists

result_uni = job_uni.result()
 
uni_state_fid_mean, uni_state_fid_std = [], []
for i in range(len(result_uni)):
    counts = get_counts_from_bitarray(result_uni[i].data)
 
    state_fid_temp = []
    for _ in range(BOOTSTRAP_SAMPLES):
        resampled_counts = resample_single_dictionary(counts)
        fid = hellinger_fidelity(resampled_counts, BELL_DIST)
        state_fid_temp.append(fid)
 
    mean, std = (
        np.mean(np.array(state_fid_temp)),
        np.std(np.array(state_fid_temp)),
    )
    uni_state_fid_mean.append(mean)
    uni_state_fid_std.append(std)
 
 
print("State fidelities:")
print(["{0:0.3f}".format(i) for i in uni_state_fid_mean])
print("State fidelities std:")
print(["{0:0.3f}".format(i) for i in uni_state_fid_std])

result_dyn = job_dyn.result()
 
dyn_state_fid_mean, dyn_state_fid_std = [], []
for i in range(len(result_dyn)):
    counts = get_counts_from_bitarray(result_dyn[i].data)
 
    state_fid_temp = []
    for _ in range(BOOTSTRAP_SAMPLES):
        resampled_counts = resample_single_dictionary(counts)
        fid = hellinger_fidelity(resampled_counts, BELL_DIST)
        state_fid_temp.append(fid)
 
    mean, std = (
        np.mean(np.array(state_fid_temp)),
        np.std(np.array(state_fid_temp)),
    )
    dyn_state_fid_mean.append(mean)
    dyn_state_fid_std.append(std)
 
 
print("State fidelities:")
print(["{0:0.3f}".format(i) for i in dyn_state_fid_mean])
print("State fidelities std:")
print(["{0:0.3f}".format(i) for i in dyn_state_fid_std])

fig, ax = plt.subplots()
 
ax.errorbar(
    range(MIN_NUMBER_QUBITS, MAX_NUMBER_QUBITS),
    uni_state_fid_mean,
    uni_state_fid_std,
    fmt="o-.",
    color="c",
    label="Unitary",
)
 
ax.errorbar(
    range(MIN_NUMBER_QUBITS, MAX_NUMBER_QUBITS),
    dyn_state_fid_mean,
    dyn_state_fid_std,
    fmt="o-.",
    color="m",
    label="Dynamic",
)
ax.axhline(y=1 / 4, color="g", linestyle="--", label="Random gate")
legend = ax.legend(frameon=True)
for text in legend.get_texts():
    text.set_color("black")  # Set the legend text color to black
legend.get_frame().set_facecolor(
    "white"
)  # Set the legend background color to white
legend.get_frame().set_edgecolor(
    "black"
)  # Optional: set the legend border color to black
ax.set_xlabel("Number of qubits between control and target", color="black")
ax.set_ylabel("Bell state fidelity", color="black")
ax.grid(linestyle=":", linewidth=0.6, alpha=0.4, color="gray")
ax.set_ylim((0.2, 1))
ax.set_facecolor("white")  # Set the background color of the axes
fig.patch.set_facecolor("white")  # Set the background color of the figure
 
# Ensure the axis lines and ticks are visible
for spine in ax.spines.values():
    spine.set_visible(True)
    spine.set_color("black")  # Set the color of the axis lines to black
ax.tick_params(axis="x", colors="black")
ax.tick_params(axis="y", colors="black")
 
plt.show()