Real-Time Benchmarking for Qubit Selection

import pyscf
import pyscf.cc
import pyscf.mcscf
from pyscf.data.elements import chemcore
 
import rustworkx as rx
from rustworkx import NoEdgeBetweenNodes, PyGraph
from typing import List, Tuple
import numpy as np
import pandas as pd
import ffsim
from datetime import datetime, timedelta
import random
from collections import Counter
import matplotlib.pyplot as plt
import re
 
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.transpiler import CouplingMap
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passmanager import PassManager
from qiskit.transpiler.passes import RemoveIdentityEquivalent
from qiskit.providers import BackendV2
from qiskit.result import marginal_counts as mcts
 
from qiskit_experiments.framework import BatchExperiment, ParallelExperiment
from qiskit_experiments.library import StandardRB
from qiskit_experiments.library.randomized_benchmarking import LayerFidelity
 
import qiskit_device_benchmarking.utilities.graph_utils as gu
import qiskit_device_benchmarking.utilities.layer_fidelity_utils as lfu
 
from qiskit_ibm_runtime import QiskitRuntimeService, Session
from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime.models import BackendProperties
from qiskit_ibm_runtime.transpiler.passes.basis.fold_rzz_angle import (
    FoldRzzAngle,
)

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, use_fractional_gates=True
)
device = backend.name
# print(device)

# Create readout experiment on all qubits in parallel
num_qubits = backend.num_qubits
spam0 = QuantumCircuit(num_qubits, num_qubits)
spam0.measure_all()
spam1 = QuantumCircuit(num_qubits, num_qubits)
spam1.x(range(num_qubits))
spam1.measure_all()

# Create SQRB experiments in batches
G = backend.coupling_map.graph.to_undirected(multigraph=False)
sqrb_batches = gu.get_iso_qubit_list(G)
 
lengths = np.array([1, 50, 100, 500, 1000, 3000])  # Clifford lengths (x-axis)
num_samples = 6  # number of samples per clifford length (how many samplings are taken at each clifford length)
 
sqrb_exp_list = []
for batch in sqrb_batches:
    rb1q_exps = []
    for qubit in batch:
        rb1q_exp = StandardRB(
            physical_qubits=[int(qubit)],
            lengths=lengths,
            backend=backend,
            seed=42,
            num_samples=num_samples,
        )
        rb1q_exps.append(rb1q_exp)
    sqrb_exp_list.append(
        ParallelExperiment(rb1q_exps, backend=backend, flatten_results=True)
    )
sqrb_exp = BatchExperiment(
    sqrb_exp_list, backend=backend, flatten_results=True
)
sqrb_exp.set_experiment_options(separate_jobs=True)
samples_m = 3  # multiple of samples, determines number of circuits per job
sqrb_exp.experiment_options.max_circuits = samples_m * num_samples

# Create Layer Fidelity experiment
 
# Get two qubit gate
if "ecr" in backend.configuration().basis_gates:
    twoq_gate = "ecr"
elif "cz" in backend.configuration().basis_gates:
    twoq_gate = "cz"
else:
    twoq_gate = "cx"
 
# Get one qubit basis gates
oneq_gates = []
for i in backend.configuration().basis_gates:
    # put in a case to handle rx and rzz
    if i.casefold() == "rx" or i.casefold() == "rzz":
        continue
    if i.casefold() != twoq_gate.casefold():
        oneq_gates.append(i)
 
# First, get the grid chains, these are hard coded in the layer fidelity utilities module
grid_chains = lfu.get_grids(backend)
 
# There are two sets of chains that can be run in four disjoint experiments
coupling_map = backend.coupling_map
edges = list(backend.target[twoq_gate].keys())
layers = [[] for i in range(4)]
grid_chain_flt = [[], []]
for i in range(2):
    all_pairs = gu.path_to_edges(grid_chains[i], coupling_map)
    for j, pair_lst in enumerate(all_pairs):
        grid_chain_flt[i] += grid_chains[i][j]
        sub_pairs = [
            tuple(pair) if tuple(pair) in edges else tuple(pair)[::-1]
            for pair in pair_lst
        ]
        layers[2 * i] += sub_pairs[0::2]
        layers[2 * i + 1] += sub_pairs[1::2]
 
# LF for horizontally-trending chains
h_qubits = grid_chain_flt[0]
lf_h = LayerFidelity(
    physical_qubits=h_qubits,
    two_qubit_layers=[layers[0], layers[1]],
    lengths=[1, 10, 20, 30, 40, 60, 80, 100, 150, 200, 400],
    backend=backend,
    num_samples=12,
    seed=60,
    two_qubit_gate=twoq_gate,
    one_qubit_basis_gates=oneq_gates,
)
 
# LF for vertically-trending chains
v_qubits = grid_chain_flt[1]
lf_v = LayerFidelity(
    physical_qubits=v_qubits,
    two_qubit_layers=[layers[2], layers[3]],
    lengths=[1, 10, 20, 30, 40, 60, 80, 100, 150, 200, 400],
    backend=backend,
    num_samples=12,
    seed=60,
    two_qubit_gate=twoq_gate,
    one_qubit_basis_gates=oneq_gates,
)
 
# Set maximum number of circuits per job to avoid errors due to too large payload
lf_v.experiment_options.max_circuits = 144
lf_h.experiment_options.max_circuits = 144
 
batches_lf = [lf_v, lf_h]

# Let's retrieve the properties over a certain time range
days = 7
errors_list = []
for day_idx in range(0, days):
    calibrations_time = datetime.now() - timedelta(days=day_idx)
    targer_hist = backend.target_history(datetime=calibrations_time)
 
    t1_dict, t2_dict = {}, {}
    for qubit in range(targer_hist.num_qubits):
        t1_dict[qubit] = targer_hist.qubit_properties[qubit].t1
        t2_dict[qubit] = targer_hist.qubit_properties[qubit].t2
 
    errors_dict = {
        "1q": targer_hist["sx"],
        "2q": targer_hist[twoq_gate],
        "spam": targer_hist["measure"],
    }
 
    errors_list.append(errors_dict)

# Let's plot these properties over time
fig, axs = plt.subplots(3, 1, figsize=(10, 20), sharex=False)
 
# Plot readout values
for qubit in range(targer_hist.num_qubits):
    spams = []
    for errors_dict in errors_list:
        spam_dict = errors_dict["spam"]
        spams.append(spam_dict[tuple([qubit])].error)
 
    axs[0].plot(spams)
 
axs[0].set_title("Readout Errors")
axs[0].set_ylabel("Error Rate")
axs[0].set_xlabel("Days")
 
# Plot 1Q Gate Errors
for qubit in range(targer_hist.num_qubits):
    oneq_gates = []
    for errors_dict in errors_list:
        oneq_gate_dict = errors_dict["1q"]
        oneq_gates.append(oneq_gate_dict[tuple([qubit])].error)
 
    axs[1].plot(oneq_gates)
 
axs[1].set_title("1Q Gate Errors")
axs[1].set_ylabel("Error Rate")
axs[1].set_xlabel("Days")
axs[1].set_ylim([0, 0.01])
 
# Plot 2Q Gate Errors
for pair in edges:
    twoq_gates = []
    for errors_dict in errors_list:
        twoq_gate_dict = errors_dict["2q"]
        twoq_gates.append(twoq_gate_dict[pair].error)
 
    axs[2].plot(twoq_gates)
 
axs[2].set_title("2Q Gate Errors")
axs[2].set_ylabel("Error Rate")
axs[2].set_xlabel("Days")
axs[2].set_ylim([0, 0.15])
 
plt.subplots_adjust(hspace=0.5)
plt.show()

# Specify molecule properties
open_shell = False
spin_sq = 0
 
# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="6-31g",
    symmetry="Dooh",
)
 
# Get the active space
ncore = chemcore(mol)
naorb = mol.nao_nr() - ncore
active_space = [p for p in range(ncore, mol.nao_nr())]
 
# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2
 
# Prepare UCJ operator
n_reps = 1
alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]
alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]
nelec = (num_elec_a, num_elec_b)
 
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),
)
orbital_rotation = ucj_op.orbital_rotations[0]
diag_coulomb_mat_aa, diag_coulomb_mat_ab = ucj_op.diag_coulomb_mats[0]
 
# Create an empty quantum circuit
qubits = QuantumRegister(2 * num_orbitals, name="q")
circuit = QuantumCircuit(qubits)
 
# Prepare Hartree-Fock state as the reference state and append it to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)
 
# Apply the UCJ operator to the reference state
circuit.append(
    ffsim.qiskit.OrbitalRotationJW(ucj_op.norb, orbital_rotation.T.conj()),
    qubits,
)
circuit.append(
    ffsim.qiskit.DiagCoulombEvolutionJW(
        ucj_op.norb,
        (diag_coulomb_mat_aa, diag_coulomb_mat_ab, diag_coulomb_mat_aa),
        -1.0,
    ),
    qubits,
)
 
# Add a barrier for transpilation
circuit.barrier()
 
# Apply the inverse UCJ operator
circuit.append(
    ffsim.qiskit.DiagCoulombEvolutionJW(
        ucj_op.norb,
        (diag_coulomb_mat_aa, diag_coulomb_mat_ab, diag_coulomb_mat_aa),
        1.0,
    ),
    qubits,
)
circuit.append(
    ffsim.qiskit.OrbitalRotationJW(ucj_op.norb, orbital_rotation), qubits
)
 
# Measure qubits
circuit.measure_all()

# Let's draw the circuit
circuit.draw(fold=-1)

def create_linear_chains(num_orbitals: int) -> PyGraph:
    """
    In zig-zag layout, there are two linear chains (with connecting qubits between
    the chains). This function creates those two linear chains: a rustworkx PyGraph
    with two disconnected linear chains. Each chain contains `num_orbitals` number
    of nodes, i.e., in the final graph there are `2 * num_orbitals` number of nodes.
 
    Args:
        num_orbitals (int): Number orbitals or nodes in each linear chain. They are
            also known as alpha-alpha interaction qubits.
 
    Returns:
        A rustworkx.PyGraph with two disconnected linear chains each with `num_orbitals`
            number of nodes.
    """
    G = rx.PyGraph()
 
    for n in range(num_orbitals):
        G.add_node(n)
 
    for n in range(num_orbitals - 1):
        G.add_edge(n, n + 1, None)
 
    for n in range(num_orbitals, 2 * num_orbitals):
        G.add_node(n)
 
    for n in range(num_orbitals, 2 * num_orbitals - 1):
        G.add_edge(n, n + 1, None)
 
    return G
 
 
def create_lucj_zigzag_layout(
    num_orbitals: int, backend_coupling_graph: PyGraph
) -> tuple[PyGraph, int]:
    """
    This function creates the complete zigzag graph that 'can be mapped' to a IBM QPU with
    heavy-hex connectivity (the zigzag must be an isomorphic sub-graph to the QPU/backend
    coupling graph for it to be mapped).
    The zigzag pattern includes both linear chains (alpha-alpha interactions) and connecting
    qubits between the linear chains (alpha-beta interactions).
 
    Args:
        num_orbitals (int): Number of orbitals, i.e., number of nodes in each alpha-alpha linear chain.
        backend_coupling_graph (PyGraph): Coupling map of the backend
 
    Returns:
        G_new (PyGraph): The graph with IBM backend compliant zigzag pattern.
        num_alpha_beta_qubits (int): Number of connecting qubits between the linear chains
            in the zigzag pattern.
    """
    isomorphic = False
    G = create_linear_chains(num_orbitals=num_orbitals)
 
    num_iters = num_orbitals
    while not isomorphic:
        G_new = G.copy()
        num_alpha_beta_qubits = 0
        for n in range(num_iters):
            if n % 4 == 0:
                new_node = 2 * num_orbitals + num_alpha_beta_qubits
                G_new.add_node(new_node)
                G_new.add_edge(n, new_node, None)
                G_new.add_edge(new_node, n + num_orbitals, None)
                num_alpha_beta_qubits = num_alpha_beta_qubits + 1
        isomorphic = rx.is_subgraph_isomorphic(backend_coupling_graph, G_new)
        num_iters -= 1
 
    return G_new, num_alpha_beta_qubits
 
 
def _make_backend_cmap_pygraph(backend: BackendV2) -> PyGraph:
    graph = backend.coupling_map.graph
    if not graph.is_symmetric():
        graph.make_symmetric()
    backend_coupling_graph = graph.to_undirected()
 
    edge_list = backend_coupling_graph.edge_list()
    removed_edge = []
    for edge in edge_list:
        if set(edge) in removed_edge:
            continue
        try:
            backend_coupling_graph.remove_edge(edge[0], edge[1])
            removed_edge.append(set(edge))
        except NoEdgeBetweenNodes:
            pass
 
    return backend_coupling_graph
 
 
def get_zigzag_layouts(
    num_orbitals: int,
    backend: BackendV2,
) -> List[List[int]]:
    """
    Function that generates a list of zigzag pattern layouts from the
    coupling map and physical qubits of a given backend.
 
    Args:
        num_orbitals (int): Number of orbitals.
        backend (BackendV2): A backend.
 
    Returns:
        A list of device compliant layouts (list[list[int]]) with zigzag pattern
    """
    backend_coupling_graph = _make_backend_cmap_pygraph(backend)
 
    G, num_aplha_beta_qubits = create_lucj_zigzag_layout(
        num_orbitals=num_orbitals,
        backend_coupling_graph=backend_coupling_graph,
    )
 
    # Instead of generating many isomorphic layouts, we can generate only one
    # and return that as later we will be using `VF2PostLayout` to find and score
    # many isomorphic layouts.
    isomorphic_mappings = rx.vf2_mapping(
        backend_coupling_graph, G, subgraph=True
    )
    isomorphic_mappings = list(isomorphic_mappings)
 
    edges = list(G.edge_list())
 
    layouts = []
    for mapping in isomorphic_mappings:
        initial_layout = [None] * (2 * num_orbitals + num_aplha_beta_qubits)
        for key, value in mapping.items():
            initial_layout[value] = key
        layouts.append(initial_layout[:-num_aplha_beta_qubits])
    return layouts
 
 
def get_edges(path: List[int], backend) -> List[Tuple[int]]:
    """
    Return a list of edges from a path of qubits.
 
    Args:
    - path: List of nodes (qubits)
    - backend: a backend
    Returns:
    - List of edges
    """
 
    coupling_map = CouplingMap(backend.configuration().coupling_map)
    G = coupling_map.graph
    edges = []
    prev_node = None
    # get all prossible edges from the nodes provided
    for node in path:
        if prev_node is not None:
            # check the edge is in G and in the coupling map
            if (
                G.has_edge(prev_node, node)
                and (prev_node, node) in coupling_map
            ):
                edges.append((prev_node, node))
            elif (
                G.has_edge(node, prev_node)
                and (node, prev_node) in coupling_map
            ):
                edges.append((node, prev_node))
        prev_node = node
    return edges

# find all possible zig-zag layouts
all_layouts = get_zigzag_layouts(num_orbitals, backend)
print("Number of possible zig zag layouts:", len(all_layouts))
 
# randomly select a subset of N layouts
N = 30 if len(all_layouts) >= 30 else len(all_layouts)
subsets = random.sample(all_layouts, N)
print("Number of randomly selected layouts:", len(subsets))

# # Save randomly selected layouts locally (optional)
# today = date.today()
# print('Date:', today)
 
# subsets_dict = {}
# for i, v in enumerate(subsets):
#     subsets_dict[i] = v
 
# with open(f'{today}_{device}_{num_orbitals}_subsets.yaml', 'w') as file:
#     yaml.dump(subsets_dict, file)

transpiled_circuits = []
 
for i, layout in enumerate(subsets):
    # Transpile circuit on different layouts
    pass_manager = generate_preset_pass_manager(
        optimization_level=3,
        backend=backend,
        initial_layout=layout,
        seed_transpiler=42,
    )
    # Set the pre-initialization stage of the pass manager with passes suggested by ffsim
    pass_manager.pre_init = ffsim.qiskit.PRE_INIT
    pass_manager.post_init = PassManager([RemoveIdentityEquivalent()])
    pass_manager.post_optimization = PassManager(
        [FoldRzzAngle(), RemoveIdentityEquivalent(target=backend.target)]
    )
    transpiled = pass_manager.run(circuit)
    transpiled_circuits.append(transpiled)
    print("Counts for layout", i + 1, transpiled.count_ops())

# Helper functions to update QPU properties
 
# Make a readout error map
def get_readout_errors(num_qubits, readout_result):
    ro_error = {}
    cts_spam0 = readout_result[0].data.meas.get_counts()
    cts_spam1 = readout_result[1].data.meas.get_counts()
    spam_shots = readout_result[0].data.meas.num_shots
    for q in range(num_qubits):
        try:
            ro_error[q] = (
                1
                - (
                    (mcts(cts_spam0, [q])["0"] + mcts(cts_spam1, [q])["1"])
                    / 2
                )
                / spam_shots
            )
        except KeyError:
            ro_error[q] = 1 - ((mcts(cts_spam0, [q])["0"]) / 2) / spam_shots
    return ro_error
 
 
# Make a 1Q error map
def get_error_dict(job_results):
    values_map = {}
    for index, row in job_results.iterrows():
        if len(row.components) > 1:
            q1 = row.components[0].index
            q2 = row.components[1].index
            q = (q1, q2)
        else:
            q = row.components[0].index
        value = row.value.nominal_value
        values_map[q] = value
    return values_map
 
 
# Make a 2Q error map
def get_twoq_errors(backend, lf_result_h, lf_result_v):
    # Get two qubit gate
    if "ecr" in backend.configuration().basis_gates:
        twoq_gate = "ecr"
    elif "cz" in backend.configuration().basis_gates:
        twoq_gate = "cz"
    else:
        twoq_gate = "cx"
    lf_err_dict = lfu.make_error_dict(backend, twoq_gate)
    updated_err_dicts = []
    for i, lf_data in enumerate([lf_result_h, lf_result_v]):
        for j in range(2):
            updated_err_dicts.append(
                lfu.df_to_error_dict(lf_data, layers[2 * i + j])
            )
    lf_err_dict = lfu.update_error_dict(lf_err_dict, updated_err_dicts)
    return lf_err_dict
 
 
# Update readout errors
def update_readout_props(props, error_map):
    qubits = error_map.keys()
    for i in range(len(qubits)):
        for param in props["qubits"][i]:
            if param["name"] == "readout_error":
                param["value"] = error_map[i]
                break
    return props
 
 
# Update 1Q errors
def update_1q_props(props, error_map):
    for ix, data in enumerate(props["gates"]):
        # Determine if referenced gate is 2Q
        if len(data["qubits"]) > 1:
            # Get the gate type and hardware component acted on by the referenced gate
            x = re.search("^([a-zA-Z]+)([0-9]+_[0-9]+)", data["name"])
            gate_type, component = (
                x[1],
                (int(x[2].split("_")[0]), int(x[2].split("_")[1])),
            )
        # Determine that referenced gate is 1Q
        else:
            # Get the gate type and hardware component acted on by the referenced gate
            x = re.search("^([a-zA-Z]+)([0-9]+)", data["name"])
            gate_type, component = x[1], int(x[2])
 
        # Scan over gate parameters
        if gate_type == "sx" or gate_type == "x":
            for iy, parameter in enumerate(data["parameters"]):
                # Check if a parameter stores a gate error rate
                if parameter["name"] == "gate_error":
                    # Check if that parameter is in the characterized error rates
                    if (
                        f"EPG_{gate_type}" == "EPG_sx"
                        or f"EPG_{gate_type}" == "EPG_x"
                    ):
                        gate_error = error_map[component]
                        props["gates"][ix]["parameters"][iy]["value"] = (
                            gate_error
                        )
    return props
 
 
# Update layer fidelity 2Q errors
def update_lf_props(props, error_map):
    for ix, gate in enumerate(props["gates"]):
        if gate["gate"] in ["cz", "cx", "ecr"]:
            qubit_pair = gate["qubits"]
            pair_key = f"{qubit_pair[0]}_{qubit_pair[1]}"
            pair_key_reversed = f"{qubit_pair[1]}_{qubit_pair[0]}"
            if pair_key in error_map or pair_key_reversed in error_map:
                try:
                    gate_error = error_map[pair_key]
                except:
                    gate_error = error_map[pair_key_reversed]
                for iy, param in enumerate(gate["parameters"]):
                    if param["name"] == "gate_error":
                        props["gates"][ix]["parameters"][iy]["value"] = (
                            gate_error
                        )
                        break
    return props

# Run circuits
with Session(backend=backend) as session:
    # Run characterization experiments
    sampler = Sampler()
 
    # QPU usage ~8s
    sampler.options.default_shots = 10_000
    job_readout = sampler.run([spam0, spam1])
    print("Readout job id:", job_readout.job_id())
 
    # QPU usage 16s
    sampler.options.default_shots = 250
    job_sqrb = sqrb_exp.run(sampler=sampler)
    print("Single qubit RB job ids:", job_sqrb.job_ids)
 
    # QPU usage: ~15s each with the use of fractional gates, ~60s total
    sampler.options.default_shots = 250
    job_lf_vert = lf_v.run(sampler=sampler)
    job_lf_ho = lf_h.run(sampler=sampler)
    print("Vertical Layer Fidelity job ids:", job_lf_vert.job_ids)
    print("Horizontal Layer Fidelity job ids:", job_lf_ho.job_ids)
 
    # Retrieve characterization results
    readout_result = job_readout.result()
    sqrb_result_x = job_sqrb.analysis_results("EPG_x", dataframe=True)
    sqrb_result_sx = job_sqrb.analysis_results("EPG_sx", dataframe=True)
    lf_result_v = job_lf_vert.analysis_results(
        "ProcessFidelity", dataframe=True
    )
    lf_result_h = job_lf_ho.analysis_results(
        "ProcessFidelity", dataframe=True
    )
 
    # Create error maps from all characterization results
    ro_error = get_readout_errors(backend.num_qubits, readout_result)
    sqrb_err_x_dict = get_error_dict(sqrb_result_x)
    sqrb_err_sx_dict = get_error_dict(sqrb_result_sx)
    lf_err_dict = get_twoq_errors(backend, lf_result_h, lf_result_v)
 
    # Update QPU properties
    backend_updated = service.backend(device, use_fractional_gates=True)
    props_dict = backend_updated.properties().to_dict()
    props_dict = update_readout_props(props_dict, ro_error)
    props_dict = update_1q_props(props_dict, sqrb_err_x_dict)
    props_dict = update_1q_props(props_dict, sqrb_err_sx_dict)
    props_dict = update_lf_props(props_dict, lf_err_dict)
    props = BackendProperties.from_dict(props_dict)
    backend_updated._properties = props
 
    # # Transpile circuit using updated QPU properties to find the best qubit layout
    # pass_manager = generate_preset_pass_manager(optimization_level=3, backend=backend_updated)
    # pass_manager.pre_init = ffsim.qiskit.PRE_INIT
    # pass_manager.post_init = PassManager([RemoveIdentityEquivalent()])
    # pass_manager.post_optimization = PassManager(
    #     [
    #         FoldRzzAngle(), RemoveIdentityEquivalent(target=backend.target)
    #     ]
    # )
    # transpiled_circuits = pass_manager.run(circuit) # for some circuit of interest
 
    # Run your circuit (in this case the modified LUCJ circuit) and define sampler options to turn DD on/off
    sampler_options = {
        "dynamical_decoupling": {
            "enable": True,
            "sequence_type": "XpXm",
        }
    }
    sampler = Sampler(options=sampler_options)
    sampler.options.default_shots = 200_000
    job_sampler = sampler.run(transpiled_circuits)
    print("Job id for list of circuits:", job_sampler.job_id())

# readout error dictionary
ro_error

print("SPAM error for qubit 0:", backend.properties().readout_error(0))
print(
    "Updated SPAM error for qubit 0:",
    backend_updated.properties().readout_error(0),
)

# # Optionally, we can save these backend properties for future analysis
 
# props = backend.properties().to_dict()
# props_lf = backend_updated.properties().to_dict()
 
# today = date.today()
 
# try:
#     with open(f'{today}_{device}_properties.yaml', 'w') as file:
#         yaml.dump(props, file, sort_keys=False)
# except Exception as e:
#     print(f"An error occurred: {e}")
 
 
# try:
#     with open(f'{today}_{device}_properties_realtime_lf.yaml', 'w') as file:
#         yaml.dump(props_lf, file, sort_keys=False)
# except Exception as e:
#     print(f"An error occurred: {e}")

results = job_sampler.result()

def get_correct_bitstrings(counts, hf_bitstring):
    """This function returns the number of bitstrings with the correct
    Hartree-Fock bitstring.
 
    Inputs:
        - counts (Dict): dictionary of bitstrings with frequency counts as values
        - hf_bitstring (str): string representing a Hartree-Fock state
    Output:
        - correct_Bs (int): number of bitstrings with the correct
                            Hartree-Fock bitstring
    """
    correct_bs = 0
    for bitstring, freq in counts.items():
        if bitstring == hf_bitstring:
            correct_bs += freq
    return correct_bs

# Hartree-Fock state bitstring
hf_bitstring = (
    "0" * (num_orbitals - num_elec_b)
    + "1" * num_elec_b
    + "0" * (num_orbitals - num_elec_a)
    + "1" * num_elec_a
)
 
print("Hartree-Fock bitstring", hf_bitstring)

# Get ratios of Hartreefock states
ratios = []
for i, layout in enumerate(subsets):
    result = results[i]
    counts = result.data.meas.get_counts()
    correct_bs = get_correct_bitstrings(counts, hf_bitstring)
    shots = result.data.meas.num_shots
    ratio = correct_bs / shots
    ratios.append(ratio)
    print(f"{i+1} ratio of correct Hartree-Fock bitstrings:", ratio)

def qiskit_cost(backend, circuit):
    # Get the counts of 1q and 2q for all qubits and edges in the circuit
    gate_counts = {}
    two_q_list = []
    one_q_list = []
    for item in circuit._data:
        if (
            item.operation.num_qubits == 2
            and item.operation.name != "barrier"
        ):
            q0 = circuit.find_bit(item.qubits[0]).index
            q1 = circuit.find_bit(item.qubits[1]).index
            if item.operation.name in ["cz", "rzz", "ecr" "cx"]:
                two_q_list.append((q0, q1))
        elif item.operation.name in [
            "rx",
            "x",
            "sx",
            "measure",
            "reset",
            "id",
        ]:
            q0 = circuit.find_bit(item.qubits[0]).index
            one_q_list.append(q0)
    gate_counts = {
        "twoq": (Counter(two_q_list)),
        "oneq": (Counter(one_q_list)),
    }
 
    # Get the average 1q error for all qubits in the circuit
    props = backend.properties()
    oneq_errors = {}
    oneq_gates = list(
        set(backend.configuration().basis_gates)
        - set(["cz", "rzz", "ecr", "cx"])
    )
    for q in one_q_list:
        ave_err = 0
        for gate in oneq_gates:
            ave_err += props.gate_error(gate, q)
        ave_err += props.readout_error(q)
        ave_err = ave_err / (
            len(oneq_errors) + 1
        )  # take average err of all 1q operations
        oneq_errors[q] = ave_err
 
    # Get the average 2q error for all edges in the circuit
    twoq_errors = {}
    twoq_gates = list(
        set(backend.configuration().basis_gates) - set(oneq_gates)
    )
    for pair in two_q_list:
        q0 = pair[0]
        q1 = pair[1]
        ave_err = 0
        for gate in twoq_gates:
            try:
                ave_err += props.gate_error(gate, [q0, q1])
            except:
                ave_err += props.gate_error(gate, [q1, q0])
        ave_err = ave_err / len(twoq_gates)
        twoq_errors[pair] = ave_err
 
    # Compute the product of fidelities from 1q and 2q errors
    fidelity = 1
    for k in gate_counts.keys():
        for indx, num_gates in gate_counts[k].items():
            if k == "twoq":
                err = twoq_errors[indx]
            elif k == "oneq":
                err = oneq_errors[indx]
            fidelity *= (1 - err) ** num_gates
    return fidelity

# Find the estimated cost of each layout
scores = []
for i, layout in enumerate(subsets):
    reported = qiskit_cost(backend, circuit=transpiled_circuits[i])
    lf = qiskit_cost(backend_updated, circuit=transpiled_circuits[i])
    scores.append(
        {
            "layout": i,
            "cost_reported": reported,
            "cost_lf": lf,
            "ratio": ratios[i],
        }
    )
scores = pd.DataFrame(scores)
 
scores = scores.sort_values(by="ratio", ascending=False).reset_index(
    drop=True
)

# Function to normalize an array using Min-Max Scaling
def normalize(x):
    x = np.asarray(x)
    return (x - np.min(x)) / (np.max(x) - np.min(x))
 
 
# Get costs data
x = np.asarray(scores.index.to_list()) + 1
r = scores.ratio
repo = scores.cost_reported
lf = scores.cost_lf
 
# Normaliza data
repo = normalize(repo)
lf = normalize(lf)
 
# Create the figure and primary y-axis
fig, ax1 = plt.subplots(figsize=(12, 5))
 
# Create the secondary y-axis
ax2 = ax1.twinx()
 
# Plot the data
ax1.plot(x, repo, marker="o", linestyle="-", label="reported")
ax1.plot(x, lf, marker="o", linestyle="-", label="real-time", color="crimson")
ax2.plot(
    x,
    r,
    color="black",
    marker="x",
    linestyle="--",
    label="percetange of correct bitstrings",
)
 
# Set labels and titles
ax1.set_xlabel("Layout Index")
ax1.set_xticks(x)
ax1.set_ylabel("Layout Score")
ax2.set_ylabel("Ratio of valid HF states")
ax1.legend(loc="upper left", bbox_to_anchor=(1.15, 0.5))
ax2.legend(loc="lower left", bbox_to_anchor=(1.15, 0.5))
ax1.grid()
plt.xlabel("Layout number")
fig.tight_layout()
plt.title(
    "Modified LUCJ circuit: Correlation between layout score and hardware performance for different QPU properties"
)

valid_layouts = []
invalid_layouts = []
invalid_edges = []
invalid_qubits = []
for i, layout in enumerate(all_layouts):
    errors = []
    edges = get_edges(layout, backend)
    for edge in edges:
        error = backend.properties().gate_error(twoq_gate, edge)
        errors.append(error)
        if error == 1:
            invalid_edges.append(edge) if edge not in invalid_edges else None
    for qubit in layout:
        error = backend.properties().gate_error("sx", qubit)
        errors.append(error)
        if error == 1:
            invalid_qubits.append(
                qubit
            ) if qubit not in invalid_qubits else None
    if all(err < 1 for err in errors):
        valid_layouts.append(layout)
    else:
        invalid_layouts.append(layout)
 
print("Number of valid layouts (all errors < 1):", len(valid_layouts))
print(
    "Number of invalid layouts (one or more errors == 1):",
    len(invalid_layouts),
)
print("List of invalid 1Q gates:", invalid_qubits)
print("List of invalid 2Q gates:", invalid_edges)
if invalid_qubits:
    print("Quick example:")
    print(
        f"1Q error for {invalid_qubits[0]}:",
        backend.properties().gate_error("sx", invalid_qubits[0]),
    )
if invalid_edges:
    print(
        f"2Q error for {invalid_edges[0]}:",
        backend.properties().gate_error(twoq_gate, invalid_edges[0]),
    )

if invalid_layouts:
    pass_manager_invalid = generate_preset_pass_manager(
        optimization_level=3,
        backend=backend,
        initial_layout=invalid_layouts[8],
        seed_transpiler=42,
    )
 
    # Set the pre-initialization stage of the pass manager with passes suggested by ffsim
    pass_manager_invalid.pre_init = ffsim.qiskit.PRE_INIT
    pass_manager_invalid.post_init = PassManager([RemoveIdentityEquivalent()])
    pass_manager_invalid.post_optimization = PassManager(
        [FoldRzzAngle(), RemoveIdentityEquivalent(target=backend.target)]
    )
    transpiled_invalid = pass_manager_invalid.run(circuit)
 
    print(
        "Score from reported QPU properties:",
        qiskit_cost(backend, transpiled_invalid),
    )
    print(
        "Score from real-time QPU properties:",
        qiskit_cost(backend_updated, transpiled_invalid),
    )