Compilation methods for Hamiltonian simulation circuits

from qiskit.circuit import QuantumCircuit
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2
from qiskit.circuit.library import PauliEvolutionGate
from qiskit_ibm_transpiler import generate_ai_pass_manager
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passes.synthesis.high_level_synthesis import HLSConfig
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error
from collections import Counter
from statistics import mean, stdev
from scipy.sparse import SparseEfficiencyWarning
import time
import warnings
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import numpy as np
import json
import requests
import logging

# Suppress noisy loggers and warnings
logging.getLogger(
    "qiskit_ibm_transpiler.wrappers.ai_local_synthesis"
).setLevel(logging.ERROR)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=SparseEfficiencyWarning)

seed = 42  # Seed for reproducibility

# QiskitRuntimeService.save_account(channel="ibm_quantum_platform",
# token="<YOUR-API-KEY>", overwrite=True, set_as_default=True)
service = QiskitRuntimeService(channel="ibm_quantum_platform")
backend = service.least_busy(operational=True, simulator=False)
print(f"Using backend: {backend.name}")

# SABRE pass manager (Qiskit default at optimization level 3)
pm_sabre = generate_preset_pass_manager(
    optimization_level=3, backend=backend, seed_transpiler=seed
)

# AI transpiler pass manager (local mode)
pm_ai = generate_ai_pass_manager(
    backend=backend, optimization_level=3, ai_optimization_level=3
)

# Rustiq pass manager for PauliEvolutionGate synthesis
hls_config = HLSConfig(
    PauliEvolution=[
        (
            "rustiq",
            {
                "nshuffles": 400,
                "upto_phase": True,
                "fix_clifford": True,
                "preserve_order": False,
                "metric": "depth",
            },
        )
    ]
)
pm_rustiq = generate_preset_pass_manager(
    optimization_level=3,
    backend=backend,
    hls_config=hls_config,
    seed_transpiler=seed,
)

def capture_transpilation_metrics(
    results, pass_manager, circuits, method_name
):
    """
    Transpile circuits and append one metrics record per circuit to
    ``results``.

    Args:
        results (list): List of dicts to append the metrics records to.
        pass_manager: Pass manager used for transpilation.
        circuits (list): List of quantum circuits to transpile.
        method_name (str): Name of the transpilation method.

    Returns:
        list: List of transpiled circuits.
    """
    transpiled_circuits = []

    for i, qc in enumerate(circuits):
        start_time = time.time()
        transpiled_qc = pass_manager.run(qc)
        end_time = time.time()

        # Decompose swaps for consistency across methods
        transpiled_qc = transpiled_qc.decompose(gates_to_decompose=["swap"])

        transpilation_time = end_time - start_time
        two_qubit_depth = transpiled_qc.depth(
            lambda x: x.operation.num_qubits == 2
        )
        circuit_size = transpiled_qc.size()

        results.append(
            {
                "method": method_name,
                "qc_name": qc.name,
                "qc_index": i,
                "num_qubits": qc.num_qubits,
                "two_qubit_depth": two_qubit_depth,
                "size": circuit_size,
                "runtime": transpilation_time,
            }
        )
        transpiled_circuits.append(transpiled_qc)
        print(
            f"[{method_name}] Circuit {i} ({qc.name}): "
            f"2Q depth={two_qubit_depth}, size={circuit_size}, "
            f"time={transpilation_time:.2f}s"
        )

    return transpiled_circuits

def _method_order(results):
    """Return the distinct method names in their first-seen order."""
    order = []
    for r in results:
        if r["method"] not in order:
            order.append(r["method"])
    return order


def print_summary_table(results):
    """
    Print the mean and standard deviation of each metric per compilation
    method, followed by the mean percent improvement relative to SABRE.
    """
    metrics = [
        ("two_qubit_depth", "2Q Depth"),
        ("size", "Gate Count"),
        ("runtime", "Runtime (s)"),
    ]
    methods = _method_order(results)
    by_method = {m: [r for r in results if r["method"] == m] for m in methods}
    sabre_by_index = {r["qc_index"]: r for r in by_method.get("SABRE", [])}

    col_w = 22
    name_w = max(len(m) for m in methods)
    header = f"{'Method':<{name_w}}" + "".join(
        f"  {label:>{col_w}}" for _, label in metrics
    )

    print("Mean +/- std per compilation method")
    print(header)
    print("-" * len(header))
    for method in methods:
        cells = []
        for key, _ in metrics:
            values = [r[key] for r in by_method[method]]
            std = stdev(values) if len(values) > 1 else 0.0
            cells.append(f"{mean(values):,.1f} +/- {std:,.1f}")
        print(
            f"{method:<{name_w}}" + "".join(f"  {c:>{col_w}}" for c in cells)
        )

    others = [m for m in methods if m != "SABRE"]
    if others and sabre_by_index:
        print()
        print("Mean % improvement vs SABRE (positive = better than SABRE)")
        print(header)
        print("-" * len(header))
        for method in others:
            cells = []
            for key, _ in metrics:
                pct = [
                    (sabre_by_index[r["qc_index"]][key] - r[key])
                    / sabre_by_index[r["qc_index"]][key]
                    * 100
                    for r in by_method[method]
                    if sabre_by_index.get(r["qc_index"])
                    and sabre_by_index[r["qc_index"]][key]
                ]
                if pct:
                    std = stdev(pct) if len(pct) > 1 else 0.0
                    cells.append(f"{mean(pct):+.1f}% +/- {std:.1f}%")
                else:
                    cells.append("n/a")
            print(
                f"{method:<{name_w}}"
                + "".join(f"  {c:>{col_w}}" for c in cells)
            )

def print_per_circuit_comparison(results, num_rows=5):
    """
    Print a per-metric comparison of the compilation methods for the
    first ``num_rows`` circuits (sorted by qubit count). The best
    (lowest) value for each metric is marked with an asterisk.
    """
    metrics = [
        ("two_qubit_depth", "2Q Depth"),
        ("size", "Gate Count"),
        ("runtime", "Runtime (s)"),
    ]
    methods = _method_order(results)

    by_index = {}
    for r in results:
        by_index.setdefault(r["qc_index"], {})[r["method"]] = r
    ordered = sorted(
        by_index.items(),
        key=lambda kv: (next(iter(kv[1].values()))["num_qubits"], kv[0]),
    )[:num_rows]

    for key, label in metrics:
        print(f"{label} (first {num_rows} circuits by qubit count); * = best")
        header = f"{'Idx':>3}  {'Circuit':<16} {'Q':>3}" + "".join(
            f"{m:>9}" for m in methods
        )
        print(header)
        print("-" * len(header))
        for idx, method_map in ordered:
            any_record = next(iter(method_map.values()))
            present = {
                m: method_map[m][key] for m in methods if m in method_map
            }
            best = min(present.values())
            line = (
                f"{idx:>3}  {any_record['qc_name'][:16]:<16} "
                f"{any_record['num_qubits']:>3}"
            )
            for m in methods:
                value = method_map[m][key]
                text = f"{value:.2f}" if key == "runtime" else f"{int(value)}"
                if value == best:
                    text += "*"
                line += f"{text:>9}"
            print(line)
        print()

# Obtain the Hamiltonian JSON from the benchpress repository
url = "https://raw.githubusercontent.com/Qiskit/benchpress/e7b29ef7be4cc0d70237b8fdc03edbd698908eff/benchpress/hamiltonian/hamlib/100_representative.json"
response = requests.get(url)
response.raise_for_status()
ham_records = json.loads(response.text)

# Remove circuits that are too large for the backend
ham_records = [
    h for h in ham_records if h["ham_qubits"] <= backend.num_qubits
]

# Build PauliEvolutionGate circuits
qc_ham_list = []
for h in ham_records:
    terms = h["ham_hamlib_hamiltonian_terms"]
    coeff = h["ham_hamlib_hamiltonian_coefficients"]
    num_qubits = h["ham_qubits"]
    name = h["ham_problem"]

    evo_gate = PauliEvolutionGate(SparsePauliOp(terms, coeff))
    qc = QuantumCircuit(num_qubits)
    qc.name = name
    qc.append(evo_gate, range(num_qubits))
    qc_ham_list.append(qc)

# Remove circuits whose decomposed size exceeds 1500 gates so that transpilation completes in a reasonable time frame
qc_ham_list = [qc for qc in qc_ham_list if qc.decompose().size() <= 1500]

print(f"Total Hamiltonian circuits loaded: {len(qc_ham_list)}")
print(
    f"Qubit range: {min(qc.num_qubits for qc in qc_ham_list)} to {max(qc.num_qubits for qc in qc_ham_list)}"
)

qc_small = [qc for qc in qc_ham_list if qc.num_qubits < 20]
qc_large = [qc for qc in qc_ham_list if qc.num_qubits >= 20]

print(f"Small-scale circuits (<20 qubits): {len(qc_small)}")
print(f"Large-scale circuits (>=20 qubits): {len(qc_large)}")

# We decompose the circuit here, otherwise it would just be a PauliEvolutionGate box,
# which isn't very informative to look at!
qc_small[0].decompose().draw("mpl", fold=-1)

results_small = []

tqc_sabre_small = capture_transpilation_metrics(
    results_small, pm_sabre, qc_small, "SABRE"
)
tqc_ai_small = capture_transpilation_metrics(
    results_small, pm_ai, qc_small, "AI"
)
tqc_rustiq_small = capture_transpilation_metrics(
    results_small, pm_rustiq, qc_small, "Rustiq"
)

print_summary_table(results_small)

print_per_circuit_comparison(results_small, num_rows=8)

def plot_transpilation_comparison(results, title_prefix):
    """
    Create a three-panel figure comparing compilation methods on
    two-qubit depth, circuit size, and runtime.

    Circuits are sorted by qubit count and plotted by circuit index.
    """
    methods = _method_order(results)
    palette = {"SABRE": "#1f77b4", "AI": "#ff7f0e", "Rustiq": "#2ca02c"}
    markers = {"SABRE": "o", "AI": "^", "Rustiq": "s"}

    # Order circuits by qubit count (then index) and map to plot positions
    ref = sorted(
        [r for r in results if r["method"] == methods[0]],
        key=lambda r: (r["num_qubits"], r["qc_index"]),
    )
    pos_map = {r["qc_index"]: pos for pos, r in enumerate(ref)}
    tick_positions = [pos_map[r["qc_index"]] for r in ref]
    tick_labels = [
        f"{pos_map[r['qc_index']]} ({r['num_qubits']}q)" for r in ref
    ]

    metrics = [
        ("two_qubit_depth", "Two-Qubit Depth"),
        ("size", "Total Gate Count (Circuit Size)"),
        ("runtime", "Transpilation Runtime (s)"),
    ]

    fig, axes = plt.subplots(1, 3, figsize=(20, 5.5))
    fig.suptitle(title_prefix, fontsize=15, fontweight="bold", y=1.02)

    for ax, (metric, ylabel) in zip(axes, metrics):
        for method in methods:
            subset = sorted(
                [r for r in results if r["method"] == method],
                key=lambda r: pos_map[r["qc_index"]],
            )
            ax.plot(
                [pos_map[r["qc_index"]] for r in subset],
                [r[metric] for r in subset],
                marker=markers.get(method, "o"),
                label=method,
                color=palette.get(method, None),
                linewidth=1.5,
                markersize=6,
                alpha=0.85,
            )
        ax.set_xlabel("Circuit Index (num qubits)", fontsize=11)
        ax.set_ylabel(ylabel, fontsize=11)
        ax.legend(frameon=True, fontsize=9)
        ax.grid(True, linestyle="--", alpha=0.4)
        step = max(1, len(tick_positions) // 15)
        ax.set_xticks(tick_positions[::step])
        ax.set_xticklabels(
            [tick_labels[i] for i in range(0, len(tick_labels), step)],
            fontsize=7,
            rotation=45,
            ha="right",
        )

    plt.tight_layout()
    plt.show()

def plot_pct_improvement_vs_sabre(results, title_prefix):
    """
    Plot the per-circuit percent improvement of each non-SABRE method
    relative to SABRE, for each metric. A positive value means the
    method improved on SABRE; negative means SABRE was better.
    """
    metrics = [
        ("two_qubit_depth", "2Q Depth"),
        ("size", "Gate Count"),
        ("runtime", "Runtime"),
    ]
    palette = {"AI": "#ff7f0e", "Rustiq": "#2ca02c"}
    markers = {"AI": "^", "Rustiq": "s"}

    methods = _method_order(results)
    sabre = sorted(
        [r for r in results if r["method"] == "SABRE"],
        key=lambda r: (r["num_qubits"], r["qc_index"]),
    )
    other_methods = [m for m in methods if m != "SABRE"]

    tick_positions = list(range(len(sabre)))
    tick_labels = [
        f"{i} ({sabre[i]['num_qubits']}q)" for i in range(len(sabre))
    ]

    fig, axes = plt.subplots(1, 3, figsize=(20, 5.5))
    fig.suptitle(
        f"{title_prefix}: % Improvement over SABRE",
        fontsize=15,
        fontweight="bold",
        y=1.02,
    )

    for ax, (metric, label) in zip(axes, metrics):
        ax.axhline(0, color="#1f77b4", linewidth=2, label="SABRE (baseline)")
        for method in other_methods:
            data = sorted(
                [r for r in results if r["method"] == method],
                key=lambda r: (r["num_qubits"], r["qc_index"]),
            )
            pct = [
                (sabre[i][metric] - data[i][metric]) / sabre[i][metric] * 100
                for i in range(len(sabre))
            ]
            ax.plot(
                tick_positions,
                pct,
                marker=markers.get(method, "o"),
                label=method,
                color=palette.get(method, None),
                linewidth=1.5,
                markersize=6,
                alpha=0.85,
            )
        ax.set_xlabel("Circuit Index (num qubits)", fontsize=11)
        ax.set_ylabel(f"% Improvement ({label})", fontsize=11)
        ax.legend(frameon=True, fontsize=9)
        ax.grid(True, linestyle="--", alpha=0.4)
        step = max(1, len(tick_positions) // 15)
        ax.set_xticks(tick_positions[::step])
        ax.set_xticklabels(
            [tick_labels[i] for i in range(0, len(tick_labels), step)],
            fontsize=7,
            rotation=45,
            ha="right",
        )
        ylims = ax.get_ylim()
        ax.axhspan(0, max(ylims[1], 1), alpha=0.04, color="green")
        ax.axhspan(min(ylims[0], -1), 0, alpha=0.04, color="red")

    plt.tight_layout()
    plt.show()

plot_transpilation_comparison(
    results_small,
    "Small-Scale Hamiltonian Circuits: Compilation Comparison",
)

plot_pct_improvement_vs_sabre(
    results_small,
    "Small-Scale Hamiltonian Circuits",
)

def plot_best_method_bars(results, metrics_list=None):
    """
    Plot a grouped bar chart showing the percentage of circuits
    where each method achieved the best (lowest) value for each metric.

    Ties are counted for all tied methods, so percentages per metric
    can sum to more than 100%.
    """
    if metrics_list is None:
        metrics_list = ["two_qubit_depth", "size", "runtime"]

    labels = {
        "two_qubit_depth": "2Q Depth",
        "size": "Gate Count",
        "runtime": "Runtime",
    }
    methods = _method_order(results)
    palette = {"SABRE": "#1f77b4", "AI": "#ff7f0e", "Rustiq": "#2ca02c"}

    by_index = {}
    for r in results:
        by_index.setdefault(r["qc_index"], []).append(r)
    n_circuits = len(by_index)

    win_data = {m: [] for m in methods}
    tie_counts = []
    metric_labels = []

    for metric in metrics_list:
        metric_labels.append(
            labels.get(metric, metric.replace("_", " ").title())
        )
        counts = Counter()
        ties = 0
        for group in by_index.values():
            min_val = min(r[metric] for r in group)
            best = [r["method"] for r in group if r[metric] == min_val]
            if len(best) > 1:
                ties += 1
            counts.update(best)
        tie_counts.append(ties)
        for m in methods:
            win_data[m].append(counts.get(m, 0) / n_circuits * 100)

    x = np.arange(len(metric_labels))
    width = 0.22
    fig, ax = plt.subplots(figsize=(8, 5))

    for i, method in enumerate(methods):
        bars = ax.bar(
            x + i * width,
            win_data[method],
            width,
            label=method,
            color=palette.get(method, None),
            edgecolor="black",
            linewidth=0.5,
        )
        for bar in bars:
            height = bar.get_height()
            if height > 0:
                ax.text(
                    bar.get_x() + bar.get_width() / 2,
                    height + 1.5,
                    f"{height:.0f}%",
                    ha="center",
                    va="bottom",
                    fontsize=9,
                )

    # Annotate tie counts below each metric label
    for j, ties in enumerate(tie_counts):
        if ties > 0:
            ax.text(
                x[j] + width,
                -8,
                f"({ties} tie{'s' if ties != 1 else ''})",
                ha="center",
                va="top",
                fontsize=8,
                color="gray",
            )

    ax.set_xticks(x + width)
    ax.set_xticklabels(metric_labels, fontsize=11)
    ax.set_ylabel("Circuits with best value (%)", fontsize=11)
    ax.set_title(
        "Best-Performing Method by Metric (ties counted for all tied methods)",
        fontsize=12,
        fontweight="bold",
    )
    ax.legend(frameon=True, fontsize=10)
    ax.set_ylim(-12, 120)
    ax.yaxis.set_major_formatter(ticker.PercentFormatter())
    ax.grid(axis="y", linestyle="--", alpha=0.4)

    plt.tight_layout()
    plt.show()

plot_best_method_bars(results_small)

# Select circuit index 6 from the small-scale transpiled circuits
test_idx = 6
test_circuit = qc_small[test_idx]
print(f"Test circuit: {test_circuit.name}, {test_circuit.num_qubits} qubits")

# Get the transpiled versions
tqc_methods_small = {
    "SABRE": tqc_sabre_small[test_idx],
    "AI": tqc_ai_small[test_idx],
    "Rustiq": tqc_rustiq_small[test_idx],
}

# Show transpilation metrics for this circuit
print(f"\nTranspilation metrics for circuit index {test_idx}:")
for method, tqc in tqc_methods_small.items():
    depth_2q = tqc.depth(lambda x: x.operation.num_qubits == 2)
    size = tqc.size()
    print(f"  {method:8s}  2Q depth={depth_2q:5d}  size={size:6d}")

def remap_to_contiguous(tqc):
    """Remap a transpiled circuit to use contiguous qubit indices.

    Transpiled circuits target specific physical qubits (e.g., qubit 45, 67)
    on a large backend. This remaps them to 0, 1, 2, ... so Aer only
    simulates the active qubits.
    """
    active = sorted(
        {tqc.find_bit(q).index for inst in tqc.data for q in inst.qubits}
    )
    qubit_map = {old: new for new, old in enumerate(active)}
    new_qc = QuantumCircuit(len(active))
    for inst in tqc.data:
        old_indices = [tqc.find_bit(q).index for q in inst.qubits]
        new_qc.append(inst.operation, [qubit_map[i] for i in old_indices])
    return new_qc


def build_mirror_circuit(tqc):
    """Build a mirror circuit: U followed by U-dagger, with measurements.

    The combined circuit U-dagger @ U should be the identity, so measuring
    all zeros indicates a noise-free execution.
    """
    tqc_compact = remap_to_contiguous(tqc)
    mirror = tqc_compact.compose(tqc_compact.inverse())
    mirror.measure_all()
    return mirror


# Build a simple depolarizing noise model
noise_model = NoiseModel()
noise_model.add_all_qubit_quantum_error(
    depolarizing_error(0.001, 1),
    ["sx", "x", "rz"],  # ~0.1% per 1Q gate
)
noise_model.add_all_qubit_quantum_error(
    depolarizing_error(0.01, 2),
    ["cx", "ecr"],  # ~1% per 2Q gate
)

aer_sim = AerSimulator(noise_model=noise_model)

shots = 10000
fidelities = {}

for method, tqc in tqc_methods_small.items():
    mirror = build_mirror_circuit(tqc)

    sampler = SamplerV2(mode=aer_sim)
    job = sampler.run([mirror], shots=shots)
    result = job.result()
    counts = result[0].data.meas.get_counts()

    # Fidelity = fraction of all-zeros (error-free) outcomes
    n_qubits = mirror.num_qubits - mirror.num_clbits  # active qubits
    all_zeros = "0" * mirror.num_qubits
    fidelity = counts.get(all_zeros, 0) / shots
    fidelities[method] = fidelity
    print(
        f"{method:8s}  P(|00...0>) = {fidelity:.4f}  ({counts.get(all_zeros, 0)}/{shots})"
    )

def plot_mirror_results(tqc_methods, fidelities, circuit_name):
    """
    Plot a three-panel comparison: fidelity, 2Q depth,
    and gate count for each compilation method.
    """
    methods = list(tqc_methods.keys())
    palette = {"SABRE": "#1f77b4", "AI": "#ff7f0e", "Rustiq": "#2ca02c"}
    colors = [palette.get(m, "gray") for m in methods]

    fidelity_vals = [fidelities[m] for m in methods]
    depth_vals = [
        tqc_methods[m].depth(lambda x: x.operation.num_qubits == 2)
        for m in methods
    ]
    size_vals = [tqc_methods[m].size() for m in methods]

    fig, axes = plt.subplots(1, 3, figsize=(16, 5))
    fig.suptitle(
        f"Mirror Circuit Results: {circuit_name}",
        fontsize=14,
        fontweight="bold",
        y=1.02,
    )

    def _annotate_bars(ax, bars, values, fmt="{}"):
        ymax = ax.get_ylim()[1]
        for bar, val in zip(bars, values):
            label = fmt.format(val)
            y = val + ymax * 0.03
            ax.text(
                bar.get_x() + bar.get_width() / 2,
                y,
                label,
                ha="center",
                va="bottom",
                fontsize=10,
                fontweight="bold",
            )

    # Panel 1: Survival Probability
    bars = axes[0].bar(
        methods, fidelity_vals, color=colors, edgecolor="black", linewidth=0.5
    )
    axes[0].set_ylabel("Fidelity  P(|00...0>)", fontsize=11)
    axes[0].set_title("Fidelity (higher is better)", fontsize=12)
    axes[0].set_ylim(
        0, max(fidelity_vals) * 1.18 if max(fidelity_vals) > 0 else 1.0
    )
    axes[0].grid(axis="y", linestyle="--", alpha=0.4)
    _annotate_bars(axes[0], bars, fidelity_vals, fmt="{:.4f}")

    # Panel 2: Two-Qubit Depth
    bars = axes[1].bar(
        methods, depth_vals, color=colors, edgecolor="black", linewidth=0.5
    )
    axes[1].set_ylabel("Two-Qubit Depth", fontsize=11)
    axes[1].set_title("2Q Depth (lower is better)", fontsize=12)
    axes[1].set_ylim(0, max(depth_vals) * 1.18)
    axes[1].grid(axis="y", linestyle="--", alpha=0.4)
    _annotate_bars(axes[1], bars, depth_vals)

    # Panel 3: Gate Count
    bars = axes[2].bar(
        methods, size_vals, color=colors, edgecolor="black", linewidth=0.5
    )
    axes[2].set_ylabel("Total Gate Count", fontsize=11)
    axes[2].set_title("Gate Count (lower is better)", fontsize=12)
    axes[2].set_ylim(0, max(size_vals) * 1.18)
    axes[2].grid(axis="y", linestyle="--", alpha=0.4)
    _annotate_bars(axes[2], bars, size_vals)

    plt.tight_layout()
    plt.show()


plot_mirror_results(tqc_methods_small, fidelities, test_circuit.name)

results_large = []

tqc_sabre_large = capture_transpilation_metrics(
    results_large, pm_sabre, qc_large, "SABRE"
)
tqc_ai_large = capture_transpilation_metrics(
    results_large, pm_ai, qc_large, "AI"
)
tqc_rustiq_large = capture_transpilation_metrics(
    results_large, pm_rustiq, qc_large, "Rustiq"
)

print_summary_table(results_large)

print_per_circuit_comparison(results_large, num_rows=8)

plot_transpilation_comparison(
    results_large,
    "Large-Scale Hamiltonian Circuits: Compilation Comparison",
)

plot_pct_improvement_vs_sabre(
    results_large,
    "Large-Scale Hamiltonian Circuits",
)

plot_best_method_bars(results_large)

# Select circuit index 3 from the large-scale transpiled circuits
test_idx_large = 3
test_circuit_large = qc_large[test_idx_large]
print(
    f"Test circuit: {test_circuit_large.name}, {test_circuit_large.num_qubits} qubits"
)

tqc_methods_large = {
    "SABRE": tqc_sabre_large[test_idx_large],
    "AI": tqc_ai_large[test_idx_large],
    "Rustiq": tqc_rustiq_large[test_idx_large],
}

print(f"\nTranspilation metrics for circuit index {test_idx_large}:")
for method, tqc in tqc_methods_large.items():
    depth_2q = tqc.depth(lambda x: x.operation.num_qubits == 2)
    size = tqc.size()
    print(f"  {method:8s}  2Q depth={depth_2q:5d}  size={size:6d}")

pm_mirror = generate_preset_pass_manager(
    optimization_level=0, backend=backend
)

for method, tqc in tqc_methods_large.items():
    # print the count ops for each circuit
    mirror = tqc.copy()
    mirror.compose(tqc.inverse(), inplace=True)
    mirror.measure_all()
    mirror = pm_mirror.run(mirror)
    print(f"\n{method} transpiled circuit:")
    print(tqc.count_ops())
    print(f"{method} mirror circuit count ops:")
    print(mirror.count_ops())

# Build mirror circuits and submit to real hardware
# The inverse may introduce gates (e.g., sxdg) not in the backend's
# basis gate set, so we re-transpile the mirror circuit.
pm_mirror = generate_preset_pass_manager(
    optimization_level=0, backend=backend
)

shots_hw = 10000
hw_jobs = {}

for method, tqc in tqc_methods_large.items():
    mirror = tqc.copy()
    mirror.compose(tqc.inverse(), inplace=True)
    mirror.measure_all()

    # Re-transpile at opt level 0 to decompose into basis gates
    # without changing the layout or routing
    mirror = pm_mirror.run(mirror)

    sampler = SamplerV2(mode=backend)
    sampler.options.environment.job_tags = ["TUT_CMHSC"]
    job = sampler.run([mirror], shots=shots_hw)
    hw_jobs[method] = job
    print(f"{method}: submitted job {job.job_id()}")

# Retrieve results and compute fidelities
fidelities_large = {}

for method, job in hw_jobs.items():
    result = job.result()
    counts = result[0].data.meas.get_counts()

    n_qubits = backend.num_qubits
    all_zeros = "0" * n_qubits
    fidelity = counts.get(all_zeros, 0) / shots_hw
    fidelities_large[method] = fidelity
    print(
        f"{method:8s}  P(|00...0>) = {fidelity:.4f}  ({counts.get(all_zeros, 0)}/{shots_hw})"
    )

plot_mirror_results(
    tqc_methods_large, fidelities_large, test_circuit_large.name
)