Debug Qiskit Runtime jobs

import numpy as np
import random
 
from qiskit.circuit import QuantumCircuit
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.quantum_info import SparsePauliOp
 
from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2 as Estimator
from qiskit_ibm_runtime.debug_tools import Neat
 
from qiskit_aer.noise import NoiseModel, depolarizing_error

# Choose the least busy backend
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
 
# Generate a preset pass manager
# This will be used to convert the abstract circuit to an equivalent Instruction Set Architecture (ISA) circuit.
pm = generate_preset_pass_manager(backend=backend, optimization_level=0)
 
# Set the random seed
random.seed(10)

def generate_circuit(n_qubits, n_layers):
    r"""
    A function to generate a pseudo-random a circuit with ``n_qubits`` qubits and
    ``2*n_layers`` entangling layers of the type used in this notebook.
    """
    # An array of random angles
    angles = [
        [random.random() for q in range(n_qubits)] for s in range(n_layers)
    ]
 
    qc = QuantumCircuit(n_qubits)
    qubits = list(range(n_qubits))
 
    # do random circuit
    for layer in range(n_layers):
        # rotations
        for q_idx, qubit in enumerate(qubits):
            qc.rz(angles[layer][q_idx], qubit)
 
        # cx gates
        control_qubits = (
            qubits[::2] if layer % 2 == 0 else qubits[1 : n_qubits - 1 : 2]
        )
        for qubit in control_qubits:
            qc.cx(qubit, qubit + 1)
 
    # undo random circuit
    for layer in range(n_layers)[::-1]:
        # cx gates
        control_qubits = (
            qubits[::2] if layer % 2 == 0 else qubits[1 : n_qubits - 1 : 2]
        )
        for qubit in control_qubits:
            qc.cx(qubit, qubit + 1)
 
        # rotations
        for q_idx, qubit in enumerate(qubits):
            qc.rz(-angles[layer][q_idx], qubit)
 
    return qc
 
 
# Generate a random circuit
qc = generate_circuit(6, 3)
# Convert the abstract circuit to an equivalent ISA circuit.
isa_qc = pm.run(qc)
 
qc.draw("mpl", idle_wires=0)

# Initialize the observables
obs = ["ZIIIII", "IZIIII", "IIZIII", "IIIZII", "IIIIZI", "IIIIIZ"]
print(f"Observables: {obs}")
 
# Map the observables to the backend's layout
isa_obs = [SparsePauliOp(o).apply_layout(isa_qc.layout) for o in obs]
 
# Initialize the PUBs, which consist of six-qubit circuits with `n_layers` 1, ..., 6
all_n_layers = [1, 2, 3, 4, 5, 6]
 
pubs = [(pm.run(generate_circuit(6, n)), isa_obs) for n in all_n_layers]

# You could specify a custom `NoiseModel` here. If `None`, `Neat`
# pulls the noise model from the given backend
noise_model = None
 
# Initialize `Neat`
analyzer = Neat(backend, noise_model)

clifford_pubs = analyzer.to_clifford(pubs)
 
clifford_pubs[0].circuit.draw("mpl", idle_wires=0)

# Perform a noiseless simulation
ideal_results = analyzer.ideal_sim(clifford_pubs)
print(f"Ideal results:\n {ideal_results}\n")
 
# Perform a noisy simulation with the backend's noise model
noisy_results = analyzer.noisy_sim(clifford_pubs)
print(f"Noisy results:\n {noisy_results}\n")

# Figure of merit: Absolute difference
def rdiff(res1, re2):
    r"""The absolute difference between `res1` and re2`.
 
    --> The closer to `0`, the better.
    """
    d = abs(res1 - re2)
    return np.round(d.vals * 100, 2)
 
 
for idx, (ideal_res, noisy_res) in enumerate(
    zip(ideal_results, noisy_results)
):
    vals = rdiff(ideal_res, noisy_res)
 
    # Print the mean absolute difference for the observables
    mean_vals = np.round(np.mean(vals), 2)
    print(
        f"Mean absolute difference between ideal and noisy results for circuits with {all_n_layers[idx]} layers:\n  {mean_vals}%\n"
    )

# Set up a noise model with strength 0.02 on every two-qubit gate
noise_model = NoiseModel()
for qubits in backend.coupling_map:
    noise_model.add_quantum_error(
        depolarizing_error(0.02, 2), ["ecr", "cx"], qubits
    )
 
# Update the analyzer's noise model
analyzer.noise_model = noise_model
 
# Perform a noiseless simulation
ideal_results = analyzer.ideal_sim(clifford_pubs)
 
# Perform a noisy simulation with the backend's noise model
noisy_results = analyzer.noisy_sim(clifford_pubs)
 
# Compare the results
for idx, (ideal_res, noisy_res) in enumerate(
    zip(ideal_results, noisy_results)
):
    values = rdiff(ideal_res, noisy_res)
 
    # Print the mean absolute difference for the observables
    mean_values = np.round(np.mean(values), 2)
    print(
        f"Mean absolute difference between ideal and noisy results for circuits with {all_n_layers[idx]} layers:\n  {mean_values}%\n"
    )

# Generate a circuit with six qubits and six layers
isa_qc = pm.run(generate_circuit(6, 3))
 
# Use the same observables as previously
pubs = [(isa_qc, isa_obs)]
clifford_pubs = analyzer.to_clifford(pubs)

noise_factors = [
    [1, 1.1],
    [1, 1.1, 1.2],
    [1, 1.5, 2],
    [1, 1.5, 2, 2.5, 3],
    [1, 4],
]

# Run the PUBs on a QPU
estimator = Estimator(backend)
estimator.options.default_shots = 100000
estimator.options.twirling.enable_gates = True
estimator.options.twirling.enable_measure = True
estimator.options.twirling.shots_per_randomization = 100
estimator.options.resilience.measure_mitigation = True
estimator.options.resilience.zne_mitigation = True
estimator.options.resilience.zne.amplifier = "pea"
 
jobs = []
for factors in noise_factors:
    estimator.options.resilience.zne.noise_factors = factors
    jobs.append(estimator.run(clifford_pubs))
 
results = [job.result() for job in jobs]

# Perform a noiseless simulation
ideal_results = analyzer.ideal_sim(clifford_pubs)

# Look at the mean absolute difference to quickly tell the best choice for your options
for factors, res in zip(noise_factors, results):
    d = rdiff(ideal_results[0], res[0])
    print(
        f"Mean absolute difference for factors {factors}:\n  {np.round(np.mean(d), 2)}%\n"
    )