Learning MCMs

Utility functions and Environment setup

• Joint instrument per shot:

  \[ p_s^{(o,s')} :=\text{Prob}(\text{outcome}=o,\ \text{post-meas state}=s'\mid \text{pre-state}=s), \quad s,o,s'\in\{0,1\}, \]

  with column normalization \(\sum_{o,s'}p_s^{(o,s')}=1\) for each \(s\) (⇒ 6 Degrees of Freedom in total).

• Outcome-indexed operators (“observable operators”)

  \[ M^{(o)} \in \mathbb{R}_{\ge 0}^{2\times 2},\quad [M^{(o)}]_{s',s}=p_s^{(o,s')}. \]

  Note: columns of \(M^{(o)}\) sum to \(A^{(o)}_s=\sum_{s'}p_s^{(o,s')}\) (not 1). The two \(M\) matrices explicitly are:

  \[ M^{(o=0)} = \begin{pmatrix} p_0^{(0, 0)} & p_1^{(0, 0)} \\ p_0^{(0, 1)} & p_1^{(0, 1)} \end{pmatrix} \text{ , } M^{(o=1)} = \begin{pmatrix} p_0^{(1, 0)} & p_1^{(1, 0)} \\ p_0^{(1, 1)} & p_1^{(1, 1)} \end{pmatrix} \]

import matplotlib.pyplot as plt

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import numpy as np
import numpy.linalg as LA
from numpy.linalg import eig, inv, solve
import pandas as pd
from scipy.optimize import minimize, Bounds, LinearConstraint, basinhopping, differential_evolution, NonlinearConstraint
import os
from typing import Callable, Dict, Any, Tuple, List
import pickle
from itertools import product
from dataclasses import dataclass

from collections import Counter
from tqdm.notebook import tqdm
from IPython.display import display

rng = np.random.default_rng(7)


@dataclass
class Instrument2x2:
    """
    Holds the two symbol-operator matrices M0, M1 for a single-qubit instrument (mid-circuit measurement gadget).

    Each matrix is 2x2, with nonnegative entries, and the sum of the two matrices has columns summing to 1.

    M0 corresponds to measurement outcome 0, M1 to outcome 1.

    The entry M^o_[s',s] = p_s^{(o,s')} is the probability of obtaining outcome o and post-measurement state s'
    given starting state s.
    """

    M0: np.ndarray  # shape (2,2), nonnegative
    M1: np.ndarray  # shape (2,2), nonnegative

    def check_valid(self, atol: float = 1e-8) -> bool:
        """Validate nonnegativity and column-sum constraints."""
        assert self.M0.shape == (2, 2) and self.M1.shape == (2, 2)
        if (self.M0 < 0).any() or (self.M1 < 0).any():
            raise ValueError("Negative entries in M0/M1.")
        colsums = self.M0.sum(axis=0) + self.M1.sum(axis=0)
        if not np.allclose(colsums, np.ones(2), atol=atol):
            raise ValueError(f"Column sums of (M0+M1) must be 1, got {colsums}.")
        return True

    def reveal(self) -> int:
        """Print the instrument matrices in a readable format."""
        print("\nMCM Instrument:")
        print("\nM0 matrix (outcome 0):")
        print(
            f"  [[p_0^(0,0), p_1^(0,0)]] = [[{self.M0[0, 0]:.8f}, {self.M0[0, 1]:.8f}]]"
        )
        print(
            f"  [[p_0^(0,1), p_1^(0,1)]]   [[{self.M0[1, 0]:.8f}, {self.M0[1, 1]:.8f}]]"
        )
        print("\nM1 matrix (outcome 1):")
        print(
            f"  [[p_0^(1,0), p_1^(1,0)]] = [[{self.M1[0, 0]:.8f}, {self.M1[0, 1]:.8f}]]"
        )
        print(
            f"  [[p_0^(1,1), p_1^(1,1)]]   [[{self.M1[1, 0]:.8f}, {self.M1[1, 1]:.8f}]]"
        )
        return 0

    def show_readout_errors(self) -> List[float]:
        """
        Calculate and return the readout error rates for starting states |0> and |1>.
        Returns:
            List of readout error rates [error_rate_for_|0>, error_rate_for_|1>]
        """
        error_rate_0 = self.M1[0, 0] + self.M1[1, 0]  # Probability of getting outcome 1 when starting from |0>
        error_rate_1 = self.M0[0, 1] + self.M0[1, 1]  # Probability of getting outcome 0 when starting from |1>

        print(f"\nPrep 0 meas 1: {error_rate_0:.3e}")
        print(f"Prep 1 meas 0: {error_rate_1:.3e}")

        return [error_rate_0, error_rate_1]

    def show_backaction_errors(self) -> List[float]:
        """
        Calculate and return the back-action error rates for starting states |0> and |1>.
        Returns:
            List of back-action error rates [backaction_error_for_|0>, backaction_error_for_|1>]
        """
        backaction_error_0 = self.M0[1, 0] + self.M1[1, 0]  # Probability of ending in state |1> when starting from |0>
        backaction_error_1 = self.M0[0, 1] + self.M1[0, 1]  # Probability of ending in state |0> when starting from |1>

        print(f"\nPrep 0 ends in 1: {backaction_error_0:.3e}")
        print(f"Prep 1 ends in 0: {backaction_error_1:.3e}")

        return [backaction_error_0, backaction_error_1]


def make_instrument_from_columns(
    col0: np.ndarray,
    col1: np.ndarray
) -> Instrument2x2:
    """
    Build Instrument2x2 from two columns of 4 probabilities each:
    col_s = [p_s^{(0,0)}, p_s^{(0,1)}, p_s^{(1,0)}, p_s^{(1,1)}], sums to 1.
    Args:
        col0: Column for starting state |0>
        col1: Column for starting state |1>
    Returns:
        Instrument2x2 instance
    """
    col0 = np.asarray(col0, dtype=float).reshape(4)
    col1 = np.asarray(col1, dtype=float).reshape(4)
    if not np.isclose(col0.sum(), 1.0) or not np.isclose(col1.sum(), 1.0):
        raise ValueError("Each column-vector must sum to 1.")
    if (col0 < 0).any() or (col1 < 0).any():
        raise ValueError("Probabilities must be nonnegative.")

    M0 = np.zeros((2, 2))
    M1 = np.zeros((2, 2))

    # Column s=0, starting state |0>
    M0[0, 0], M0[1, 0], M1[0, 0], M1[1, 0] = col0
    # Column s=1, starting state |1>
    M0[0, 1], M0[1, 1], M1[0, 1], M1[1, 1] = col1

    instrument = Instrument2x2(M0=M0, M1=M1)
    instrument.check_valid()
    return instrument


def random_instrument(
    corr_strength: float = 0.2, 
    seed: int = 1, 
    fidelity: float = 0.85
) -> Instrument2x2:
    """
    Sample a random but 'physical' instrument with a tunable correlation flavor.
    We bias columns so that outcome=0 is more likely when s=0 and vice versa, and
    inject some back-action to create realistic correlations.
    For high-fidelity cases, p_0^(0,0) and p_1^(1,1) should be close to 1.

    Args:
        corr_strength: Tunes how much the Dirichlet distribution is biased.
        seed: Random seed for reproducibility.
        fidelity: The base probability for the 'correct' outcome (e.g., p_0(0,0)). Must be between 0.25 and 1.0.
    Returns:
        Instrument2x2 instance
    """
    if not (0.25 <= fidelity <= 1.0):
        raise ValueError("Fidelity must be between 0.25 and 1.0.")

    local_rng = np.random.default_rng(seed) if seed is not None else rng

    def biased_dirichlet(s: int) -> np.ndarray:
        """
        Sample a column of the instrument with bias towards correct outcomes.
        Args:
            s: The input state index (0 or 1).
        Returns:
            A column vector of 4 probabilities summing to 1.
        """
        # Base weights for [p_s(0,0), p_s(0,1), p_s(1,0), p_s(1,1)]
        off_diag_prob = (1.0 - fidelity) / 3.0
        if s == 0:  # Input state |0>
            # Expect outcome 0, state |0> (p_0(0,0) is high)
            base = np.array([fidelity, off_diag_prob, off_diag_prob, off_diag_prob])
        else:  # Input state |1>
            # Expect outcome 1, state |1> (p_1(1,1) is high)
            base = np.array([off_diag_prob, off_diag_prob, off_diag_prob, fidelity])

        bias = base ** (1.0 - corr_strength)
        return local_rng.dirichlet(1e-3 + 20 * bias)

    col0 = biased_dirichlet(0)
    col1 = biased_dirichlet(1)
    return make_instrument_from_columns(col0, col1)


def get_all_binary_strings(
    max_len: int
) -> List[str]:
    """
    Generates all binary strings up to a given maximum length.
    Args:
        max_len: Maximum length of binary strings to generate.
    Returns:
        List of binary strings.
    Example:
        get_all_binary_strings(2) returns ['0', '1', '00', '01', '10', '11']
    """
    strings = []
    for length in range(1, max_len + 1):
        for s_tuple in product("01", repeat=length):
            strings.append("".join(s_tuple))
    return strings


def calculate_exact_all_string_probabilities_from_v0_and_instrument(
    inst: Instrument2x2, 
    v0: np.ndarray, 
    max_len: int
) -> Dict[str, float]:
    """
    Calculate exact probabilities for all binary strings up to max_len (L)
    given an initial state distribution v0 and a MCM instrument.
    Exact Prob(w) = [1,1]^T M^{(w_L)} ... M^{(w_1)} v0 for all |w|<=L.
    Args:
        inst: Instrument2x2 instance representing the MCM instrument.
        v0: Initial state distribution as a numpy array of shape (2,).
        max_len: Maximum length of binary strings to consider.
    Returns:
        Dictionary mapping binary strings to their exact probabilities.
    """
    all_strings = get_all_binary_strings(max_len)
    all_probs = {}

    for single_string in all_strings:
        M_product = np.eye(2)
        for outcome_char in single_string:
            M_product = (inst.M0 if outcome_char == "0" else inst.M1) @ M_product
        all_probs[single_string] = np.sum(M_product @ v0)

    return all_probs


def monte_carlo_empirical_probabilities(
    inst: Instrument2x2, 
    v0: np.ndarray, 
    L: int, 
    shots: int, 
    seed: int = 1234
) -> Dict[str, float]:
    """
    Simulate the MCM process to obtain empirical probabilities for all binary strings up to length L.
    Args:
        inst: Instrument2x2 instance representing the MCM instrument.
        v0: Initial state distribution as a numpy array of shape (2,).
        L: Maximum length of binary strings to observe.
        shots: Number of Monte Carlo shots to simulate.
        seed: Random seed for reproducibility.
    Returns:
        Dictionary mapping binary strings to their empirical probabilities.
    """
    rng_local = np.random.default_rng(seed) if seed is not None else rng
    counts = Counter()
    total_prefixes = Counter({ell: 0 for ell in range(1, L + 1)})

    # Precompute per-column cumulative distributions over (o,s')
    cum_cdfs = {}
    for s in [0, 1]:
        c = np.array([inst.M0[0, s], inst.M0[1, s], inst.M1[0, s], inst.M1[1, s]])
        cum_cdfs[s] = np.cumsum(c)

    for _ in tqdm(range(shots), desc="Monte Carlo on strings"):

        # Z-twirl: flip v_initial with 50% probability
        v_after_twirling = v0 if rng_local.random() < 0.5 else v0[::-1]

        # sample the initial state from {0,1} according to v_after_twirling
        s = 0 if rng_local.random() < v_after_twirling[0] else 1

        # sample L steps of the MCM process
        readout_sequence = []
        for ell in range(1, L + 1):
            r = rng_local.random()
            cdf = cum_cdfs[s]
            idx = int(np.searchsorted(cdf, r, side="right"))
            if idx == 0:
                o, post_measurement_state = 0, 0
            elif idx == 1:
                o, post_measurement_state = 0, 1
            elif idx == 2:
                o, post_measurement_state = 1, 0
            else:
                o, post_measurement_state = 1, 1
            readout_sequence.append(str(o))
            s = post_measurement_state
            w = "".join(readout_sequence)
            counts[w] += 1
            total_prefixes[ell] += 1

    empirical_probabilities = {w: c / total_prefixes[len(w)] for w, c in counts.items()}
    return empirical_probabilities


def generate_and_cache_simulated_data(
    GT: Instrument2x2,
    num_seeds: int,
    shots_per_seed: int = 100_000,
    L: int = 4,
    use_same_v0_for_all_seeds: None | np.ndarray = np.array([0.5, 0.5]),
    seed_for_reproduce: int = 123,
    chosen_dir: str = "simulated_data"
) -> Tuple[Dict[str, Any], List[np.ndarray], List[Dict[str, float]]]:
    """
    Generate and cache simulated data for multiple seeds.
    Args:
        GT: Ground truth Instrument2x2 instance.
        num_seeds: Number of different random seeds to simulate.
        shots_per_seed: Number of Monte Carlo shots per seed.
        L: Maximum length of binary strings to simulate.
        use_same_v0_for_all_seeds: If provided, use this initial state distribution for all seeds.
        seed_for_reproduce: Seed for reproducibility of the entire simulation.
        chosen_dir: Directory to save the simulated data.
    Returns:
        Tuple of:
            - Experiment info dictionary.
            - List of initial state distributions (v0) for each seed.
            - List of empirical probability dictionaries for each seed.
    """

    v0_list = []
    empirical_probs_list = []
    experiment_info = {
        "GT": GT,
        "number_of_seeds": num_seeds,
        "shots_per_seed": shots_per_seed,
        "max_string_length_L": L,
        "use_same_v0_for_all_seeds": use_same_v0_for_all_seeds,
        "seed_for_reproduce": seed_for_reproduce
    }

    path_to_info = "experiment_info.pkl"
    path_to_v0_and_emp_probs = "v0_list_and_emp_probs_list.pkl"

    if os.path.exists(chosen_dir):
        print(f"Loading cached simulated data from {chosen_dir}...")
        print("Override experiment info with cached version.")
        with open(os.path.join(chosen_dir, path_to_info), "rb") as f:
            experiment_info = pickle.load(f)
        with open(os.path.join(chosen_dir, path_to_v0_and_emp_probs), "rb") as f:
            v0_list, empirical_probs_list = pickle.load(f)
        return experiment_info, v0_list, empirical_probs_list
    else:
        pass

    os.makedirs(chosen_dir)
    main_rng = np.random.default_rng(seed_for_reproduce)
    v_initial = main_rng.dirichlet([1.0, 1.0]) if use_same_v0_for_all_seeds is None else use_same_v0_for_all_seeds

    for i in tqdm(range(num_seeds), desc="Generating data for seeds"):
        monte_carlo_seed = main_rng.integers(0, 2**32 - 1)

        v0_list.append(v_initial)

        emp_probs = monte_carlo_empirical_probabilities(
            inst=GT,
            v0=v0_list[-1],
            L=L,
            shots=shots_per_seed,
            seed=monte_carlo_seed
        )
        empirical_probs_list.append(emp_probs)

    with open(os.path.join(chosen_dir, path_to_info), "wb") as f:
        pickle.dump(experiment_info, f)

    with open(os.path.join(chosen_dir, path_to_v0_and_emp_probs), "wb") as f:
        pickle.dump((v0_list, empirical_probs_list), f)

    return experiment_info, v0_list, empirical_probs_list


def calculate_average_probs_np(
    empirical_probs_list
) -> Dict[str, float]:
    """
    Calculates the average probabilities for each key across a list of dictionaries.

    Args:
        empirical_probs_list (list[dict]): 
            A list of dictionaries, where each
            dictionary maps string keys to
            probability floats.

    Returns:
        dict: 
            A single dictionary mapping each string key to its
            average probability across all input dictionaries.
    """
    # Return an empty dict if the input list is empty
    if not empirical_probs_list:
        return {}

    # 1. Get a fixed, sorted list of all keys from the first dictionary
    ordered_keys = sorted(empirical_probs_list[0].keys())

    # 2. Create a 2D NumPy array from the values
    data_array = np.array([
        [prob_dict[key] for key in ordered_keys] 
        for prob_dict in empirical_probs_list
    ])

    # 3. Calculate the mean of each column (axis=0)
    average_values = np.mean(data_array, axis=0)

    # 4. Zip the keys back with the average values to create the final dict
    average_probs_np = dict(zip(ordered_keys, average_values))

    return average_probs_np


def derived_constraints_from_empirical_probs(
    emp_probs: Dict[str, float]
) -> Tuple[float, float, float, float, None | float, None | float, float, float]:
    """
    Derive the constraint values from empirical probabilities.
    Args:
        emp_probs: Dictionary of empirical probabilities for binary strings.
    Returns:
        Tuple of derived constraint values:
            - trace_M0
            - det_M0
            - trace_M1
            - det_M1
            - trace_M0M1 and trace_M1M0 (they are the same but we can use both to check consistency)
            - prob_string_0 (need assumption of v0 being in maximally-mixed state)
            - prob_string_1 (need assumption of v0 being in maximally-mixed state)
    """

    # trace_M0 = (Prob('0') * Prob('00') - Prob('000')) / (Prob('0')^2 - Prob('00'))
    trace_M0 = (emp_probs.get('0', 0) * emp_probs.get('00', 0) - emp_probs.get('000', 0)) / \
                (emp_probs.get('0', 0)**2 - emp_probs.get('00', 0))

    # det_M0 = trace_M0 * Prob('0') - Prob('00')
    det_M0 = trace_M0 * emp_probs.get('0', 0) - emp_probs.get('00', 0)

    # trace_M1 = (Prob('1') * Prob('11') - Prob('111')) / (Prob('1')^2 - Prob('11'))
    trace_M1 = (emp_probs.get('1', 0) * emp_probs.get('11', 0) - emp_probs.get('111', 0)) / \
                (emp_probs.get('1', 0)**2 - emp_probs.get('11', 0))

    # det_M1 = trace_M1 * Prob('1') - Prob('11')
    det_M1 = trace_M1 * emp_probs.get('1', 0) - emp_probs.get('11', 0)

    trace_M0M1 = None
    trace_M1M0 = None

    # trace_M0M1 = (Prob('0101') + det_M1 * det_M0) / Prob('01')
    if '0101' in emp_probs and '01' in emp_probs and emp_probs['01'] != 0:
        trace_M0M1 = (emp_probs['0101'] + det_M1 * det_M0) / emp_probs['01']

    # trace_M1M0 = (Prob('1010') + det_M0 * det_M1) / Prob('10')
    if '1010' in emp_probs and '10' in emp_probs and emp_probs['10'] != 0:
        trace_M1M0 = (emp_probs['1010'] + det_M0 * det_M1) / emp_probs['10']

    return (trace_M0, 
            det_M0, 
            trace_M1, 
            det_M1, 
            trace_M0M1,
            trace_M1M0,
            emp_probs.get('0', 0),
            emp_probs.get('1', 0))

🟢 Gauge Continuum Characterization

Symbolic and Numeric Gauge Transformation

import sympy as sp

def construct_instrument_symbolically(
    a: sp.Symbol,
    b: sp.Symbol,
    c: sp.Symbol,
    e: sp.Symbol,
    f: sp.Symbol,
    g: sp.Symbol
) -> Tuple[sp.Matrix, sp.Matrix]:
    """
    Constructs symbolic instrument matrices M0 and M1 from 6 independent parameters.
    The parameters are probabilities for each column, which sum to 1.
    p_s^{(o,s')} are the probabilities.
    col0 = [p_0^{(0,0)}, p_0^{(0,1)}, p_0^{(1,0)}, p_0^{(1,1)}]
    col1 = [p_1^{(0,0)}, p_1^{(0,1)}, p_1^{(1,0)}, p_1^{(1,1)}]

    Args:
        a, b, c: Parameters for column 0. The fourth, p1_10, is 1 - a - b - c.
        e, f, g: Parameters for column 1. The fourth, p1_11, is 1 - e - f - g.

    Returns:
        A tuple of symbolic matrices (M0, M1).
    """
    # Column 0:
    p0_00 = a
    p0_10 = b
    p1_00 = c
    p1_10 = 1 - a - b - c

    # Column 1:
    p0_01 = e
    p0_11 = f
    p1_01 = g
    p1_11 = 1 - e - f - g

    M0 = sp.Matrix([[p0_00, p0_01], [p0_10, p0_11]])
    M1 = sp.Matrix([[p1_00, p1_01], [p1_10, p1_11]])

    return M0, M1


def gauge_transform_instrument_symbolically(
    M0: sp.Matrix,
    M1: sp.Matrix,
    t: sp.Symbol,
    ) -> Tuple[sp.Matrix, sp.Matrix]:
    """
    Apply a symbolic gauge transformation to the instrument matrices M0 and M1.
    M' = R^{-1} @ M @ R
    where R is a 2x2 matrix [[1 - t, t], [t, 1 - t]]

    Args:
        M0: Symbolic matrix for M0.
        M1: Symbolic matrix for M1.
        t: Symbolic variable for the gauge parameter.
    Returns:
        Transformed matrices M0' and M1'.
    """
    R = sp.Matrix([[1 - t, t], [t, 1 - t]])
    R_inv = R.inv() # Symbolic inverse
    M0_hat = R_inv * M0 * R
    M1_hat = R_inv * M1 * R
    return M0_hat, M1_hat


def gauge_transform_instrument_numerically(
    M0: np.ndarray,
    M1: np.ndarray,
    t: float
    ) -> Tuple[np.ndarray, np.ndarray]:
    """
    Apply a numerical gauge transformation to the instrument matrices M0 and M1.
    M' = R^{-1} @ M @ R
    where R is a 2x2 matrix [[1 - t, t], [t, 1 - t]]

    Args:
        M0: Numpy array for M0.
        M1: Numpy array for M1.
        t: Float value for the gauge parameter in [0,1].
    Returns:
        Transformed matrices M0' and M1'.
    """
    R = np.array([[1 - t, t], [t, 1 - t]])
    R_inv = LA.inv(R)
    M0_hat = R_inv @ M0 @ R
    M1_hat = R_inv @ M1 @ R
    return M0_hat, M1_hat

# This cell is to analytically check the gauge transformation don't change invariants.

# Define symbolic variables with assumptions
a = sp.Symbol('a', real=True)
b = sp.Symbol('b', real=True)
c = sp.Symbol('c', real=True)
e = sp.Symbol('e', real=True)
f = sp.Symbol('f', real=True)
g = sp.Symbol('g', real=True)
t = sp.Symbol('t', real=True)

# Add assumptions about the sums and the range of t
# These are not added to the symbols directly but used in subsequent symbolic manipulations.
# sp.assume(a + b + c <= 1) # sp.assume is for new-style assumptions, let's handle this logically.
# sp.assume(e + f + g <= 1)
# sp.assume(t <= 1)

# Construct the symbolic instrument
test_M0, test_M1 = construct_instrument_symbolically(a, b, c, e, f, g)

# Apply the gauge transformation
test_m0_prime, test_m1_prime = gauge_transform_instrument_symbolically(
    test_M0, test_M1, t
)

# Simplify the expressions for the transformed matrices
simplified_m0_prime = sp.simplify(test_m0_prime)
simplified_m1_prime = sp.simplify(test_m1_prime)

print("Original M0 and M1:")
display(test_M0, test_M1)

print("Original 2 matrices' det and trace:")
print(f"det(M0) = {sp.simplify(test_M0.det())}, trace(M0) = {sp.simplify(test_M0.trace())}")
print(f"det(M1) = {sp.simplify(test_M1.det())}, trace(M1) = {sp.simplify(test_M1.trace())}")

print("The original trace of M0M1 and M1M0:")
print(f"trace(M0M1) = {sp.simplify((test_M0 * test_M1).trace())}")
print(f"trace(M1M0) = {sp.simplify((test_M1 * test_M0).trace())}")

print("\nTransformed M0' and M1' (before simplification):")
display(test_m0_prime, test_m1_prime)

print("\nTransformed M0' and M1' (after simplification):")
display(simplified_m0_prime, simplified_m1_prime)

print("\nTransformed 2 matrices' det and trace:")
print(f"det(M0') = {sp.simplify(simplified_m0_prime.det())}, trace(M0') = {sp.simplify(simplified_m0_prime.trace())}")
print(f"det(M1') = {sp.simplify(simplified_m1_prime.det())}, trace(M1') = {sp.simplify(simplified_m1_prime.trace())}")

print("The transformed trace of M0'P1' and P1'M0':")
print(f"trace(M0'P1') = {sp.simplify((simplified_m0_prime * simplified_m1_prime).trace())}")
print(f"trace(M1'P0') = {sp.simplify((simplified_m1_prime * simplified_m0_prime).trace())}")

# Check the detailed forms of each element
print("\nDetailed forms of each element in M0' and M1' after transformation:")
for i in range(2):
    for j in range(2):
        elem_m0_prime = sp.simplify(simplified_m0_prime[i, j])
        elem_m1_prime = sp.simplify(simplified_m1_prime[i, j])
        print(f"M0'[{i},{j}]: {elem_m0_prime}")
        print(f"M1'[{i},{j}]: {elem_m1_prime}")
        print("-"*40)

Original M0 and M1:

\(\displaystyle \left[\begin{matrix}a & e\\b & f\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}c & g\\- a - b - c + 1 & - e - f - g + 1\end{matrix}\right]\)

Original 2 matrices' det and trace:
det(M0) = a*f - b*e, trace(M0) = a + f
det(M1) = a*g + b*g - c*e - c*f + c - g, trace(M1) = c - e - f - g + 1
The original trace of M0M1 and M1M0:
trace(M0M1) = a*c + b*g - e*(a + b + c - 1) - f*(e + f + g - 1)
trace(M1M0) = a*c + b*g - e*(a + b + c - 1) - f*(e + f + g - 1)

Transformed M0' and M1' (before simplification):

\(\displaystyle \left[\begin{matrix}t \left(\frac{e \left(t - 1\right)}{2 t - 1} + \frac{f t}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{a \left(t - 1\right)}{2 t - 1} + \frac{b t}{2 t - 1}\right) & t \left(\frac{a \left(t - 1\right)}{2 t - 1} + \frac{b t}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{e \left(t - 1\right)}{2 t - 1} + \frac{f t}{2 t - 1}\right)\\t \left(\frac{e t}{2 t - 1} + \frac{f \left(t - 1\right)}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{a t}{2 t - 1} + \frac{b \left(t - 1\right)}{2 t - 1}\right) & t \left(\frac{a t}{2 t - 1} + \frac{b \left(t - 1\right)}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{e t}{2 t - 1} + \frac{f \left(t - 1\right)}{2 t - 1}\right)\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}t \left(\frac{g \left(t - 1\right)}{2 t - 1} + \frac{t \left(- e - f - g + 1\right)}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{c \left(t - 1\right)}{2 t - 1} + \frac{t \left(- a - b - c + 1\right)}{2 t - 1}\right) & t \left(\frac{c \left(t - 1\right)}{2 t - 1} + \frac{t \left(- a - b - c + 1\right)}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{g \left(t - 1\right)}{2 t - 1} + \frac{t \left(- e - f - g + 1\right)}{2 t - 1}\right)\\t \left(\frac{g t}{2 t - 1} + \frac{\left(t - 1\right) \left(- e - f - g + 1\right)}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{c t}{2 t - 1} + \frac{\left(t - 1\right) \left(- a - b - c + 1\right)}{2 t - 1}\right) & t \left(\frac{c t}{2 t - 1} + \frac{\left(t - 1\right) \left(- a - b - c + 1\right)}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{g t}{2 t - 1} + \frac{\left(t - 1\right) \left(- e - f - g + 1\right)}{2 t - 1}\right)\end{matrix}\right]\)

Transformed M0' and M1' (after simplification):

\(\displaystyle \left[\begin{matrix}\frac{t \left(e \left(t - 1\right) + f t\right) - \left(t - 1\right) \left(a \left(t - 1\right) + b t\right)}{2 t - 1} & \frac{t \left(a \left(t - 1\right) + b t\right) - \left(t - 1\right) \left(e \left(t - 1\right) + f t\right)}{2 t - 1}\\\frac{t \left(e t + f \left(t - 1\right)\right) - \left(t - 1\right) \left(a t + b \left(t - 1\right)\right)}{2 t - 1} & \frac{t \left(a t + b \left(t - 1\right)\right) - \left(t - 1\right) \left(e t + f \left(t - 1\right)\right)}{2 t - 1}\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}\frac{t \left(g \left(t - 1\right) - t \left(e + f + g - 1\right)\right) - \left(t - 1\right) \left(c \left(t - 1\right) - t \left(a + b + c - 1\right)\right)}{2 t - 1} & \frac{t \left(c \left(t - 1\right) - t \left(a + b + c - 1\right)\right) - \left(t - 1\right) \left(g \left(t - 1\right) - t \left(e + f + g - 1\right)\right)}{2 t - 1}\\\frac{t \left(g t - \left(t - 1\right) \left(e + f + g - 1\right)\right) - \left(t - 1\right) \left(c t - \left(t - 1\right) \left(a + b + c - 1\right)\right)}{2 t - 1} & \frac{t \left(c t - \left(t - 1\right) \left(a + b + c - 1\right)\right) - \left(t - 1\right) \left(g t - \left(t - 1\right) \left(e + f + g - 1\right)\right)}{2 t - 1}\end{matrix}\right]\)

Transformed 2 matrices' det and trace:
det(M0') = a*f - b*e, trace(M0') = a + f
det(M1') = a*g + b*g - c*e - c*f + c - g, trace(M1') = c - e - f - g + 1
The transformed trace of M0'P1' and P1'M0':
trace(M0'P1') = a*c - a*e - b*e + b*g - c*e - e*f + e - f**2 - f*g + f
trace(M1'P0') = a*c - a*e - b*e + b*g - c*e - e*f + e - f**2 - f*g + f

Detailed forms of each element in M0' and M1' after transformation:
M0'[0,0]: (t*(e*(t - 1) + f*t) - (t - 1)*(a*(t - 1) + b*t))/(2*t - 1)
M1'[0,0]: (t*(g*(t - 1) - t*(e + f + g - 1)) - (t - 1)*(c*(t - 1) - t*(a + b + c - 1)))/(2*t - 1)
----------------------------------------
M0'[0,1]: (t*(a*(t - 1) + b*t) - (t - 1)*(e*(t - 1) + f*t))/(2*t - 1)
M1'[0,1]: (t*(c*(t - 1) - t*(a + b + c - 1)) - (t - 1)*(g*(t - 1) - t*(e + f + g - 1)))/(2*t - 1)
----------------------------------------
M0'[1,0]: (t*(e*t + f*(t - 1)) - (t - 1)*(a*t + b*(t - 1)))/(2*t - 1)
M1'[1,0]: (t*(g*t - (t - 1)*(e + f + g - 1)) - (t - 1)*(c*t - (t - 1)*(a + b + c - 1)))/(2*t - 1)
----------------------------------------
M0'[1,1]: (t*(a*t + b*(t - 1)) - (t - 1)*(e*t + f*(t - 1)))/(2*t - 1)
M1'[1,1]: (t*(c*t - (t - 1)*(a + b + c - 1)) - (t - 1)*(g*t - (t - 1)*(e + f + g - 1)))/(2*t - 1)
----------------------------------------

Close to ibm_pittsburgh's MCM collection

ibm_pittsburgh_mcm = [
    random_instrument(corr_strength=0.005, seed=3, fidelity=0.98),
    random_instrument(corr_strength=0.08, seed=37, fidelity=0.96), # this might be a little little bit off
    random_instrument(corr_strength=0.2, seed=1, fidelity=0.99), # this is good for demo of a ideally near-perfect instrument
    random_instrument(corr_strength=0.1, seed=98, fidelity=0.97), # this is a good demo for ibm_pittsburgh in the morning of Nov 6, 2025.
]

# 😃 Find seeds that produce GT instruments with readout errors in the desired range
target_min = 5e-4
target_max = 4e-3

good_seeds = []

for seed in range(100):
    GT_test = random_instrument(corr_strength=0.1, seed=seed, fidelity=0.97)
    readout_errors = GT_test.show_readout_errors()

    # Check if both readout errors are within the target range
    if (target_min <= readout_errors[0] <= target_max and 
        target_min <= readout_errors[1] <= target_max):
        good_seeds.append({
            'seed': seed,
            'prep_0_meas_1': readout_errors[0],
            'prep_1_meas_0': readout_errors[1]
        })

print(f"\nFound {len(good_seeds)} seeds with readout errors in range [{target_min:.1e}, {target_max:.1e}]:\n")
print("="*70)

for item in good_seeds:
    print(f"Seed {item['seed']:3d}: prep_0_meas_1 = {item['prep_0_meas_1']:.6f}, "
          f"prep_1_meas_0 = {item['prep_1_meas_0']:.6f}")

# Create a summary DataFrame
if good_seeds:
    df_good_seeds = pd.DataFrame(good_seeds)
    print("\n" + "="*70)
    print("Summary Statistics:")
    print("="*70)
    display(df_good_seeds.describe())

Prep 0 meas 1: 3.954e-02
Prep 1 meas 0: 5.291e-02

Prep 0 meas 1: 3.216e-02
Prep 1 meas 0: 2.313e-02

Prep 0 meas 1: 9.910e-03
Prep 1 meas 0: 8.432e-03

Prep 0 meas 1: 3.295e-03
Prep 1 meas 0: 5.632e-03

Prep 0 meas 1: 8.134e-02
Prep 1 meas 0: 5.731e-02

Prep 0 meas 1: 3.745e-03
Prep 1 meas 0: 3.655e-03

Prep 0 meas 1: 1.185e-02
Prep 1 meas 0: 2.531e-02

Prep 0 meas 1: 1.157e-03
Prep 1 meas 0: 2.906e-02

Prep 0 meas 1: 6.697e-02
Prep 1 meas 0: 8.373e-04

Prep 0 meas 1: 6.891e-02
Prep 1 meas 0: 4.451e-03

Prep 0 meas 1: 4.473e-02
Prep 1 meas 0: 4.349e-02

Prep 0 meas 1: 9.686e-02
Prep 1 meas 0: 1.365e-02

Prep 0 meas 1: 1.654e-02
Prep 1 meas 0: 6.799e-02

Prep 0 meas 1: 4.431e-02
Prep 1 meas 0: 3.722e-02

Prep 0 meas 1: 2.976e-02
Prep 1 meas 0: 1.988e-02

Prep 0 meas 1: 1.738e-02
Prep 1 meas 0: 2.324e-03

Prep 0 meas 1: 5.305e-02
Prep 1 meas 0: 6.418e-04

Prep 0 meas 1: 3.145e-03
Prep 1 meas 0: 7.141e-04

Prep 0 meas 1: 6.877e-02
Prep 1 meas 0: 3.752e-02

Prep 0 meas 1: 2.885e-02
Prep 1 meas 0: 5.163e-02

Prep 0 meas 1: 3.447e-03
Prep 1 meas 0: 6.137e-03

Prep 0 meas 1: 1.067e-01
Prep 1 meas 0: 5.064e-02

Prep 0 meas 1: 6.313e-03
Prep 1 meas 0: 7.208e-03

Prep 0 meas 1: 1.226e-02
Prep 1 meas 0: 2.436e-02

Prep 0 meas 1: 4.061e-02
Prep 1 meas 0: 3.393e-02

Prep 0 meas 1: 6.340e-05
Prep 1 meas 0: 1.964e-02

Prep 0 meas 1: 6.610e-02
Prep 1 meas 0: 1.379e-02

Prep 0 meas 1: 6.097e-02
Prep 1 meas 0: 3.602e-03

Prep 0 meas 1: 7.018e-02
Prep 1 meas 0: 1.482e-03

Prep 0 meas 1: 3.115e-03
Prep 1 meas 0: 3.088e-02

Prep 0 meas 1: 1.939e-02
Prep 1 meas 0: 2.923e-03

Prep 0 meas 1: 1.063e-02
Prep 1 meas 0: 2.176e-02

Prep 0 meas 1: 2.897e-02
Prep 1 meas 0: 1.569e-02

Prep 0 meas 1: 2.730e-02
Prep 1 meas 0: 8.904e-04

Prep 0 meas 1: 5.401e-02
Prep 1 meas 0: 8.231e-02

Prep 0 meas 1: 7.207e-02
Prep 1 meas 0: 7.186e-02

Prep 0 meas 1: 2.724e-02
Prep 1 meas 0: 1.543e-03

Prep 0 meas 1: 1.384e-03
Prep 1 meas 0: 3.894e-02

Prep 0 meas 1: 8.157e-02
Prep 1 meas 0: 4.408e-03

Prep 0 meas 1: 9.377e-03
Prep 1 meas 0: 7.330e-02

Prep 0 meas 1: 1.973e-02
Prep 1 meas 0: 6.033e-05

Prep 0 meas 1: 1.075e-02
Prep 1 meas 0: 1.736e-02

Prep 0 meas 1: 2.048e-02
Prep 1 meas 0: 2.122e-03

Prep 0 meas 1: 1.235e-02
Prep 1 meas 0: 8.039e-02

Prep 0 meas 1: 1.381e-03
Prep 1 meas 0: 1.287e-02

Prep 0 meas 1: 3.622e-02
Prep 1 meas 0: 1.845e-02

Prep 0 meas 1: 2.039e-02
Prep 1 meas 0: 1.941e-02

Prep 0 meas 1: 1.496e-01
Prep 1 meas 0: 1.872e-02

Prep 0 meas 1: 2.454e-02
Prep 1 meas 0: 2.842e-02

Prep 0 meas 1: 2.178e-02
Prep 1 meas 0: 1.982e-02

Prep 0 meas 1: 5.102e-04
Prep 1 meas 0: 4.857e-03

Prep 0 meas 1: 6.100e-02
Prep 1 meas 0: 6.834e-02

Prep 0 meas 1: 8.164e-04
Prep 1 meas 0: 1.268e-02

Prep 0 meas 1: 3.096e-02
Prep 1 meas 0: 1.077e-02

Prep 0 meas 1: 9.210e-02
Prep 1 meas 0: 7.279e-03

Prep 0 meas 1: 1.877e-02
Prep 1 meas 0: 1.567e-03

Prep 0 meas 1: 1.298e-02
Prep 1 meas 0: 9.034e-02

Prep 0 meas 1: 5.746e-03
Prep 1 meas 0: 1.529e-03

Prep 0 meas 1: 4.415e-02
Prep 1 meas 0: 1.281e-02

Prep 0 meas 1: 3.364e-04
Prep 1 meas 0: 4.238e-03

Prep 0 meas 1: 1.677e-01
Prep 1 meas 0: 2.087e-02

Prep 0 meas 1: 2.944e-02
Prep 1 meas 0: 1.529e-02

Prep 0 meas 1: 2.757e-02
Prep 1 meas 0: 3.409e-02

Prep 0 meas 1: 2.196e-03
Prep 1 meas 0: 1.356e-02

Prep 0 meas 1: 4.640e-02
Prep 1 meas 0: 1.075e-01

Prep 0 meas 1: 1.141e-01
Prep 1 meas 0: 9.404e-02

Prep 0 meas 1: 4.400e-03
Prep 1 meas 0: 3.722e-02

Prep 0 meas 1: 1.308e-02
Prep 1 meas 0: 6.963e-04

Prep 0 meas 1: 7.421e-03
Prep 1 meas 0: 8.472e-02

Prep 0 meas 1: 1.042e-02
Prep 1 meas 0: 2.376e-02

Prep 0 meas 1: 4.594e-02
Prep 1 meas 0: 1.437e-04

Prep 0 meas 1: 6.883e-03
Prep 1 meas 0: 2.595e-02

Prep 0 meas 1: 2.412e-02
Prep 1 meas 0: 7.091e-02

Prep 0 meas 1: 5.190e-02
Prep 1 meas 0: 2.750e-02

Prep 0 meas 1: 8.021e-03
Prep 1 meas 0: 2.017e-01

Prep 0 meas 1: 1.773e-01
Prep 1 meas 0: 3.356e-02

Prep 0 meas 1: 9.497e-03
Prep 1 meas 0: 5.684e-04

Prep 0 meas 1: 2.711e-02
Prep 1 meas 0: 3.854e-02

Prep 0 meas 1: 5.583e-03
Prep 1 meas 0: 3.309e-02

Prep 0 meas 1: 4.179e-03
Prep 1 meas 0: 7.091e-03

Prep 0 meas 1: 1.783e-02
Prep 1 meas 0: 4.603e-02

Prep 0 meas 1: 7.258e-03
Prep 1 meas 0: 6.195e-02

Prep 0 meas 1: 2.481e-03
Prep 1 meas 0: 4.681e-02

Prep 0 meas 1: 7.330e-02
Prep 1 meas 0: 6.523e-03

Prep 0 meas 1: 2.089e-03
Prep 1 meas 0: 4.543e-03

Prep 0 meas 1: 3.731e-02
Prep 1 meas 0: 3.813e-04

Prep 0 meas 1: 7.714e-03
Prep 1 meas 0: 1.226e-03

Prep 0 meas 1: 9.860e-03
Prep 1 meas 0: 1.639e-02

Prep 0 meas 1: 3.073e-03
Prep 1 meas 0: 1.018e-01

Prep 0 meas 1: 2.898e-02
Prep 1 meas 0: 8.594e-02

Prep 0 meas 1: 1.144e-02
Prep 1 meas 0: 2.055e-02

Prep 0 meas 1: 3.395e-02
Prep 1 meas 0: 3.625e-03

Prep 0 meas 1: 5.848e-06
Prep 1 meas 0: 6.457e-03

Prep 0 meas 1: 1.315e-01
Prep 1 meas 0: 4.339e-02

Prep 0 meas 1: 7.851e-03
Prep 1 meas 0: 5.244e-02

Prep 0 meas 1: 3.493e-03
Prep 1 meas 0: 5.407e-02

Prep 0 meas 1: 8.774e-03
Prep 1 meas 0: 9.827e-02

Prep 0 meas 1: 2.518e-04
Prep 1 meas 0: 8.628e-02

Prep 0 meas 1: 1.262e-03
Prep 1 meas 0: 1.506e-03

Prep 0 meas 1: 1.191e-02
Prep 1 meas 0: 7.654e-02

Found 3 seeds with readout errors in range [5.0e-04, 4.0e-03]:

======================================================================
Seed   5: prep_0_meas_1 = 0.003745, prep_1_meas_0 = 0.003655
Seed  17: prep_0_meas_1 = 0.003145, prep_1_meas_0 = 0.000714
Seed  98: prep_0_meas_1 = 0.001262, prep_1_meas_0 = 0.001506

======================================================================
Summary Statistics:
======================================================================

	seed	prep_0_meas_1	prep_1_meas_0
count	3.000000	3.000000	3.000000
mean	40.000000	0.002717	0.001958
std	50.586559	0.001295	0.001522
min	5.000000	0.001262	0.000714
25%	11.000000	0.002204	0.001110
50%	17.000000	0.003145	0.001506
75%	57.500000	0.003445	0.002580
max	98.000000	0.003745	0.003655

GT = random_instrument(corr_strength=0.1, seed=98, fidelity=0.97)
GT.reveal()

GT.show_readout_errors()
GT.show_backaction_errors()

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.99872032, 0.00132291]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00001741, 0.00018317]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00003101, 0.00078360]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00123126, 0.99771031]]

Prep 0 meas 1: 1.262e-03
Prep 1 meas 0: 1.506e-03

Prep 0 ends in 1: 1.249e-03
Prep 1 ends in 0: 2.107e-03





[np.float64(0.0012486678967334145), np.float64(0.0021065159600563614)]

Parametrizing the MCM matrices based on learnable values & Plotting Gauge Transformations

def reconstruct_instrument_from_invariants_mixed_det(
    trM0: float,
    detM0: float,
    trM1: float,
    detM1: float,
    S0: float,               # S0 = 2 * 1^T M0 v0 with v0 = (1,1)/2, S0 is the sum of all elements of M0
    gauge_p00: float = 1.0,  # a = M0[0,0]
    noise_tol: float = 1e-12
) -> List[Tuple[np.ndarray, np.ndarray]]:
    """
    Inputs: trM0, detM0, trM1, detM1, S0, and gauge M0[0,0]=a.
    the parametrization is:
        M0 = [[a, b],
              [c, d]]
        M1 = [[e, f],
              [g, h]]
    Solve M0 from:
        d = trM0 - a
        b + c = S0 - trM0
        b * c = a*d - detM0
    Generally we will have two set of solutions for (b,c).
    Then solve M1 from:
        e+g = 1 - a - c   (= u)
        f+h = 1 - b - d   (= v)
        e+h = trM1        (= T1)
        eh - fg = detM1
    The last equation gives h linearly: h = [detM1 + v(u - T1)]/(u - v) unless u=v. (this is unlikely for general MCMs)
    If u≈v, fall back to tr(M0M1) computed from detM1 by identity.
    """
    a  = float(gauge_p00)
    T0 = float(trM0)
    d  = T0 - a
    D0 = float(detM0)
    T1 = float(trM1)
    D1 = float(detM1)

    # M0 via quadratic in c
    S = float(S0) - T0           # S = b + c
    K = a*d - D0                 # K = b*c
    disc = S*S - 4.0*K
    if disc < -noise_tol:
        raise ValueError(f"Negative discriminant {disc:.3e}. Check inputs.")
    r = np.sqrt(max(disc, 0.0))
    c_roots = [(S + r)/2.0, (S - r)/2.0]

    sols: List[Tuple[np.ndarray, np.ndarray]] = []
    seen_c: List[float] = []

    for c in c_roots:
        if any(abs(c - cc) < noise_tol for cc in seen_c):
            continue
        seen_c.append(c)
        b = S - c

        u = 1.0 - a - c
        v = 1.0 - b - d

        if abs(u - v) > noise_tol:
            # Linear solve for h from det:
            h = (D1 + v*(u - T1)) / (u - v)
            e = T1 - h
            g = u - e
            f = v - h
        else:
            # Degenerate case: use tr(M0M1) from the identity and do a 4x4 linear solve
            X = D1 + D0 + (T0 - 1.0)*(T1 - 1.0)  # tr(M0M1)
            A = np.array([[1.0, 0.0, 1.0, 0.0],
                          [0.0, 1.0, 0.0, 1.0],
                          [1.0, 0.0, 0.0, 1.0],
                          [a,   c,   b,   d  ]], dtype=float)
            y = np.array([u, v, T1, X], dtype=float)
            try:
                e, f, g, h = np.linalg.solve(A, y)
            except np.linalg.LinAlgError:
                sol, *_ = np.linalg.lstsq(A, y, rcond=None)
                e, f, g, h = sol
                if np.linalg.norm(A @ sol - y, ord=np.inf) > 1e-9:
                    continue  # reject

        M0 = np.array([[a, b],
                       [c, d]], dtype=float)
        M1 = np.array([[e, f],
                       [g, h]], dtype=float)
        sols.append((M0, M1))

    return sols


def summarize_instrument(M0: np.ndarray, M1: np.ndarray) -> Dict[str, Any]:
    S = M0 + M1
    return dict(
        trM0=float(np.trace(M0)),
        detM0=float(np.linalg.det(M0)),
        trM1=float(np.trace(M1)),
        detM1=float(np.linalg.det(M1)),
        trM0M1=float(np.trace(M0 @ M1)),
        colsum0=float(S[:,0].sum()),
        colsum1=float(S[:,1].sum()),
        S0=float(M0.sum()),
    )

# Analytically compute the allowed ranges for $t$ starting from one of the solution pairs

def _abgd_from_M(M):
    a, b = M[0][0], M[0][1]
    c, d = M[1][0], M[1][1]
    alpha = 0.5*(a+b+c+d)
    beta  = 0.5*(a+c - b - d)
    gamma = 0.5*(a+b - c - d)
    delta = 0.5*(a - b - c + d)
    return alpha, beta, gamma, delta

def _intervals_quad(a, b, c, rel, tol=1e-24):
    """
    Solve { x : a x^2 + b x + c (rel) 0 }, rel in {'ge','le'}.
    Returns a list of (lo, hi) with lo/hi allowed to be ±inf.
    """
    inf = float('inf')
    # Linear or constant
    if abs(a) < tol:
        if abs(b) < tol:
            ok = (c >= -tol) if rel == 'ge' else (c <= tol)
            return [(-inf, inf)] if ok else []
        x0 = -c/b
        if b > 0:
            return [(x0, inf)] if rel == 'ge' else [(-inf, x0)]
        else:
            return [(-inf, x0)] if rel == 'ge' else [(x0, inf)]

    # Quadratic
    D = b*b - 4*a*c
    if D < -tol:
        # No real roots
        if a > 0:
            return [(-inf, inf)] if rel == 'ge' else []
        else:
            return [] if rel == 'ge' else [(-inf, inf)]
    if D < 0:  # treat tiny negatives as zero
        D = 0.0
    sqrtD = np.sqrt(D)
    r1 = (-b - sqrtD)/(2*a)
    r2 = (-b + sqrtD)/(2*a)
    if r1 > r2:
        r1, r2 = r2, r1

    if a > 0:
        if rel == 'ge':
            return [(-inf, r1), (r2, inf)] if D > tol else [(-inf, r1), (r2, inf)]
        else:
            return [(r1, r2)] if D > tol else [(r1, r2)]
    else:  # a < 0
        if rel == 'ge':
            return [(r1, r2)] if D > tol else [(r1, r2)]
        else:
            return [(-inf, r1), (r2, inf)] if D > tol else [(-inf, r1), (r2, inf)]

def _intersect_interval_lists(A, B, tol=1e-24):
    C = []
    for lo1, hi1 in A:
        for lo2, hi2 in B:
            lo = max(lo1, lo2)
            hi = min(hi1, hi2)
            if lo <= hi + tol:
                C.append((lo, hi))
    if not C:
        return []
    C.sort(key=lambda x: x[0])
    merged = []
    cur_lo, cur_hi = C[0]
    for lo, hi in C[1:]:
        if lo <= cur_hi + tol:
            cur_hi = max(cur_hi, hi)
        else:
            merged.append((cur_lo, cur_hi))
            cur_lo, cur_hi = lo, hi
    merged.append((cur_lo, cur_hi))
    return merged

def _intersect_many(list_of_interval_lists, tol=1e-24):
    if not list_of_interval_lists:
        return []
    out = list_of_interval_lists[0]
    for L in list_of_interval_lists[1:]:
        out = _intersect_interval_lists(out, L, tol=tol)
        if not out:
            break
    return out

def _D_intervals_for_M(M, D_positive=True, tol=1e-24, margin_tol=0.0):
    """
    Return feasible D-intervals for a single M.
    For D>0: enforce -margin_tol <= a',b',c',d' <= 1+margin_tol.
    For D<0: same but inequality directions flip (because 2D < 0).
    """
    alpha, beta, gamma, delta = _abgd_from_M(M)

    # Adjust linear coefficients for margin
    # Lower bound: M' >= -eps => Q/(2D) >= -eps => Q + 2*eps*D >= 0 (for D>0)
    # Upper bound: M' <= 1+eps => Q/(2D) <= 1+eps => Q - 2*(1+eps)*D <= 0 => Q - 2D - 2*eps*D <= 0 (for D>0)
    shift = 2.0 * margin_tol

    # Build the eight quadratic constraints (lower/upper for a',b',c',d').
    # For D>0:
    #  a':  β D^2 + (α+δ+shift)D + γ ≥ 0;   β D^2 + (α+δ-2-shift)D + γ ≤ 0
    #  b': -β D^2 + (α-δ+shift)D + γ ≥ 0;  -β D^2 + (α-δ-2-shift)D + γ ≤ 0
    #  c':  β D^2 + (α-δ+shift)D - γ ≥ 0;   β D^2 + (α-δ-2-shift)D - γ ≤ 0
    #  d': -β D^2 + (α+δ+shift)D - γ ≥ 0;  -β D^2 + (α+δ-2-shift)D - γ ≤ 0
    if D_positive:
        polys = [
            ( beta,  alpha+delta+shift,     gamma, 'ge'),
            ( beta,  alpha+delta-2-shift,   gamma, 'le'),
            (-beta,  alpha-delta+shift,     gamma, 'ge'),
            (-beta,  alpha-delta-2-shift,   gamma, 'le'),
            ( beta,  alpha-delta+shift,    -gamma, 'ge'),
            ( beta,  alpha-delta-2-shift,  -gamma, 'le'),
            (-beta,  alpha+delta+shift,    -gamma, 'ge'),
            (-beta,  alpha+delta-2-shift,  -gamma, 'le'),
        ]
        domain = (0.0, float('inf'))
    else:
        # For D<0, reverse directions
        polys = [
            ( beta,  alpha+delta+shift,     gamma, 'le'),
            ( beta,  alpha+delta-2-shift,   gamma, 'ge'),
            (-beta,  alpha-delta+shift,     gamma, 'le'),
            (-beta,  alpha-delta-2-shift,   gamma, 'ge'),
            ( beta,  alpha-delta+shift,    -gamma, 'le'),
            ( beta,  alpha-delta-2-shift,  -gamma, 'ge'),
            (-beta,  alpha+delta+shift,    -gamma, 'le'),
            (-beta,  alpha+delta-2-shift,  -gamma, 'ge'),
        ]
        domain = (-float('inf'), 0.0)

    all_sets = []
    for a,b,c,rel in polys:
        S = _intervals_quad(a,b,c, rel, tol=tol)
        if not S:
            return []
        all_sets.append(S)

    inter = _intersect_many(all_sets, tol=tol)
    if not inter:
        return []
    return _intersect_interval_lists(inter, [domain], tol=tol)

def _merge_intervals(intervals, tol=1e-24):
    if not intervals:
        return []
    intervals = sorted(intervals, key=lambda x: x[0])
    merged = []
    lo, hi = intervals[0]
    for L, H in intervals[1:]:
        if L <= hi + tol:
            hi = max(hi, H)
        else:
            merged.append((lo, hi))
            lo, hi = L, H
    merged.append((lo, hi))
    return merged

def _map_D_to_t(D_intervals, tol=1e-24):
    """
    t = (1 - D)/2. Monotone decreasing map.
    Map each [D_lo, D_hi] to [t_lo, t_hi] with t_lo=(1-D_hi)/2, t_hi=(1-D_lo)/2.
    Excludes t=1/2 automatically because D-intervals never include D=0 (domain split).
    """
    out = []
    for D_lo, D_hi in D_intervals:
        t_lo = (1 - D_hi)/2
        t_hi = (1 - D_lo)/2
        out.append((t_lo, t_hi))
    return _merge_intervals(out, tol=tol)

def allowed_t_regions_for_M(M, tol=1e-24, margin_tol=0.0):
    """
    Input: M (2x2 numpy array, real).
    Output: list of (t_min, t_max) intervals such that R(t)^{-1} @ M @ R(t) has all entries in [-margin_tol, 1+margin_tol].
            Endpoints may be ±np.inf. t=1/2 is excluded by construction.
    """
    D_pos = _D_intervals_for_M(M, D_positive=True,  tol=tol, margin_tol=margin_tol)
    D_neg = _D_intervals_for_M(M, D_positive=False, tol=tol, margin_tol=margin_tol)
    t_sets = []
    if D_pos:
        t_sets += _map_D_to_t(D_pos, tol=tol)
    if D_neg:
        t_sets += _map_D_to_t(D_neg, tol=tol)
    return _merge_intervals(t_sets, tol=tol)

# Optional helper to intersect across multiple matrices at once:
def allowed_t_regions_for_list(M_list, tol=1e-24, margin_tol=0.0):
    """
    Intersect allowed regions across several matrices.
    """
    regions = None
    for M in M_list:
        r = allowed_t_regions_for_M(M, tol=tol, margin_tol=margin_tol)
        if regions is None:
            regions = r
        else:
            regions = _intersect_interval_lists(regions, r, tol=tol)
        if not regions:
            return []
    return _merge_intervals(regions, tol=tol)

def _R(t):
    """R(t) = [[1-t, t],[t, 1-t]], valid for t != 1/2."""
    return np.array([[1.0 - t, t],
                    [t, 1.0 - t]], dtype=float)

def _R_inv(t):
    """Closed-form inverse of R(t)."""
    D = 1.0 - 2.0*t
    if abs(D) < 1e-16:
        raise ValueError("t is too close to 1/2; R(t) is nearly singular.")
    return (1.0/D) * np.array([[1.0 - t, -t],
                            [-t, 1.0 - t]], dtype=float)

def _gauge_transform_pair(M0, M1, t):
    """Return (M0', M1') = (R^-1 M0 R, R^-1 M1 R)."""
    Rin = _R_inv(t); R = _R(t)
    return Rin @ M0 @ R, Rin @ M1 @ R

def rebase_and_anchor_instrument(M0, M1, t_regions, p00_min=0.5, tol=1e-24):
    """
    For each allowed t-interval [t_lo, t_hi], anchor at t_lo, transform (M0,M1),
    and return the anchored pair with its adjusted local-gauge interval [0, s_max],
    where s_max = (t_hi - t_lo)/(1 - 2*t_lo).

    Parameters
    ----------
    M0, M1 : (2,2) np.ndarray
        Trial instrument matrices.
    t_regions : list of (t_lo, t_hi)
        Allowed t-intervals from allowed_t_regions_for_list([M0, M1]).
    p00_min : float
        Required lower bound for anchored [M0']_{00}.
    tol : float
        Numerical tolerance for boundary checks.

    Returns
    -------
    results : list of dict
        Each dict contains:
            - 'anchor_t' : float
            - 'anchored_M0' : (2,2) np.ndarray
            - 'anchored_M1' : (2,2) np.ndarray
            - 'adjusted_region' : (0.0, s_max)  # can have s_max < 0 when anchor_t > 1/2
            - 'ok_p00' : bool  # whether [anchored_M0]_{00} >= p00_min (within tol)
            - 'p00' : float
    """
    results = []
    for (t_lo, t_hi) in t_regions:
        # Anchor at the left endpoint
        t1 = float(t_lo)
        M0a, M1a = _gauge_transform_pair(M0, M1, t1)

        # Check p_0^(0,0) >= p00_min
        p00 = float(M0a[0, 0])
        ok = (p00 >= p00_min - tol)

        # Compute s_max via group law: t3 = t1 + s - 2*t1*s
        denom = 1.0 - 2.0*t1
        if abs(denom) < 1e-16:
            # This should not occur if t_regions came from the feasibility solver
            s_max = np.sign(t_hi - t_lo) * np.inf
        else:
            s_max = (float(t_hi) - t1) / denom

        # Report adjusted local-gauge interval as [0, s_max]
        results.append({
            'anchor_t': t1,
            'anchored_M0': M0a,
            'anchored_M1': M1a,
            'adjusted_region': (0.0, s_max),
            'ok_p00': ok,
            'p00': p00
        })
    return results

def plot_gauge_transformation_effects(
    MCM_to_transform, 
    t_width_factor: float = 1.20, 
    verbose: bool = False, 
    MCM_reference: List[Instrument2x2] = None,
    resolution: float = 1e5,
    p00_min: float = 0.5,
    margin_tol: float = 0.0
):
    """
    Analyzes and plots the effect of a gauge transformation on an instrument.

    This function takes an instrument, applies a gauge transformation over a range
    of the gauge parameter 't', and plots how each of the 8 matrix entries evolves.
    It uses analytical methods to determine valid gauge parameter regions.

    Note: t=0.5 is excluded from the analysis as the gauge transformation matrix
    is non-invertible at that point.

    Args:
        MCM_to_transform: An Instrument2x2 object to be transformed.
        t_width_factor: Factor to scale the plotting range around valid regions. Default 1.20.
                        1.00 means plot exactly the valid regions, >1.00 adds padding.
        verbose: If True, print detailed information about valid intervals. Default False.
        MCM_reference: List of reference Instrument2x2 objects to compare RMSE against. Default [].
        resolution: Number of points to sample per unit t-range.
        p00_min: Minimum value for M0[0,0] when identifying focus regions. Default 0.5.
        margin_tol: Margin tolerance for allowed regions. Default 0.0.

    Returns:
        Tuple containing:
            - Instrument2x2 object constructed from center points of valid entry ranges, or None if no valid ranges.
            - List of valid t-regions for the center_instrument (relative to itself).
    """
    if MCM_reference is None:
        MCM_reference = []

    # Helper for formatting value column
    def fmt_val_err(min_v, max_v):
        c = (min_v + max_v) / 2
        h = (max_v - min_v) / 2

        # Determine exponent from the larger of abs(c) or abs(h) to avoid tiny numbers if c is near zero
        ref = abs(c) if abs(c) > 0 else abs(h)
        if ref == 0:
            return "(0.000 ± 0.000)e+0"

        exponent = int(np.floor(np.log10(ref)))
        scale = 10.0 ** (-exponent)

        c_s = c * scale
        h_s = h * scale

        return f"({c_s:.3f} ± {h_s:.3f})e{exponent:+d}"

    if verbose:
        print("Original Instrument to be transformed:")
        MCM_to_transform.reveal()

    # Use analytical method to find valid t-regions
    valid_t_regions = allowed_t_regions_for_list(
        [MCM_to_transform.M0, MCM_to_transform.M1], 
        tol=1e-24,
        margin_tol=margin_tol
    )

    if not valid_t_regions:
        raise ValueError(f"No valid gauge parameter regions found for this instrument (margin_tol={margin_tol}).")

    if verbose:
        print(f"\nAnalytically determined valid t-regions (total: {len(valid_t_regions)}) with margin {margin_tol}:")
        for i, (t_lo, t_hi) in enumerate(valid_t_regions, 1):
            print(f"  Region {i}: t ∈ [{t_lo:.6f}, {t_hi:.6f}] (width: {t_hi - t_lo:.6f})")

    # Determine plotting range based on valid regions and width factor
    all_t_mins = [r[0] for r in valid_t_regions if not np.isinf(r[0])]
    all_t_maxs = [r[1] for r in valid_t_regions if not np.isinf(r[1])]

    if all_t_mins and all_t_maxs:
        t_plot_min = min(all_t_mins)
        t_plot_max = max(all_t_maxs)
        t_center = (t_plot_min + t_plot_max) / 2
        t_half_span = (t_plot_max - t_plot_min) / 2

        # Apply width factor
        t_plot_min = t_center - t_half_span * t_width_factor
        t_plot_max = t_center + t_half_span * t_width_factor
    else:
        # Fallback if regions are unbounded
        t_plot_min = -0.5
        t_plot_max = 1.5

    # Ensure we don't include t=0.5 in our sampling
    if abs(t_plot_min - 0.5) < 1e-6:
        t_plot_min = 0.5 - 1e-6
    if abs(t_plot_max - 0.5) < 1e-6:
        t_plot_max = 0.5 + 1e-6

    # Generate t values for plotting, excluding t=0.5
    n_points = int(resolution * (t_plot_max - t_plot_min))
    if t_plot_min < 0.5 < t_plot_max:
        t_values_left = np.linspace(t_plot_min, 0.5 - 1e-6, n_points // 2)
        t_values_right = np.linspace(0.5 + 1e-6, t_plot_max, n_points // 2)
        t_values = np.concatenate([t_values_left, t_values_right])
    elif t_plot_max < 0.5:
        t_values = np.linspace(t_plot_min, t_plot_max, n_points)
    else:
        t_values = np.linspace(t_plot_min, t_plot_max, n_points)

    # Store the 8 entries of the transformed instrument for each value of t
    transformed_entries = []

    for t_val in t_values:
        M0_prime, M1_prime = gauge_transform_instrument_numerically(
            MCM_to_transform.M0, MCM_to_transform.M1, t_val
        )
        entries = np.concatenate((M0_prime.flatten(), M1_prime.flatten()))
        transformed_entries.append(entries)

    transformed_entries = np.array(transformed_entries)

    # --- Generate dedicated statistics samples from the first valid region ---
    # This ensures table values are independent of plotting width factor.
    stats_entries = None
    if len(valid_t_regions) > 0:
        t_lo_stat, t_hi_stat = valid_t_regions[0]

        # Handle potential infinite bounds for stats (clip to reasonable range if needed)
        t_start_stat = t_lo_stat if not np.isinf(t_lo_stat) else -5.0
        t_end_stat = t_hi_stat if not np.isinf(t_hi_stat) else 5.0

        width_stat = t_end_stat - t_start_stat
        # Ensure sufficient points for statistics
        n_stats = max(200, int(resolution * width_stat))
        t_stats = np.linspace(t_start_stat, t_end_stat, n_stats)

        stats_entries_list = []
        for t_val in t_stats:
            M0_p, M1_p = gauge_transform_instrument_numerically(
                MCM_to_transform.M0, MCM_to_transform.M1, t_val
            )
            stats_entries_list.append(np.concatenate((M0_p.flatten(), M1_p.flatten())))
        stats_entries = np.array(stats_entries_list)
    # -----------------------------------------------------------------------

    # Use rebase_and_anchor_instrument to find focus regions with ok_p00
    rebased_results = rebase_and_anchor_instrument(
        MCM_to_transform.M0, 
        MCM_to_transform.M1, 
        valid_t_regions, 
        p00_min=p00_min
    )

    # Find first region with ok_p00 for focus plot
    focus_region_info = None
    for res in rebased_results:
        if res['ok_p00']:
            focus_region_info = res
            break

    # Determine number of subplots
    has_focus_plot = focus_region_info is not None or len(MCM_reference) > 0

    if has_focus_plot:
        fig = plt.figure(figsize=(14, 28))
        gs = fig.add_gridspec(4, 1, height_ratios=[1.2, 1.2, 1.2, 2.4], hspace=0.15)

        ax1 = fig.add_subplot(gs[0])
        ax2 = fig.add_subplot(gs[1], sharex=ax1)
        axes = [ax1, ax2]

        if len(MCM_reference) > 0:
            ax3 = fig.add_subplot(gs[2], sharex=ax1)
            axes.append(ax3)

        ax4 = fig.add_subplot(gs[3])
        axes.append(ax4)
    else:
        fig, axes = plt.subplots(2, 1, figsize=(14, 14), sharex=True)
        if isinstance(axes, plt.Axes):
            axes = [axes]

    ax1 = axes[0]
    ax2 = axes[1]

    # ===== First subplot: Individual matrix entries =====
    labels = [
        r"$M^0$[0,0] = $p_0^{(0,0)}$", 
        r"$M^0$[0,1] = $p_1^{(0,0)}$", 
        r"$M^0$[1,0] = $p_0^{(0,1)}$",
        r"$M^0$[1,1] = $p_1^{(0,1)}$",
        r"$M^1$[0,0] = $p_0^{(1,0)}$",
        r"$M^1$[0,1] = $p_1^{(1,0)}$",
        r"$M^1$[1,0] = $p_0^{(1,1)}$",
        r"$M^1$[1,1] = $p_1^{(1,1)}$"
    ]

    for i in range(8):
        ax1.plot(t_values, transformed_entries[:, i], linewidth=1.0, label=labels[i])

    ax1.axvline(x=0.5, color='red', linestyle=':', linewidth=1.0, label='Singularity (t=0.5)')
    ax1.axhline(y=0, color='k', linestyle='--', linewidth=1.0, label='Prob Boundary (0,1)')
    ax1.axhline(y=1, color='k', linestyle='--', linewidth=1.0)

    if margin_tol > 0:
        ax1.axhline(y=-margin_tol, color='gray', linestyle=':', linewidth=1.0, label=f'Margin (±{margin_tol})')
        ax1.axhline(y=1+margin_tol, color='gray', linestyle=':', linewidth=1.0)

    # Highlight valid regions using analytical results
    for i, (t_lo, t_hi) in enumerate(valid_t_regions):
        # Clip to plotting range
        t_lo_plot = max(t_lo, t_plot_min) if not np.isinf(t_lo) else t_plot_min
        t_hi_plot = min(t_hi, t_plot_max) if not np.isinf(t_hi) else t_plot_max

        if t_lo_plot < t_hi_plot:
            label = 'Valid Gauge Region (t)' if i == 0 else ""
            ax1.axvspan(t_lo_plot, t_hi_plot, color='green', alpha=0.2, label=label)

    ax1.set_ylabel("Value of Instrument Matrix Entry")
    ax1.set_title(f"Evolution of Instrument Entries (Valid Interval: [{-margin_tol}, {1+margin_tol}])")
    ax1.legend(loc='center left', bbox_to_anchor=(1, 0.5))
    ax1.grid(True, linestyle=':', alpha=0.6)
    ax1.set_ylim(max(np.min(transformed_entries), -0.1 - margin_tol), min(np.max(transformed_entries), 1.1 + margin_tol))

    # ===== Compute valid ranges for matrix entries =====
    entry_labels = [
        "M^0[0,0] = p_0^(0,0)",
        "M^0[0,1] = p_1^(0,0)",
        "M^0[1,0] = p_0^(0,1)",
        "M^0[1,1] = p_1^(0,1)",
        "M^1[0,0] = p_0^(1,0)",
        "M^1[0,1] = p_1^(1,0)",
        "M^1[1,0] = p_0^(1,1)",
        "M^1[1,1] = p_1^(1,1)"
    ]

    entry_ranges_data = []
    center_values = []  # Store center values for constructing return instrument

    for entry_idx, entry_label in enumerate(entry_labels):
        if stats_entries is not None:
            block_values = stats_entries[:, entry_idx]
            min_val = np.min(block_values)
            max_val = np.max(block_values)
            width = max_val - min_val
            center = (min_val + max_val) / 2

            center_values.append(center)

            row_data = {
                'Entry': entry_label,
                'Min': f"{min_val:.8f}",
                'Max': f"{max_val:.8f}",
                'Width': f"{width:.5e}",
                'Value': fmt_val_err(min_val, max_val)
            }

            # Add comparison columns for each reference instrument
            for ref_idx, ref_inst in enumerate(MCM_reference):
                ref_entries = np.concatenate((ref_inst.M0.flatten(), ref_inst.M1.flatten()))
                ref_val = ref_entries[entry_idx]

                # Calculate absolute discrepancy from center
                abs_disc = ref_val - center

                # Check if reference value is within range
                if min_val <= ref_val <= max_val:
                    row_data[f'Ref{ref_idx+1}'] = f"IN ({abs_disc:+.5e})"
                else:
                    direction = "above" if ref_val > max_val else "below"
                    row_data[f'Ref{ref_idx+1}'] = f"OUT {direction} ({abs_disc:+.5e})"

            entry_ranges_data.append(row_data)

    # Construct Instrument2x2 from center values
    center_instrument = None
    if len(center_values) == 8:
        M0_center = np.array([[center_values[0], center_values[1]],
                              [center_values[2], center_values[3]]], dtype=float)
        M1_center = np.array([[center_values[4], center_values[5]],
                              [center_values[6], center_values[7]]], dtype=float)
        center_instrument = Instrument2x2(M0=M0_center, M1=M1_center)

    # ===== Second subplot: Derived quantities =====
    prep0_meas1 = transformed_entries[:, 4] + transformed_entries[:, 6]
    prep1_meas0 = transformed_entries[:, 1] + transformed_entries[:, 3]
    prep0_excite = transformed_entries[:, 2] + transformed_entries[:, 6]
    prep1_decay = transformed_entries[:, 1] + transformed_entries[:, 5]

    derived_quantities = {
        "prep 0 meas 1": prep0_meas1,
        "prep 1 meas 0": prep1_meas0,
        "prep 0 excite to 1": prep0_excite,
        "prep 1 decay to 0": prep1_decay
    }

    # Calculate derived quantities for stats
    stats_derived = {}
    if stats_entries is not None:
        s_prep0_meas1 = stats_entries[:, 4] + stats_entries[:, 6]
        s_prep1_meas0 = stats_entries[:, 1] + stats_entries[:, 3]
        s_prep0_excite = stats_entries[:, 2] + stats_entries[:, 6]
        s_prep1_decay = stats_entries[:, 1] + stats_entries[:, 5]

        stats_derived = {
            "prep 0 meas 1": s_prep0_meas1,
            "prep 1 meas 0": s_prep1_meas0,
            "prep 0 excite to 1": s_prep0_excite,
            "prep 1 decay to 0": s_prep1_decay
        }
    # print(derived_quantities)
    # print(stats_derived)

    quantity_valid_ranges = {}
    quantity_ranges_data = []
    for quantity_name, quantity_values in derived_quantities.items():
        if quantity_name in stats_derived:
            block_values = stats_derived[quantity_name]

            min_val = np.min(block_values)
            max_val = np.max(block_values)
            width = max_val - min_val
            center = (min_val + max_val) / 2
            quantity_valid_ranges[quantity_name] = (min_val, max_val)

            row_data = {
                'Quantity': quantity_name,
                'Min': f"{min_val:.8f}",
                'Max': f"{max_val:.8f}",
                'Width': f"{width:.5e}",
                'Value': fmt_val_err(min_val, max_val)
            }

            # Add comparison columns for each reference instrument
            for ref_idx, ref_inst in enumerate(MCM_reference):
                # Compute reference quantity value
                if quantity_name == "prep 0 meas 1":
                    ref_val = ref_inst.M1[0, 0] + ref_inst.M1[1, 0]
                elif quantity_name == "prep 1 meas 0":
                    ref_val = ref_inst.M0[0, 1] + ref_inst.M0[1, 1]
                elif quantity_name == "prep 0 excite to 1":
                    ref_val = ref_inst.M0[1, 0] + ref_inst.M1[1, 0]
                elif quantity_name == "prep 1 decay to 0":
                    ref_val = ref_inst.M0[0, 1] + ref_inst.M1[0, 1]
                else:
                    ref_val = 0.0

                # Calculate absolute discrepancy from center
                abs_disc = ref_val - center

                # Check if reference value is within range
                if min_val <= ref_val <= max_val:
                    row_data[f'Ref{ref_idx+1}'] = f"IN ({abs_disc:+.5e})"
                else:
                    direction = "above" if ref_val > max_val else "below"
                    row_data[f'Ref{ref_idx+1}'] = f"OUT {direction} ({abs_disc:+.5e})"

            quantity_ranges_data.append(row_data)

    # Plot the derived quantities
    labels_with_ranges = [
        (r"prep 0 meas 1: $p_0^{(1,0)} + p_0^{(1,1)}$", "prep 0 meas 1"),
        (r"prep 1 meas 0: $p_1^{(0,0)} + p_1^{(0,1)}$", "prep 1 meas 0"),
        (r"prep 0 excite to 1: $p_0^{(0,1)} + p_0^{(1,1)}$", "prep 0 excite to 1"),
        (r"prep 1 decay to 0: $p_1^{(0,0)} + p_1^{(1,0)}$", "prep 1 decay to 0")
    ]

    quantity_list = list(derived_quantities.items())
    for idx, ((base_label, quantity_key), (quantity_name, quantity_values)) in enumerate(zip(labels_with_ranges, quantity_list)):
        if quantity_key in quantity_valid_ranges:
            min_val, max_val = quantity_valid_ranges[quantity_key]
            label = f"{base_label}\n∈ [{min_val:.6f}, {max_val:.6f}]"
        else:
            label = base_label
        ax2.plot(t_values, quantity_values, label=label, linewidth=1)

    ax2.axvline(x=0.5, color='red', linestyle=':', linewidth=1.0, label='Singularity (t=0.5)')
    ax2.axhline(y=0, color='k', linestyle='--', linewidth=1.0, label='Prob Boundary (0,1)')
    ax2.axhline(y=1, color='k', linestyle='--', linewidth=1.0)

    if margin_tol > 0:
        ax2.axhline(y=-margin_tol, color='gray', linestyle=':', linewidth=1.0, label=f'Margin (±{margin_tol})')
        ax2.axhline(y=1+margin_tol, color='gray', linestyle=':', linewidth=1.0)

    for i, (t_lo, t_hi) in enumerate(valid_t_regions):
        t_lo_plot = max(t_lo, t_plot_min) if not np.isinf(t_lo) else t_plot_min
        t_hi_plot = min(t_hi, t_plot_max) if not np.isinf(t_hi) else t_plot_max

        if t_lo_plot < t_hi_plot:
            label = 'Valid Gauge Region (t)' if i == 0 else ""
            ax2.axvspan(t_lo_plot, t_hi_plot, color='green', alpha=0.2, label=label)

    ax2.set_xlabel(r"Gauge Parameter $(t)$")
    ax2.set_ylabel("Derived Quantity Value")
    ax2.set_title(f"Derived Quantities (Valid Interval: [{-margin_tol}, {1+margin_tol}])")

    legend = ax2.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize=9, 
                        labelspacing=1.2, handlelength=2)
    ax2.grid(True, linestyle=':', alpha=0.6)
    ax2.set_ylim(-0.1 - margin_tol, 1.1 + margin_tol)

    # ===== Display DataFrames for valid intervals =====
    print("="*80)
    print("Valid Intervals for Derived Quantities (Readout & Back-action Errors)")
    print("="*80)
    if quantity_ranges_data:
        df_quantities = pd.DataFrame(quantity_ranges_data)
        display(df_quantities)
    else:
        print("No valid intervals found for derived quantities.")

    print("\n" + "="*80)
    print("Valid Intervals for Matrix Entries")
    print("="*80)
    if entry_ranges_data:
        df_entries = pd.DataFrame(entry_ranges_data)
        display(df_entries)
    else:
        print("No valid intervals found for matrix entries.")

    if len(MCM_reference) > 0:
        print("\nNote: Reference comparison format:")
        print("  'IN (±X.XXXe±YY)' - value is within range, absolute discrepancy from center")
        print("  'OUT above/below (±X.XXXe±YY)' - value is outside range, absolute discrepancy from center")
    print("="*80 + "\n")

    # ===== Third subplot: RMSE to reference instruments =====
    best_match_instruments = []

    if len(MCM_reference) > 0:
        ax3 = axes[2]

        for ref_idx, ref_inst in enumerate(MCM_reference):
            ref_entries = np.concatenate((ref_inst.M0.flatten(), ref_inst.M1.flatten()))
            rmse_values = np.sqrt(np.mean((transformed_entries - ref_entries)**2, axis=1))

            min_rmse_idx = np.argmin(rmse_values)
            min_rmse = rmse_values[min_rmse_idx]
            t_min_rmse = t_values[min_rmse_idx]

            label = f"Ref {ref_idx+1}: min RMSE={min_rmse:.6e} at t={t_min_rmse:.4f}"
            ax3.plot(t_values, rmse_values, label=label, linewidth=1.5)

            ax3.plot(t_min_rmse, min_rmse, 'o', markersize=8)

            M0_best, M1_best = gauge_transform_instrument_numerically(
                MCM_to_transform.M0, MCM_to_transform.M1, t_min_rmse
            )
            best_match_instruments.append(Instrument2x2(M0=M0_best, M1=M1_best))

        ax3.axvline(x=0.5, color='red', linestyle=':', linewidth=1.0, label='Singularity (t=0.5)')

        for i, (t_lo, t_hi) in enumerate(valid_t_regions):
            t_lo_plot = max(t_lo, t_plot_min) if not np.isinf(t_lo) else t_plot_min
            t_hi_plot = min(t_hi, t_plot_max) if not np.isinf(t_hi) else t_plot_max

            if t_lo_plot < t_hi_plot:
                label = 'Valid Gauge Region (t)' if i == 0 else ""
                ax3.axvspan(t_lo_plot, t_hi_plot, color='green', alpha=0.2, label=label)

        ax3.set_ylabel("RMSE to Reference")
        ax3.set_title("RMSE Between Gauge-Transformed and Reference Instruments")
        ax3.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize=9)
        ax3.grid(True, linestyle=':', alpha=0.6)
        ax3.set_yscale('log')

    # ===== Fourth subplot: Focused view =====
    if has_focus_plot and focus_region_info is not None:
        ax4 = axes[-1]

        # Use the adjusted region from rebase_and_anchor_instrument
        anchor_t = focus_region_info['anchor_t']
        s_min, s_max = focus_region_info['adjusted_region']

        # Convert back to absolute t values
        # t3 = t1 + s - 2*t1*s => given t1=anchor_t, s in [s_min, s_max]
        t_focus_exact_min = anchor_t + s_min - 2*anchor_t*s_min
        t_focus_exact_max = anchor_t + s_max - 2*anchor_t*s_max

        if t_focus_exact_min > t_focus_exact_max:
            t_focus_exact_min, t_focus_exact_max = t_focus_exact_max, t_focus_exact_min

        # Apply t_width_factor to the focus region
        t_focus_center = (t_focus_exact_min + t_focus_exact_max) / 2
        t_focus_half_span = (t_focus_exact_max - t_focus_exact_min) / 2

        t_focus_min = t_focus_center - t_focus_half_span * t_width_factor
        t_focus_max = t_focus_center + t_focus_half_span * t_width_factor

        if verbose:
            print(f"\nFocus region: anchor_t={anchor_t:.6f}, "
                  f"local s∈[{s_min:.6f}, {s_max:.6f}]")
            print(f"  Exact valid t∈[{t_focus_exact_min:.6f}, {t_focus_exact_max:.6f}]")
            print(f"  Plotted t∈[{t_focus_min:.6f}, {t_focus_max:.6f}] (with factor {t_width_factor:.2f})")

        focus_mask = (t_values >= t_focus_min) & (t_values <= t_focus_max)
        t_focus = t_values[focus_mask]
        entries_focus = transformed_entries[focus_mask]

        for i in range(8):
            ax4.plot(t_focus, entries_focus[:, i], linewidth=0.8, alpha=0.5, label=labels[i])

        for quantity_name, quantity_values in derived_quantities.items():
            ax4.plot(t_focus, quantity_values[focus_mask], linewidth=1.5, label=quantity_name)

        if len(MCM_reference) > 0:
            for ref_idx, ref_inst in enumerate(MCM_reference):
                ref_entries = np.concatenate((ref_inst.M0.flatten(), ref_inst.M1.flatten()))
                rmse_focus = np.sqrt(np.mean((entries_focus - ref_entries)**2, axis=1))
                rmse_normalized = rmse_focus / (rmse_focus.max() + 1e-18)
                ax4.plot(t_focus, rmse_normalized, linewidth=2, linestyle='--', 
                        label=f"Ref {ref_idx+1} RMSE (normalized)\n original max={rmse_focus.max():.2e}")

        if 0.5 >= t_focus_min and 0.5 <= t_focus_max:
            ax4.axvline(x=0.5, color='red', linestyle=':', linewidth=1.0, label='Singularity (t=0.5)')

        ax4.axhline(y=0, color='k', linestyle='--', linewidth=1.0, alpha=0.5)
        ax4.axhline(y=1, color='k', linestyle='--', linewidth=1.0, alpha=0.5)

        # Plot valid regions - only those overlapping with focus window
        for i, (t_lo, t_hi) in enumerate(valid_t_regions):
            if t_hi >= t_focus_min and t_lo <= t_focus_max:
                plot_start = max(t_lo, t_focus_min)
                plot_end = min(t_hi, t_focus_max)
                label = 'Valid Gauge Region (t)' if i == 0 else ""
                ax4.axvspan(float(plot_start), float(plot_end), color='green', alpha=0.2, label=label)

        ax4.set_xlabel(r"Gauge Parameter $(t)$ [Focused View]")
        ax4.set_ylabel("Quantity Values")
        ax4.set_title(f"Focused View: t ∈ [{t_focus_min:.4f}, {t_focus_max:.4f}]\n"
                     f"Valid region: [{t_focus_exact_min:.4f}, {t_focus_exact_max:.4f}] "
                     f"(width factor: {t_width_factor:.2f})")
        ax4.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize=8, ncol=2)
        ax4.grid(True, linestyle=':', alpha=0.6)
        ax4.set_ylim(-0.1, 1.1)
        ax4.set_xlim(t_focus_min, t_focus_max)
    else:
        axes[-1].set_xlabel(r"Gauge Parameter $(t)$")

    plt.tight_layout()
    plt.show()

    center_t_regions = []
    if center_instrument is not None:
        center_t_regions = allowed_t_regions_for_list(
            [center_instrument.M0, center_instrument.M1], 
            tol=1e-24,
            margin_tol=margin_tol
        )
    # filter the center_t_regions such that the regions that has the smallest sum of absolute values for 2 end points is chosen
    if center_t_regions:
        center_t_regions = sorted(center_t_regions, key=lambda x: abs(x[0]) + abs(x[1]))
        center_t_regions = center_t_regions[0]

    return center_instrument, center_t_regions

def refine_bounds_with_sp_error(
    reconstructed_MCM: Instrument2x2 | Any,
    t_bound: Tuple[float, float] | Any,
    prob_meas_0_val: float | Any,
    margin_tol: float = 1e-4,
    verbose: bool = False
) -> Tuple[Tuple[float, float], Tuple[float, float]]:
    """
    Refines the valid gauge parameter interval (t_bound) by incorporating the constraint
    that the state preparation error (epsilon) must be physical (within [0, 1]).

    Logic:
        P(outcome=1) = P(1|0_noisy)
                     = P(1|0_ideal) * (1 - epsilon) + P(1|1_ideal) * epsilon

        Where P(1|s_ideal) are functions of the gauge parameter t derived from reconstructed_MCM.
        We solve for epsilon(t) and find the subset of t_bound where 0 <= epsilon(t) <= 1.

    Args:
        reconstructed_MCM: The Instrument2x2 object at the center of the gauge (t=0).
        t_bound: The current valid interval (t_min, t_max) derived from MCM positivity.
        prob_meas_0_val: The empirical probability of measuring '0' (from prob dictionary).
        margin_tol: Tolerance for the physical bounds of epsilon (allows -tol to 1+tol).
        verbose: If True, prints detailed intermediate values.

    Returns:
        Tuple containing:
            - refined_t_bound: (new_t_min, new_t_max)
            - epsilon_bound: (eps_min, eps_max) over the refined t interval.
    """

    # 1. Extract baseline readout probabilities (at t=0) for outcome 1
    # M1 is the matrix for outcome 1.
    # Column 0 sum is P(1|0) at t=0.
    # Column 1 sum is P(1|1) at t=0.
    M1 = reconstructed_MCM.M1
    c10 = M1[0, 0] + M1[1, 0]  # P(1|0) baseline
    c11 = M1[0, 1] + M1[1, 1]  # P(1|1) baseline

    # We use P(outcome=1) for the calculation
    P_obs = 1.0 - prob_meas_0_val

    # Delta c = P(1|1) - P(1|0). Usually close to 1 for a good measurement.
    delta_c = c11 - c10

    if abs(delta_c) < 1e-9:
        if verbose:
            print("Warning: Instrument has negligible distinguishing power (delta_c ~ 0). Cannot bound epsilon.")
        return t_bound, (0.0, 0.0) # Return original bounds

    # 2. Define the mapping functions
    # Derived formula: epsilon(t) = [P_obs - c10 - t*delta_c] / [delta_c * (1 - 2t)]
    def get_epsilon(t):
        denom = delta_c * (1.0 - 2.0*t)
        if abs(denom) < 1e-15:
            return np.inf # Singularity at t=0.5
        num = P_obs - c10 - t * delta_c
        return num / denom

    # Inverse mapping: t(epsilon)
    def get_t(eps):
        denom = delta_c * (1.0 - 2.0*eps)
        if abs(denom) < 1e-15:
            return np.inf
        num = P_obs - c10 - eps * delta_c
        return num / denom

    # 3. Analyze bounds
    t_min_in, t_max_in = t_bound

    # Check for singularity crossing (unlikely for valid physical regions)
    if (t_min_in < 0.5 and t_max_in > 0.5):
        if verbose:
            print("Warning: Input t_bound crosses singularity t=0.5. Skipping refinement.")
        return t_bound, (get_epsilon(t_min_in), get_epsilon(t_max_in))

    # Calculate epsilon at the current t boundaries
    eps_at_min = get_epsilon(t_min_in)
    eps_at_max = get_epsilon(t_max_in)

    # Determine current epsilon range (order matters for min/max)
    curr_eps_min = min(eps_at_min, eps_at_max)
    curr_eps_max = max(eps_at_min, eps_at_max)

    # Define target physical range for epsilon
    target_eps_min = 0.0 - margin_tol
    target_eps_max = 1.0 + margin_tol

    if verbose:
        print(f"--- Refine Bounds Debug ---")
        print(f"P(outcome=1): {P_obs:.6f}")
        print(f"Baseline P(1|0): {c10:.6f}, P(1|1): {c11:.6f}")
        print(f"Input t_bound: [{t_min_in:.6f}, {t_max_in:.6f}]")
        print(f"Epsilon range on input t: [{curr_eps_min:.6f}, {curr_eps_max:.6f}]")
        print(f"Target epsilon range: [{target_eps_min:.6f}, {target_eps_max:.6f}]")

    # 4. Intersect ranges
    # We clip the current epsilon range to the target physical range
    refined_eps_min = max(curr_eps_min, target_eps_min)
    refined_eps_max = min(curr_eps_max, target_eps_max)

    # Check for empty intersection
    if refined_eps_min > refined_eps_max:
        print("Error: No valid solution found. State preparation error constraints are incompatible with MCM constraints.")
        return t_bound, (np.nan, np.nan)

    # 5. Map back to t
    # Since the mapping is monotonic (hyperbola branch), we map the refined epsilon endpoints back to t
    t_lim_1 = get_t(refined_eps_min)
    t_lim_2 = get_t(refined_eps_max)

    new_t_min = min(t_lim_1, t_lim_2)
    new_t_max = max(t_lim_1, t_lim_2)

    # Numerical safety clip to original bounds
    new_t_min = max(new_t_min, t_min_in)
    new_t_max = min(new_t_max, t_max_in)

    # Re-evaluate epsilon at exact new t bounds for consistency
    final_eps_min = min(get_epsilon(new_t_min), get_epsilon(new_t_max))
    final_eps_max = max(get_epsilon(new_t_min), get_epsilon(new_t_max))

    if verbose:
        print(f"Refined t_bound: [{new_t_min:.6f}, {new_t_max:.6f}]")
        print(f"Refined epsilon: [{final_eps_min:.6f}, {final_eps_max:.6f}]")
        print(f"---------------------------")

    return (new_t_min, new_t_max), (final_eps_min, final_eps_max)

GT = random_instrument(corr_strength=0.1, seed=59, fidelity=0.8)
GT.reveal()

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.63686009, 0.14318407]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.26631059, 0.15142864]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.02935295, 0.03553478]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.06747637, 0.66985250]]





0

ibm_pittsburgh_mcm[3].reveal()

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.99872032, 0.00132291]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00001741, 0.00018317]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00003101, 0.00078360]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00123126, 0.99771031]]





0

chosen_sol_index = 0

GT = ibm_pittsburgh_mcm[3]
GT.reveal()

prob_dict_GT = calculate_exact_all_string_probabilities_from_v0_and_instrument(
    inst = GT,
    v0=np.array([0.5, 0.5]),
    max_len=3
)

M0t, M1t = GT.M0, GT.M1
invrnts = summarize_instrument(M0t, M1t)
sols = reconstruct_instrument_from_invariants_mixed_det(
    trM0=invrnts["trM0"], detM0=invrnts["detM0"],
    trM1=invrnts["trM1"], detM1=invrnts["detM1"],
    S0=invrnts["S0"], gauge_p00=1.0
)

display('')

reconstruct_instrument, _ = plot_gauge_transformation_effects(
    Instrument2x2(M0=sols[chosen_sol_index][0], M1=sols[chosen_sol_index][1]), 
    t_width_factor=1.1, 
    verbose=False, 
    MCM_reference=[GT]
    )

if reconstruct_instrument is not None:
    print("Reconstructed Instrument from center points of valid entry ranges:")
    reconstruct_instrument.reveal()

print("\nGround Truth Instrument:")
GT.reveal()

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.99872032, 0.00132291]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00001741, 0.00018317]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00003101, 0.00078360]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00123126, 0.99771031]]



''


================================================================================
Valid Intervals for Derived Quantities (Readout & Back-action Errors)
================================================================================

	Quantity	Min	Max	Width	Value	Ref1
0	prep 0 meas 1	0.00003005	0.00127966	1.24961e-03	(6.549 ± 6.248)e-4	IN (+6.07420e-04)
1	prep 1 meas 0	0.00027386	0.00152347	1.24961e-03	(8.987 ± 6.248)e-4	IN (+6.07420e-04)
2	prep 0 excite to 1	0.00124865	0.00124973	1.07234e-06	(1.249 ± 0.001)e-3	IN (-5.21213e-07)
3	prep 1 decay to 0	0.00210546	0.00210653	1.07234e-06	(2.106 ± 0.001)e-3	IN (+5.21213e-07)

================================================================================
Valid Intervals for Matrix Entries
================================================================================

	Entry	Min	Max	Width	Value	Ref1
0	M^0[0,0] = p_0^(0,0)	0.99871989	0.99872034	4.49770e-07	(9.987 ± 0.000)e-1	IN (+2.01820e-07)
1	M^0[0,1] = p_1^(0,0)	0.00009060	0.00134032	1.24973e-03	(7.155 ± 6.249)e-4	IN (+6.07453e-04)
2	M^0[1,0] = p_0^(0,1)	0.00000000	0.00124973	1.24973e-03	(6.249 ± 6.249)e-4	IN (-6.07453e-04)
3	M^0[1,1] = p_1^(0,1)	0.00018315	0.00018360	4.49770e-07	(1.834 ± 0.002)e-4	IN (-2.01820e-07)
4	M^1[0,0] = p_0^(1,0)	0.00003005	0.00003106	1.01727e-06	(3.056 ± 0.051)e-5	IN (+4.58425e-07)
5	M^1[0,1] = p_1^(1,0)	0.00076621	0.00201486	1.24865e-03	(1.391 ± 0.624)e-3	IN (-6.06932e-04)
6	M^1[1,0] = p_0^(1,1)	-0.00000000	0.00124865	1.24865e-03	(6.243 ± 6.243)e-4	IN (+6.06932e-04)
7	M^1[1,1] = p_1^(1,1)	0.99771026	0.99771128	1.01727e-06	(9.977 ± 0.000)e-1	IN (-4.58425e-07)

Note: Reference comparison format:
  'IN (±X.XXXe±YY)' - value is within range, absolute discrepancy from center
  'OUT above/below (±X.XXXe±YY)' - value is outside range, absolute discrepancy from center
================================================================================



/var/folders/dk/_dd6ng8n3yv_crq6vrdt8qzc0000gn/T/ipykernel_62248/724936655.py:540: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()

Reconstructed Instrument from center points of valid entry ranges:

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.99872012, 0.00071546]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00062486, 0.00018337]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00003056, 0.00139053]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00062433, 0.99771077]]

Ground Truth Instrument:

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.99872032, 0.00132291]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00001741, 0.00018317]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00003101, 0.00078360]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00123126, 0.99771031]]





0

ibm_pittsburgh_mcm = [
    random_instrument(corr_strength=0.005, seed=3, fidelity=0.98),
    random_instrument(corr_strength=0.08, seed=37, fidelity=0.96), # this might be a little little bit off
    random_instrument(corr_strength=0.2, seed=1, fidelity=0.99), # this is good for demo of a ideally near-perfect instrument
    random_instrument(corr_strength=0.1, seed=98, fidelity=0.97), # this is a good demo for ibm_pittsburgh in the morning of Nov 6, 2025.
]

GT = random_instrument(corr_strength=0.1, seed=98, fidelity=0.97)
GT.reveal()

GT.show_readout_errors()
GT.show_backaction_errors()

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.99872032, 0.00132291]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00001741, 0.00018317]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00003101, 0.00078360]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00123126, 0.99771031]]

Prep 0 meas 1: 1.262e-03
Prep 1 meas 0: 1.506e-03

Prep 0 ends in 1: 1.249e-03
Prep 1 ends in 0: 2.107e-03





[np.float64(0.0012486678967334145), np.float64(0.0021065159600563614)]

chosen_sol_index = 0

reconstructed_MCM, _ = plot_gauge_transformation_effects(
    Instrument2x2(M0=sols[chosen_sol_index][0], M1=sols[chosen_sol_index][1]), 
    t_width_factor=1.1, 
    verbose=False, 
    MCM_reference=[GT],
    resolution=1e5,
    )

print("Original Instrument:")
GT.reveal()

print("\nReconstructed Instrument from center points of valid entry ranges:")
if reconstructed_MCM is not None:
    reconstructed_MCM.reveal()

================================================================================
Valid Intervals for Derived Quantities (Readout & Back-action Errors)
================================================================================

	Quantity	Min	Max	Width	Value	Ref1
0	prep 0 meas 1	0.00003005	0.00127966	1.24961e-03	(6.549 ± 6.248)e-4	IN (+6.07420e-04)
1	prep 1 meas 0	0.00027386	0.00152347	1.24961e-03	(8.987 ± 6.248)e-4	IN (+6.07420e-04)
2	prep 0 excite to 1	0.00124865	0.00124973	1.07234e-06	(1.249 ± 0.001)e-3	IN (-5.21213e-07)
3	prep 1 decay to 0	0.00210546	0.00210653	1.07234e-06	(2.106 ± 0.001)e-3	IN (+5.21213e-07)

================================================================================
Valid Intervals for Matrix Entries
================================================================================

	Entry	Min	Max	Width	Value	Ref1
0	M^0[0,0] = p_0^(0,0)	0.99871989	0.99872034	4.49770e-07	(9.987 ± 0.000)e-1	IN (+2.01820e-07)
1	M^0[0,1] = p_1^(0,0)	0.00009060	0.00134032	1.24973e-03	(7.155 ± 6.249)e-4	IN (+6.07453e-04)
2	M^0[1,0] = p_0^(0,1)	0.00000000	0.00124973	1.24973e-03	(6.249 ± 6.249)e-4	IN (-6.07453e-04)
3	M^0[1,1] = p_1^(0,1)	0.00018315	0.00018360	4.49770e-07	(1.834 ± 0.002)e-4	IN (-2.01820e-07)
4	M^1[0,0] = p_0^(1,0)	0.00003005	0.00003106	1.01727e-06	(3.056 ± 0.051)e-5	IN (+4.58425e-07)
5	M^1[0,1] = p_1^(1,0)	0.00076621	0.00201486	1.24865e-03	(1.391 ± 0.624)e-3	IN (-6.06932e-04)
6	M^1[1,0] = p_0^(1,1)	-0.00000000	0.00124865	1.24865e-03	(6.243 ± 6.243)e-4	IN (+6.06932e-04)
7	M^1[1,1] = p_1^(1,1)	0.99771026	0.99771128	1.01727e-06	(9.977 ± 0.000)e-1	IN (-4.58425e-07)

Note: Reference comparison format:
  'IN (±X.XXXe±YY)' - value is within range, absolute discrepancy from center
  'OUT above/below (±X.XXXe±YY)' - value is outside range, absolute discrepancy from center
================================================================================



/var/folders/dk/_dd6ng8n3yv_crq6vrdt8qzc0000gn/T/ipykernel_62248/724936655.py:540: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()

Original Instrument:

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.99872032, 0.00132291]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00001741, 0.00018317]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00003101, 0.00078360]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00123126, 0.99771031]]

Reconstructed Instrument from center points of valid entry ranges:

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.99872012, 0.00071546]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00062486, 0.00018337]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00003056, 0.00139053]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00062433, 0.99771077]]

🔴 Main

selected_dir = "simulated_data_406_twirl"
selected_dir = "simulated_data_405"

GT_instrument = random_instrument(corr_strength=0.08, seed=2, fidelity=0.97)
# Generate and cache simulated data
experiment_info, v0_list, empirical_probs_list = generate_and_cache_simulated_data(
    GT=GT_instrument,
    num_seeds=50,
    shots_per_seed=int(10_000_000 / 0.1),
    L=4,
    use_same_v0_for_all_seeds=np.array([0.2, 0.8]),
    seed_for_reproduce=123,
    chosen_dir=selected_dir
)
for key in experiment_info.keys():
    print(key)
    # print(f"{key}: {experiment_info[key]}")

# GT_instrument = experiment_info["GT"]
GT_instrument = Instrument2x2(M0=experiment_info["GT"].P0, M1=experiment_info["GT"].P1)
GT_instrument.reveal()


GT_inv = summarize_instrument(GT_instrument.M0, GT_instrument.M1)
# print(GT_inv)
emp_inv = derived_constraints_from_empirical_probs(calculate_average_probs_np(empirical_probs_list))
# print(emp_inv)


sols = reconstruct_instrument_from_invariants_mixed_det(
    trM0 = emp_inv[0],
    detM0 = emp_inv[1],
    trM1 = emp_inv[2],
    detM1 = emp_inv[3],
    S0 = emp_inv[6] * 2, # S0 = sum of all elements of M0 = 2 * probability of observing 0 from v0=(0.5,0.5)
    gauge_p00 = 1.0
)

# # cheating
# sols = reconstruct_instrument_from_invariants_mixed_det(
#     trM0 = GT_instrument.M0.trace(),
#     detM0 = np.linalg.det(GT_instrument.M0),
#     trM1 = GT_instrument.M1.trace(),
#     detM1 = np.linalg.det(GT_instrument.M1),
#     S0 = GT_instrument.M0.sum(),
#     gauge_p00 = 1.0
# )

print(f"\nthere are {len(sols)} solutions from empirical data.")

# display a table to compare the args used for reconstruction from both GT and empirical data
comparison_table = pd.DataFrame({
    "Parameter": ["trM0", "detM0", "trM1", "detM1", "S0"],
    "From GT Instrument": [GT_instrument.M0.trace(), np.linalg.det(GT_instrument.M0), GT_instrument.M1.trace(), np.linalg.det(GT_instrument.M1), GT_instrument.M0.sum()],
    "From Empirical Data": [emp_inv[0], emp_inv[1], emp_inv[2], emp_inv[3], emp_inv[6] * 2],
})
print("\nComparison of parameters used for reconstruction:")
print(comparison_table)

Loading cached simulated data from simulated_data_405...
Override experiment info with cached version.
GT
number_of_seeds
shots_per_seed
max_string_length_L
use_same_v0_for_all_seeds
seed_for_reproduce

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.96510057, 0.00054589]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.02642910, 0.00642704]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00831841, 0.00260395]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00015192, 0.99042312]]

there are 2 solutions from empirical data.

Comparison of parameters used for reconstruction:
  Parameter  From GT Instrument  From Empirical Data
0      trM0            0.971528             0.971565
1     detM0            0.006188             0.006208
2      trM1            0.998742             0.998695
3     detM1            0.008238             0.008211
4        S0            0.998503             0.998559

chosen_sol_index = 0  # Choose which solution to analyze

reconstructed_MCM, _ = plot_gauge_transformation_effects(
    Instrument2x2(M0=sols[chosen_sol_index][0], M1=sols[chosen_sol_index][1]), 
    t_width_factor=1.1, 
    verbose=False,
    MCM_reference=[GT_instrument],
    resolution=1e5
    )

print("Original Instrument:")
GT_instrument.reveal()

print("\nReconstructed Instrument from invariants (empirical data):")
if reconstructed_MCM is not None:
    reconstructed_MCM.reveal()

================================================================================
Valid Intervals for Derived Quantities (Readout & Back-action Errors)
================================================================================

	Quantity	Min	Max	Width	Value	Ref1
0	prep 0 meas 1	0.00829104	0.01103694	2.74591e-03	(9.664 ± 1.373)e-3	IN (-1.19366e-03)
1	prep 1 meas 0	0.00684990	0.00959581	2.74591e-03	(8.223 ± 1.373)e-3	IN (-1.24993e-03)
2	prep 0 excite to 1	0.02658707	0.02665277	6.57042e-05	(2.662 ± 0.003)e-2	OUT below (-3.88947e-05)
3	prep 1 decay to 0	0.00308699	0.00315270	6.57042e-05	(3.120 ± 0.033)e-3	IN (+2.99946e-05)

================================================================================
Valid Intervals for Matrix Entries
================================================================================

	Entry	Min	Max	Width	Value	Ref1
0	M^0[0,0] = p_0^(0,0)	0.96505619	0.96512190	6.57042e-05	(9.651 ± 0.000)e-1	IN (+1.15256e-05)
1	M^0[0,1] = p_1^(0,0)	0.00040679	0.00308699	2.68020e-03	(1.747 ± 1.340)e-3	IN (-1.20100e-03)
2	M^0[1,0] = p_0^(0,1)	0.02390687	0.02658707	2.68020e-03	(2.525 ± 0.134)e-2	IN (+1.18214e-03)
3	M^0[1,1] = p_1^(0,1)	0.00644311	0.00650882	6.57042e-05	(6.476 ± 0.033)e-3	OUT below (-4.89266e-05)
4	M^1[0,0] = p_0^(1,0)	0.00829104	0.00829296	1.91931e-06	(8.292 ± 0.001)e-3	OUT above (+2.64095e-05)
5	M^1[0,1] = p_1^(1,0)	0.00000000	0.00274591	2.74591e-03	(1.373 ± 1.373)e-3	IN (+1.23100e-03)
6	M^1[1,0] = p_0^(1,1)	-0.00000000	0.00274591	2.74591e-03	(1.373 ± 1.373)e-3	IN (-1.22103e-03)
7	M^1[1,1] = p_1^(1,1)	0.99040227	0.99040419	1.91931e-06	(9.904 ± 0.000)e-1	OUT above (+1.98916e-05)

Note: Reference comparison format:
  'IN (±X.XXXe±YY)' - value is within range, absolute discrepancy from center
  'OUT above/below (±X.XXXe±YY)' - value is outside range, absolute discrepancy from center
================================================================================



/var/folders/dk/_dd6ng8n3yv_crq6vrdt8qzc0000gn/T/ipykernel_62248/724936655.py:540: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()

Original Instrument:

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.96510057, 0.00054589]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.02642910, 0.00642704]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00831841, 0.00260395]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00015192, 0.99042312]]

Reconstructed Instrument from invariants (empirical data):

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.96508904, 0.00174689]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.02524697, 0.00647597]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00829200, 0.00137295]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00137295, 0.99040323]]

regions = allowed_t_regions_for_list([sols[0][0], sols[0][1]], tol=1e-24)
regions

[(np.float64(0.14868186307103115), np.float64(0.15064089805685127)),
 (np.float64(0.8493591019431488), np.float64(0.8513181369289688))]

# this is demo of if we start from the "cheated" solution using true GT directly, so you will see that the RMSE is very small and happens at t=0

GT_instrument = Instrument2x2(M0=experiment_info["GT"].P0, M1=experiment_info["GT"].P1)
reconstructed_MCM, _ = plot_gauge_transformation_effects(
    GT_instrument,
    t_width_factor=1.1, 
    verbose=False,
    MCM_reference=[GT_instrument],
    resolution=1e5
    )

print("Original Instrument:")
GT_instrument.reveal()

print("Reconstructed Instrument from invariants (best match to original):")
if reconstructed_MCM is not None:
    reconstructed_MCM.reveal()

================================================================================
Valid Intervals for Derived Quantities (Readout & Back-action Errors)
================================================================================

	Quantity	Min	Max	Width	Value	Ref1
0	prep 0 meas 1	0.00831800	0.01107388	2.75587e-03	(9.696 ± 1.378)e-3	IN (-1.22561e-03)
1	prep 1 meas 0	0.00682060	0.00957648	2.75587e-03	(8.199 ± 1.378)e-3	IN (-1.22561e-03)
2	prep 0 excite to 1	0.02657740	0.02664332	6.59146e-05	(2.661 ± 0.003)e-2	IN (-2.93333e-05)
3	prep 1 decay to 0	0.00308755	0.00315346	6.59146e-05	(3.121 ± 0.033)e-3	IN (+2.93333e-05)

================================================================================
Valid Intervals for Matrix Entries
================================================================================

	Entry	Min	Max	Width	Value	Ref1
0	M^0[0,0] = p_0^(0,0)	0.96503868	0.96510460	6.59146e-05	(9.651 ± 0.000)e-1	IN (+2.89305e-05)
1	M^0[0,1] = p_1^(0,0)	0.00039759	0.00308755	2.68996e-03	(1.743 ± 1.345)e-3	IN (-1.19668e-03)
2	M^0[1,0] = p_0^(0,1)	0.02388744	0.02657740	2.68996e-03	(2.523 ± 0.134)e-2	IN (+1.19668e-03)
3	M^0[1,1] = p_1^(0,1)	0.00642301	0.00648893	6.59146e-05	(6.456 ± 0.033)e-3	IN (-2.89305e-05)
4	M^1[0,0] = p_0^(1,0)	0.00831800	0.00831994	1.93330e-06	(8.319 ± 0.001)e-3	IN (-5.63847e-07)
5	M^1[0,1] = p_1^(1,0)	0.00000000	0.00275587	2.75587e-03	(1.378 ± 1.378)e-3	IN (+1.22601e-03)
6	M^1[1,0] = p_0^(1,1)	0.00000000	0.00275587	2.75587e-03	(1.378 ± 1.378)e-3	IN (-1.22601e-03)
7	M^1[1,1] = p_1^(1,1)	0.99042159	0.99042352	1.93330e-06	(9.904 ± 0.000)e-1	IN (+5.63847e-07)

Note: Reference comparison format:
  'IN (±X.XXXe±YY)' - value is within range, absolute discrepancy from center
  'OUT above/below (±X.XXXe±YY)' - value is outside range, absolute discrepancy from center
================================================================================



/var/folders/dk/_dd6ng8n3yv_crq6vrdt8qzc0000gn/T/ipykernel_62248/724936655.py:540: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()

Original Instrument:

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.96510057, 0.00054589]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.02642910, 0.00642704]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00831841, 0.00260395]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00015192, 0.99042312]]
Reconstructed Instrument from invariants (best match to original):

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.96507164, 0.00174257]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.02523242, 0.00645597]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00831897, 0.00137794]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00137794, 0.99042256]]

import sympy as sp

# Define symbolic variables
a, b, c, d, t = sp.symbols('a b c d t', real=True)

# Define the original matrix M
M = sp.Matrix([[a, b], [c, d]])

# Define the gauge transformation matrix R(t)
R = sp.Matrix([[1-t, t], [t, 1-t]])

# Compute R^(-1)
R_inv = R.inv()

# Compute the gauge-transformed matrix M' = R^(-1) @ M @ R
M_prime = R_inv @ M @ R

# Simplify the result
M_prime_simplified = sp.simplify(M_prime)

print("Original matrix M:")
display(M)

print("\nGauge transformation matrix R(t):")
display(R)

print("\nInverse R^(-1):")
display(sp.simplify(R_inv))

print("\nGauge-transformed matrix M' = R^(-1) M R (before simplification):")
display(M_prime)

print("\nGauge-transformed matrix M' = R^(-1) M R (after simplification):")
display(M_prime_simplified)

# Let's also expand and collect terms for each element
print("\nElement-wise simplified expressions:")
for i in range(2):
    for j in range(2):
        element = sp.simplify(M_prime_simplified[i, j])
        element_expanded = sp.expand(element)
        print(f"M'[{i},{j}] = {element_expanded}")

Original matrix M:

\(\displaystyle \left[\begin{matrix}a & b\\c & d\end{matrix}\right]\)

Gauge transformation matrix R(t):

\(\displaystyle \left[\begin{matrix}1 - t & t\\t & 1 - t\end{matrix}\right]\)

Inverse R^(-1):

\(\displaystyle \left[\begin{matrix}\frac{t - 1}{2 t - 1} & \frac{t}{2 t - 1}\\\frac{t}{2 t - 1} & \frac{t - 1}{2 t - 1}\end{matrix}\right]\)

Gauge-transformed matrix M' = R^(-1) M R (before simplification):

\(\displaystyle \left[\begin{matrix}t \left(\frac{b \left(t - 1\right)}{2 t - 1} + \frac{d t}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{a \left(t - 1\right)}{2 t - 1} + \frac{c t}{2 t - 1}\right) & t \left(\frac{a \left(t - 1\right)}{2 t - 1} + \frac{c t}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{b \left(t - 1\right)}{2 t - 1} + \frac{d t}{2 t - 1}\right)\\t \left(\frac{b t}{2 t - 1} + \frac{d \left(t - 1\right)}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{a t}{2 t - 1} + \frac{c \left(t - 1\right)}{2 t - 1}\right) & t \left(\frac{a t}{2 t - 1} + \frac{c \left(t - 1\right)}{2 t - 1}\right) + \left(1 - t\right) \left(\frac{b t}{2 t - 1} + \frac{d \left(t - 1\right)}{2 t - 1}\right)\end{matrix}\right]\)

Gauge-transformed matrix M' = R^(-1) M R (after simplification):

\(\displaystyle \left[\begin{matrix}\frac{t \left(b \left(t - 1\right) + d t\right) - \left(t - 1\right) \left(a \left(t - 1\right) + c t\right)}{2 t - 1} & \frac{t \left(a \left(t - 1\right) + c t\right) - \left(t - 1\right) \left(b \left(t - 1\right) + d t\right)}{2 t - 1}\\\frac{t \left(b t + d \left(t - 1\right)\right) - \left(t - 1\right) \left(a t + c \left(t - 1\right)\right)}{2 t - 1} & \frac{t \left(a t + c \left(t - 1\right)\right) - \left(t - 1\right) \left(b t + d \left(t - 1\right)\right)}{2 t - 1}\end{matrix}\right]\)

Element-wise simplified expressions:
M'[0,0] = -a*t**2/(2*t - 1) + 2*a*t/(2*t - 1) - a/(2*t - 1) + b*t**2/(2*t - 1) - b*t/(2*t - 1) - c*t**2/(2*t - 1) + c*t/(2*t - 1) + d*t**2/(2*t - 1)
M'[0,1] = a*t**2/(2*t - 1) - a*t/(2*t - 1) - b*t**2/(2*t - 1) + 2*b*t/(2*t - 1) - b/(2*t - 1) + c*t**2/(2*t - 1) - d*t**2/(2*t - 1) + d*t/(2*t - 1)
M'[1,0] = -a*t**2/(2*t - 1) + a*t/(2*t - 1) + b*t**2/(2*t - 1) - c*t**2/(2*t - 1) + 2*c*t/(2*t - 1) - c/(2*t - 1) + d*t**2/(2*t - 1) - d*t/(2*t - 1)
M'[1,1] = a*t**2/(2*t - 1) - b*t**2/(2*t - 1) + b*t/(2*t - 1) + c*t**2/(2*t - 1) - c*t/(2*t - 1) - d*t**2/(2*t - 1) + 2*d*t/(2*t - 1) - d/(2*t - 1)

# Express the numerators as polynomials in t of degree 2
print("="*70)
print("Numerators of M0' elements as polynomials in t (degree 2):")
print("="*70)

for i in range(2):
    for j in range(2):
        element = simplified_m0_prime[i, j]

        # Get numerator
        numerator = sp.numer(element)

        # Expand and collect terms by powers of t
        numerator_expanded = sp.expand(numerator)
        numerator_poly = sp.collect(numerator_expanded, t)

        # Convert to Poly object to extract coefficients
        poly = sp.Poly(numerator_poly, t)
        coeffs = poly.all_coeffs()

        # Ensure we have 3 coefficients (degree 2)
        while len(coeffs) < 3:
            coeffs.append(sp.sympify(0))

        print(f"\nM0'[{i},{j}] numerator:")
        print(f"  Polynomial: ({coeffs[0]})*t^2 + ({coeffs[1]})*t + ({coeffs[2]})")
        print(f"  Expanded form: {numerator_poly}")

print("\n" + "="*70)
print("Numerators of M1' elements as polynomials in t (degree 2):")
print("="*70)

for i in range(2):
    for j in range(2):
        element = simplified_m1_prime[i, j]

        # Get numerator
        numerator = sp.numer(element)

        # Expand and collect terms by powers of t
        numerator_expanded = sp.expand(numerator)
        numerator_poly = sp.collect(numerator_expanded, t)

        # Convert to Poly object to extract coefficients
        poly = sp.Poly(numerator_poly, t)
        coeffs = poly.all_coeffs()

        # Ensure we have 3 coefficients (degree 2)
        while len(coeffs) < 3:
            coeffs.append(sp.sympify(0))

        print(f"\nM1'[{i},{j}] numerator:")
        print(f"  Polynomial: ({coeffs[0]})*t^2 + ({coeffs[1]})*t + ({coeffs[2]})")
        print(f"  Expanded form: {numerator_poly}")

print("\n" + "="*70)
print("Note: All elements have the common denominator (2*t - 1)")
print("="*70)

======================================================================
Numerators of M0' elements as polynomials in t (degree 2):
======================================================================

M0'[0,0] numerator:
  Polynomial: (-a - b + e + f)*t^2 + (2*a + b - e)*t + (-a)
  Expanded form: -a + t**2*(-a - b + e + f) + t*(2*a + b - e)

M0'[0,1] numerator:
  Polynomial: (a + b - e - f)*t^2 + (-a + 2*e + f)*t + (-e)
  Expanded form: -e + t**2*(a + b - e - f) + t*(-a + 2*e + f)

M0'[1,0] numerator:
  Polynomial: (-a - b + e + f)*t^2 + (a + 2*b - f)*t + (-b)
  Expanded form: -b + t**2*(-a - b + e + f) + t*(a + 2*b - f)

M0'[1,1] numerator:
  Polynomial: (a + b - e - f)*t^2 + (-b + e + 2*f)*t + (-f)
  Expanded form: -f + t**2*(a + b - e - f) + t*(-b + e + 2*f)

======================================================================
Numerators of M1' elements as polynomials in t (degree 2):
======================================================================

M1'[0,0] numerator:
  Polynomial: (a + b - e - f)*t^2 + (-a - b + c - g + 1)*t + (-c)
  Expanded form: -c + t**2*(a + b - e - f) + t*(-a - b + c - g + 1)

M1'[0,1] numerator:
  Polynomial: (-a - b + e + f)*t^2 + (-c - e - f + g + 1)*t + (-g)
  Expanded form: -g + t**2*(-a - b + e + f) + t*(-c - e - f + g + 1)

M1'[1,0] numerator:
  Polynomial: (a + b - e - f)*t^2 + (-2*a - 2*b - c + e + f + g + 1)*t + (a + b + c - 1)
  Expanded form: a + b + c + t**2*(a + b - e - f) + t*(-2*a - 2*b - c + e + f + g + 1) - 1

M1'[1,1] numerator:
  Polynomial: (-a - b + e + f)*t^2 + (a + b + c - 2*e - 2*f - g + 1)*t + (e + f + g - 1)
  Expanded form: e + f + g + t**2*(-a - b + e + f) + t*(a + b + c - 2*e - 2*f - g + 1) - 1

======================================================================
Note: All elements have the common denominator (2*t - 1)
======================================================================

# Extract numerators and denominators for each element of M'
print("Numerators of M'_simplified elements as polynomials in t:\n")
print("="*70)

for i in range(2):
    for j in range(2):
        element = simplified_m0_prime[i, j]

        # Get numerator and denominator
        numerator = sp.numer(element)
        denominator = sp.denom(element)

        # Expand and collect terms by powers of t
        numerator_expanded = sp.expand(numerator)
        numerator_poly = sp.collect(numerator_expanded, t)

        print(f"\nM0'[{i},{j}]:")
        print(f"  Numerator: {numerator_poly}")
        print(f"  Denominator: {denominator}")

        # Extract coefficients
        poly = sp.Poly(numerator_poly, t)
        coeffs = poly.all_coeffs()
        print(f"  Polynomial form: ", end="")
        for idx, coeff in enumerate(coeffs):
            power = len(coeffs) - 1 - idx
            if power == 0:
                print(f"({coeff})", end="")
            elif power == 1:
                print(f"({coeff})*t + ", end="")
            else:
                print(f"({coeff})*t^{power} + ", end="")
        print()

print("\n" + "="*70)
print("\nM1 matrix elements:\n")
print("="*70)

for i in range(2):
    for j in range(2):
        element = simplified_m1_prime[i, j]

        # Get numerator and denominator
        numerator = sp.numer(element)
        denominator = sp.denom(element)

        # Expand and collect terms by powers of t
        numerator_expanded = sp.expand(numerator)
        numerator_poly = sp.collect(numerator_expanded, t)

        print(f"\nM1'[{i},{j}]:")
        print(f"  Numerator: {numerator_poly}")
        print(f"  Denominator: {denominator}")

        # Extract coefficients
        poly = sp.Poly(numerator_poly, t)
        coeffs = poly.all_coeffs()
        print(f"  Polynomial form: ", end="")
        for idx, coeff in enumerate(coeffs):
            power = len(coeffs) - 1 - idx
            if power == 0:
                print(f"({coeff})", end="")
            elif power == 1:
                print(f"({coeff})*t + ", end="")
            else:
                print(f"({coeff})*t^{power} + ", end="")
        print()

# Summary: All elements have the common denominator (2*t - 1)
print("\n" + "="*70)
print("Note: All elements share the common denominator (2*t - 1)")
print("="*70)

Numerators of M'_simplified elements as polynomials in t:

======================================================================

M0'[0,0]:
  Numerator: -a + t**2*(-a - b + e + f) + t*(2*a + b - e)
  Denominator: 2*t - 1
  Polynomial form: (-a - b + e + f)*t^2 + (2*a + b - e)*t + (-a)

M0'[0,1]:
  Numerator: -e + t**2*(a + b - e - f) + t*(-a + 2*e + f)
  Denominator: 2*t - 1
  Polynomial form: (a + b - e - f)*t^2 + (-a + 2*e + f)*t + (-e)

M0'[1,0]:
  Numerator: -b + t**2*(-a - b + e + f) + t*(a + 2*b - f)
  Denominator: 2*t - 1
  Polynomial form: (-a - b + e + f)*t^2 + (a + 2*b - f)*t + (-b)

M0'[1,1]:
  Numerator: -f + t**2*(a + b - e - f) + t*(-b + e + 2*f)
  Denominator: 2*t - 1
  Polynomial form: (a + b - e - f)*t^2 + (-b + e + 2*f)*t + (-f)

======================================================================

M1 matrix elements:

======================================================================

M1'[0,0]:
  Numerator: -c + t**2*(a + b - e - f) + t*(-a - b + c - g + 1)
  Denominator: 2*t - 1
  Polynomial form: (a + b - e - f)*t^2 + (-a - b + c - g + 1)*t + (-c)

M1'[0,1]:
  Numerator: -g + t**2*(-a - b + e + f) + t*(-c - e - f + g + 1)
  Denominator: 2*t - 1
  Polynomial form: (-a - b + e + f)*t^2 + (-c - e - f + g + 1)*t + (-g)

M1'[1,0]:
  Numerator: a + b + c + t**2*(a + b - e - f) + t*(-2*a - 2*b - c + e + f + g + 1) - 1
  Denominator: 2*t - 1
  Polynomial form: (a + b - e - f)*t^2 + (-2*a - 2*b - c + e + f + g + 1)*t + (a + b + c - 1)

M1'[1,1]:
  Numerator: e + f + g + t**2*(-a - b + e + f) + t*(a + b + c - 2*e - 2*f - g + 1) - 1
  Denominator: 2*t - 1
  Polynomial form: (-a - b + e + f)*t^2 + (a + b + c - 2*e - 2*f - g + 1)*t + (e + f + g - 1)

======================================================================
Note: All elements share the common denominator (2*t - 1)
======================================================================

🖥️ Qiskit Experiments

Utilities for Qiskit Experiments

from qiskit import __version__

print(__version__)

2.2.3

from qiskit_ibm_runtime import QiskitRuntimeService

# # Save credentials locally (can be run once)
# QiskitRuntimeService.save_account(
#     token="E3zZk7CbDNC3EeiR8Cpr-UcI4gKPOPHo425CRKtJIEK6", # Use the 44-character API_KEY you created and saved from the IBM Quantum Platform Home dashboard
#     instance="crn:v1:bluemix:public:quantum-computing:us-east:a/ed5d7d2fb3b249c6baea6864058116cb:e46b4031-d4ee-48a6-88cd-1ae20a49ec21::", # Optional
# )

# Run every time you need the service
service = QiskitRuntimeService()

data = []
for backend in service.backends():
    config = backend.configuration()
    if "simulator" in config.backend_name:
        continue
    data.append({
        "Backend": config.backend_name,
        "Processor Type": config.processor_type,
        "Supported Instructions (measure_2 is mid-circ meas.)": ", ".join(config.supported_instructions)
    })

df_backends = pd.DataFrame(data)

# Set pandas option to display full column width
pd.set_option('display.max_colwidth', None)

display(df_backends)

	Backend	Processor Type	Supported Instructions (measure_2 is mid-circ meas.)
0	ibm_pittsburgh	{'family': 'Heron', 'revision': '3'}	cz, id, delay, measure, measure_2, reset, rz, sx, x, if_else
1	ibm_fez	{'family': 'Heron', 'revision': '2'}	cz, id, delay, measure, reset, rz, sx, x, if_else
2	ibm_marrakesh	{'family': 'Heron', 'revision': '2'}	cz, id, delay, measure, reset, rz, sx, x, if_else
3	ibm_torino	{'family': 'Heron', 'revision': '1'}	cz, id, delay, measure, reset, rz, sx, x, if_else
4	ibm_kingston	{'family': 'Heron', 'revision': '2'}	cz, id, delay, measure, measure_2, reset, rz, sx, x, if_else

from qiskit import QuantumCircuit
from qiskit import transpile
from qiskit.visualization import plot_histogram
from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from qiskit_ibm_runtime.circuit import MidCircuitMeasure
from qiskit.circuit import Measure
from qiskit.transpiler import generate_preset_pass_manager
from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime.fake_provider import FakeManilaV2, FakeFez
from qiskit_aer import AerSimulator
from qiskit_aer.noise import (
    NoiseModel,
    QuantumError,
    ReadoutError,
    depolarizing_error,
    pauli_error,
    thermal_relaxation_error,
)
from qiskit.quantum_info import Clifford, random_clifford
from qiskit.synthesis import OneQubitEulerDecomposer
import copy

# === Define registers and circuit ===

def create_list_of_circuits(max_word_length, num_qubits=2, visualize_first=4):

    circuits_for_local_test = []
    circuits_for_QPU = []
    qreg_q = QuantumRegister(num_qubits, 'q')
    # Create one classical register per qubit, each of size max_word_length
    cregs = [ClassicalRegister(max_word_length, f'c{i}') for i in range(num_qubits)]

    # Iterate over each qubit to create its specific set of circuits
    for q in range(num_qubits):
        # Iterate over the 4 Pauli basis preparations
        for gate_name in ['I', 'X', 'Y', 'Z']:

            # Create fresh circuits including all registers
            qc_local = QuantumCircuit(qreg_q, *cregs)
            qc_qpu = QuantumCircuit(qreg_q, *cregs)

            # Apply the gate to the current qubit 'q'
            if gate_name == 'X':
                qc_local.x(q)
                qc_qpu.x(q)
            elif gate_name == 'Y':
                qc_local.y(q)
                qc_qpu.y(q)
            elif gate_name == 'Z':
                qc_local.z(q)
                qc_qpu.z(q)
            # 'I' is identity

            # Add measurements to the current qubit 'q' into its corresponding classical register
            for i in range(max_word_length-1, -1, -1):
                # Local test uses standard Measure
                qc_local.append(Measure(), [q], [cregs[q][i]])
                # QPU uses MidCircuitMeasure
                qc_qpu.append(MidCircuitMeasure(), [q], [cregs[q][i]])

            circuits_for_local_test.append(qc_local)
            circuits_for_QPU.append(qc_qpu)

    # === Visualize the circuit ===
    print(f"Visualizing circuit for local test (first {visualize_first}):")
    for circ in circuits_for_local_test[:visualize_first]:
        circ.draw('mpl', style="iqp")

    print(f"Visualizing circuit for QPU (first {visualize_first}):")
    for circ in circuits_for_QPU[:visualize_first]:
        circ.draw('mpl', style="iqp")

    return circuits_for_local_test, circuits_for_QPU




def create_list_of_circuits_v2(max_word_length):  # refer to https://quantum.cloud.ibm.com/docs/en/api/qiskit/qiskit.synthesis.OneQubitEulerDecomposer
    circuits_for_local_test = []
    circuits_for_QPU = []
    qreg_q = QuantumRegister(1, 'q')
    creg_c = ClassicalRegister(max_word_length, 'c')

    # === 1. Generate the 24 Single-Qubit Clifford Unitaries ===
    # "Full twirling" requires averaging over the full Clifford group (size 24).
    # We generate them by sampling until we have all 24 unique elements.
    clifford_set = set()
    clifford_ops = []

    while len(clifford_ops) < 24:
        c = random_clifford(1)
        # We use the matrix representation to check for uniqueness and decomposition
        c_matrix = tuple(c.to_matrix().flatten()) 
        # Note: Phase differences matter for gates, but often we just need the group 
        # modulo phase for twirling. Here we keep distinct operations.
        # To be safe and simple, we collect unique Cliffords.
        if str(c) not in clifford_set:
            clifford_set.add(str(c))
            clifford_ops.append(c)

    # === 2. Decompose into Basis Gates ===
    # Your colleague recommended OneQubitEulerDecomposer.
    # We use 'ZSX' basis (common for IBM hardware) to decompose the Cliffords
    # into sequences of Rz and SX gates.
    decomposer = OneQubitEulerDecomposer(basis='ZSX')

    # Create the base circuits
    twirling_circuits = []
    for c_op in clifford_ops:
        # Convert the Clifford operator to a unitary matrix
        unitary = c_op.to_matrix()

        # Synthesize the circuit from the unitary
        # simplify=True (default) will optimize Identity to empty. 
        # If you STRICTLY need a gate for Identity, you'd need to manually handle it,
        # but for standard twirling, the idle identity is correct.
        decomposed_qc = decomposer(unitary)

        # Create a fresh circuit container
        qc = QuantumCircuit(qreg_q, creg_c)
        qc.compose(decomposed_qc, qubits=[0], inplace=True)
        twirling_circuits.append(qc)

    # === 3. Expand lists for Local/QPU usage ===
    # Extend the lists with the generated 24 circuits
    circuits_for_local_test.extend(twirling_circuits)
    circuits_for_QPU.extend(twirling_circuits)

    circuits_for_local_test = copy.deepcopy(circuits_for_local_test)
    circuits_for_QPU = copy.deepcopy(circuits_for_QPU)

    # === 4. Add Measurements ===
    for i in range(max_word_length-1, -1, -1):
        for circ in circuits_for_local_test:
            circ.append(Measure(), [0], [i])
        for circ in circuits_for_QPU:
            circ.append(MidCircuitMeasure(), [0], [i])

    # === Visualize ===
    print(f"Generated {len(circuits_for_local_test)} circuits (Clifford Group Size).")
    print("Visualizing a few sample circuits:")
    # Draw the first few distinct circuits to verify decomposition
    for k in range(min(24, len(circuits_for_local_test))):
        print(f"Circuit {k}:")
        circuits_for_local_test[k].draw('mpl', style="iqp")

    return circuits_for_local_test, circuits_for_QPU

def create_non_trivial_list_of_circuits(family_names = ['repeated-X-MCM'], max_word_length=4, num_qubits=2):
    """
    Generates a list of circuits based on specific non-trivial families of MCM sequences.

    Args:
        family_names: List of strings specifying which circuit families to generate.
                    Currently supports: 'repeated-X-MCM'.
        max_word_length: Length of the classical register (number of measurements).
        num_qubits: Number of qubits to generate circuits for.

    Returns:
        Tuple of (circuits_for_local_test, circuits_for_QPU)
    """
    circuits_for_local_test = []
    circuits_for_QPU = []

    qreg_q = QuantumRegister(num_qubits, 'q')
    # Create one classical register per qubit, each of size max_word_length
    cregs = [ClassicalRegister(max_word_length, f'c{i}') for i in range(num_qubits)]

    for family in family_names:
        if family == 'repeated-X-MCM':
            # Iterate over each qubit to create its specific set of circuits
            for q in range(num_qubits):
                qc_local = QuantumCircuit(qreg_q, *cregs)
                qc_qpu = QuantumCircuit(qreg_q, *cregs)

                # Apply X gate then Measure, repeated max_word_length times
                # This creates a toggling sequence: |0> -> X -> |1> -> Meas(1) -> X -> |0> -> Meas(0) ...
                for i in range(max_word_length-1, -1, -1):
                    # Apply X gate
                    qc_local.x(q)
                    qc_qpu.x(q)

                    # Measurement
                    # Local test uses standard Measure
                    qc_local.append(Measure(), [q], [cregs[q][i]])
                    # QPU uses MidCircuitMeasure
                    qc_qpu.append(MidCircuitMeasure(), [q], [cregs[q][i]])

                circuits_for_local_test.append(qc_local)
                circuits_for_QPU.append(qc_qpu)
        else:
            print(f"Warning: Family '{family}' not implemented.")

    return circuits_for_local_test, circuits_for_QPU

def plot_job_results(job_result, num_qubits, backend_name=None):
    """
    Plot histogram comparing results across multiple circuits (PUBs), grouped by classical register index.

    Args:
        job_result: PrimitiveResult object from Qiskit Runtime
        num_qubits: Number of qubits (groups) to split the results into.
        backend_name: Optional name of the backend used for execution

    Returns:
        list: A list of length num_qubits, where each element is a list of count dictionaries 
              for that classical register index.
    """
    # Initialize list of lists to hold counts for each qubit/classical register
    results_by_qubit = [[] for _ in range(num_qubits)]

    # Calculate how many PUBs correspond to each qubit group
    pubs_per_qubit = len(job_result) / num_qubits

    # Collect all bitstrings and counts from each PUB, grouped by register index
    for i, pub_result in enumerate(job_result):
        # Determine which classical register index this PUB targets
        c_reg_index = int(i // pubs_per_qubit)

        # Extract the specific classical register data using the provided logic
        # We iterate through keys to find the one at index c_reg_index
        for j, key in enumerate(pub_result.data.keys()):
            if j == c_reg_index:
                counts = pub_result.data[key].get_counts()
                results_by_qubit[c_reg_index].append(counts)
                print(f"\nCounts for classical register index {c_reg_index} from PUB {i}:")
                print(counts)
                break

    # Plot histograms for each qubit group
    for q_idx, pub_counts_list in enumerate(results_by_qubit):
        if not pub_counts_list:
            continue

        all_bitstrings = set()
        for counts in pub_counts_list:
            all_bitstrings.update(counts.keys())

        # Sort bitstrings for consistent ordering
        sorted_bitstrings = sorted(all_bitstrings)

        # Prepare data for grouped bar chart
        num_pubs_in_group = len(pub_counts_list)
        colors = ['blue', 'orange', 'green', 'red', 'purple', 'brown', 'cyan', 'magenta']
        # Ensure we have enough colors by cycling if necessary
        colors = [colors[k % len(colors)] for k in range(num_pubs_in_group)]

        bar_width = 0.8 / num_pubs_in_group
        x = np.arange(len(sorted_bitstrings))

        # Create plot for this qubit group
        fig, ax = plt.subplots(figsize=(10, 6))

        # Plot bars for each PUB in this group
        for k, counts in enumerate(pub_counts_list):
            values = [counts.get(bitstring, 0) for bitstring in sorted_bitstrings]
            offset = (k - num_pubs_in_group/2 + 0.5) * bar_width
            bars = ax.bar(x + offset, values, bar_width, 
                           label=f'Circuit {k}', color=colors[k])

            # Add count labels above each bar
            for bar, val in zip(bars, values):
                if val > 0:  # Only show label if count > 0
                    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 50,
                           str(val), ha='center', va='bottom', fontsize=8, rotation=90)

        ax.set_xlabel('Bitstring')
        ax.set_ylabel('Counts')
        ax.set_title(f'Comparison of Results for Classical Register Index {q_idx}')
        ax.set_xticks(x)
        ax.set_xticklabels(sorted_bitstrings, rotation=45)
        ax.legend()
        ax.grid(True, alpha=0.3, axis='y')

        plt.tight_layout()
        plt.show()

    if backend_name:
        print(f"\nReal hardware run on {backend_name}")

    return results_by_qubit

def counts_to_probabilities(
    counts_input: Dict[str, int] | List[Dict[str, int]],
    max_len: int = None
) -> Dict[str, float]:
    """
    Convert count dictionary/dictionaries to a single probability dictionary.

    This function takes either a single dictionary or a list of dictionaries containing
    bit-string counts, aggregates all counts, and normalizes them to probabilities.
    Probabilities for shorter strings are computed by marginalizing (summing) over
    the longer strings that extend them.

    Args:
        counts_input: Either a single dict mapping bit-strings to counts, or a list of such dicts.
        max_len: Maximum length of binary strings to include. If None, inferred from input data.

    Returns:
        Dictionary mapping all binary strings (up to max_len) to their probabilities.
        Shorter strings' probabilities are computed by summing over extensions.

    Example:
        >>> counts1 = {'00': 80, '01': 20}
        >>> counts2 = {'00': 70, '10': 30}
        >>> counts_to_probabilities([counts1, counts2])
        {'0': 0.75, '1': 0.25, '00': 0.75, '01': 0.1, '10': 0.15, '11': 0.0}
    """
    # Normalize input to list of dictionaries
    if isinstance(counts_input, dict):
        counts_list = [counts_input]
    else:
        counts_list = counts_input

    if not counts_list:
        return {}

    # Determine max_len if not provided
    if max_len is None:
        max_len = max(len(bitstring) for counts_dict in counts_list 
                     for bitstring in counts_dict.keys())

    # Aggregate all counts
    aggregated_counts = {}
    total_counts = 0

    for counts_dict in counts_list:
        for bitstring, count in counts_dict.items():
            if bitstring in aggregated_counts:
                aggregated_counts[bitstring] += count
            else:
                aggregated_counts[bitstring] = count
            total_counts += count

    # Initialize probability dictionary
    prob_dict = {}

    # First, normalize the max_len strings
    for bitstring in get_all_binary_strings(max_len):
        if bitstring in aggregated_counts:
            prob_dict[bitstring] = aggregated_counts[bitstring] / total_counts
        else:
            prob_dict[bitstring] = 0.0

    # Now compute probabilities for shorter strings by marginalization
    for length in range(max_len - 1, 0, -1):
        for bitstring in get_all_binary_strings(length):
            # Sum probabilities of all extensions (bitstring + '0' and bitstring + '1')
            prob = 0.0
            for bit in ['0', '1']:
                extended = bitstring + bit
                if extended in prob_dict:
                    prob += prob_dict[extended]
            prob_dict[bitstring] = prob

    return prob_dict

def compare_instrument_to_empirical(
    instrument: Instrument2x2,
    empirical_probs: Dict[str, float],
    v0: np.ndarray = np.array([0.5, 0.5]),
    max_len: int = 4
) -> Tuple[Dict[str, float], float, pd.DataFrame]:
    """
    Calculate theoretical probabilities from an instrument and compare to empirical data.

    Args:
        instrument: Instrument2x2 object to evaluate
        empirical_probs: Dictionary of empirical probabilities from hardware
        v0: Initial state distribution (default: maximally mixed [0.5, 0.5])
        max_len: Maximum length of binary strings to consider

    Returns:
        Tuple of:
            - Dictionary of theoretical probabilities
            - Overall RMSE between theoretical and empirical
            - DataFrame with detailed comparison
    """
    # Calculate theoretical probabilities
    theoretical_probs = calculate_exact_all_string_probabilities_from_v0_and_instrument(
        inst=instrument,
        v0=v0,
        max_len=max_len
    )

    # Get common keys (strings present in both dicts)
    common_keys = sorted(set(theoretical_probs.keys()) & set(empirical_probs.keys()))

    # Calculate differences
    comparison_data = []
    squared_errors = []

    for key in common_keys:
        theo_val = theoretical_probs[key]
        emp_val = empirical_probs[key]
        diff = theo_val - emp_val
        abs_diff = abs(diff)
        rel_error = abs_diff / emp_val if emp_val != 0 else np.inf

        squared_errors.append(diff**2)

        comparison_data.append({
            'String': key,
            'Length': len(key),
            'Theoretical': f"{theo_val:.8f}",
            'Empirical (ibm)': f"{emp_val:.8f}",
            'Difference': f"{diff:+.6e}",
            'Abs Diff': f"{abs_diff:.6e}",
            'Rel Error': f"{rel_error:.4%}" if rel_error != np.inf else "inf"
        })

    # Calculate RMSE
    rmse = np.sqrt(np.mean(squared_errors))

    # Create DataFrame
    df_comparison = pd.DataFrame(comparison_data)

    return theoretical_probs, rmse, df_comparison

def display_job_details(job):
    """
    Display detailed information about a Qiskit Runtime job.

    Args:
        job: Qiskit Runtime Job object
    """
    print(f"Job ID: {job.job_id()}")
    print(f"Backend: {job.backend()}")
    print(f"Status: {job.status()}")
    print(f"Creation Time: {job.creation_date}")
    print(f"Tags: {job.tags}")
    print(f"Metadata: {job.result().metadata['execution']}")

Check Usage (quota) left

# print("Available to run on:")
# display(service.backends())

# print("\n✅ But the ones equipped with mid-circuit measurement are:")
# display(service.backends(filters=lambda b: "measure_2" in b.supported_instructions))

# print("Service Usage:")
# display(service.usage())

When you initialize the Sampler, use the mode parameter to specify the mode you want it to run in. Possible values are batch, session, or backend objects for batch, session, and job execution mode, respectively. For more information, see Introduction to Qiskit Runtime execution modes. Note that Open Plan users cannot submit session jobs.

Create Jobs

should change the I and Z circuits according to the instruction in https://quantum.cloud.ibm.com/docs/en/api/qiskit/qiskit.synthesis.OneQubitEulerDecomposer

options = {
    "default_shots": int(1e5),
    # "simulator": {"seed_simulator": 42},
    }

max_word_length = 6

# roughly on Nov 26, 2025, the calibration data for ibm_kingston is:
chosen_qubits = [
    # 0,   # MEASURE_2 error: 6.616e-2, Pauli-X error: 2.447e-4
    # 5,   # MEASURE_2 error: 1.245e-2, Pauli-X error: 4.338e-4
    # 10,  # MEASURE_2 error: 2.197e-3, Pauli-X error: 2.371e-4
    # 15,  # MEASURE_2 error: 3.418e-3, Pauli-X error: 1.590e-4
    # 56,  # MEASURE_2 error: 5.859e-3, Pauli-X error: 5.146e-4
    # 64,  # MEASURE_2 error: 6.104e-3, Pauli-X error: 1.935e-4
    # 69,  # MEASURE_2 error: 3.906e-3, Pauli-X error: 1.882e-4
    # 74,  # MEASURE_2 error: 4.150e-3, Pauli-X error: 2.935e-4
    103, # MEASURE_2 error: 3.662e-2, Pauli-X error: 2.072e-4
    117, # MEASURE_2 error: 3.418e-3, Pauli-X error: 1.617e-4
    ]

num_qubits = len(chosen_qubits)

circuits_for_local, circuits_for_QPU = create_list_of_circuits(max_word_length, num_qubits, visualize_first=0)

accuracy_circuits_local, accuracy_circuits_QPU = create_non_trivial_list_of_circuits(
    family_names = ['repeated-X-MCM'],
    max_word_length=max_word_length,
    num_qubits=num_qubits
)

Visualizing circuit for local test (first 0):
Visualizing circuit for QPU (first 0):

# # noisy simulation on local simulator

# backend = FakeManilaV2()

# noise_model = NoiseModel.from_backend(backend)
# sim_noisy = AerSimulator(method="statevector", noise_model=noise_model)
# sim_noisy.set_options(shots=options["default_shots"])
# sim_result_noisy = sim_noisy.run(circuits_for_local).result()

# plot_histogram(sim_result_noisy.get_counts(), figsize=(14,6))
# print(f"noisy simulation run on backend '{backend.name}' results:")

Run Jobs

# backend = service.least_busy(operational=True, simulator=False, min_num_qubits=1)
backend = service.backend(
    # "ibm_torino", # Heron r1
    # "ibm_fez", # Heron r2
    "ibm_kingston", # Heron r2, supports MCM
    # "ibm_marrakesh", # Heron r2, many jobs pending
    # "ibm_pittsburgh", # Heron r3, supports MCM, but many jobs pending (i queue)
    )

# Transpile the circuits for execution
pm = generate_preset_pass_manager(backend=backend, optimization_level=0, initial_layout=chosen_qubits)

circuits_after_transpile = [pm.run(circuit) for circuit in circuits_for_QPU + accuracy_circuits_QPU] # change from `circuits_for_local` to `circuits_for_QPU` for QPU with MCM support

# # take a look at the transpiled circuits to run on QPU
# for circ in circuits_after_transpile:
#     circ.draw('mpl', style="iqp")

# 🔴 !! Only run this cell when ready to submit to real hardware

if "you_really_wanna_run" == "not_yet":

    print(f"🖥️ Running on {backend.name},\non qubit(s) with index: {chosen_qubits},\na total of {len(circuits_after_transpile)} circuits,\neach ciruit with {options['default_shots']} shots.")

    # Initialize the Sampler and Run
    sampler = Sampler(mode=backend, options=options)
    job = sampler.run(circuits_after_transpile)
    job_result = job.result()

    print("done 😃😃😃.")

Retrieve jobs

jobs_list = [
    'd4ekjg8lslhc73cuq4f0', # good for demo
    'd4el3rkcdebc73ev7eh0', # demo for low fidelity MCM
    'd4f1u5h2bisc73a22400', # margin_tol = 3e-4 might be tight enough
    'd4f21892bisc73a22740', # margin_tol = 5e-5 might be tight enough
    'd4f23rolslhc73cv8fhg', # margin_tol = 
    'd4f60fh2bisc73a268ag', # margin_tol = 0.0000126 will give a single-point valid t region
    'd4f78fglslhc73cvdk8g', # margin_tol = 0, 
    'd4f7g0glslhc73cvdsk0', # margin_tol = 0, a perfect reconstruction

    'd4f8vbccdebc73evs2r0', # margin_tol = 0.0000037, Kingston qubit 105, trial 1, MEASURE_2 error from ibm: 3.418e-3
    'd4f90aulo8as739oe4qg', # margin_tol = 0.00027, Kingston qubit 105, trial 2, seems bad and inconsistent with trial 1, maybe I should let it rest

    'd4f9570lslhc73cvffeg', # margin_tol = 0.00008, forgot which qubit

    'd4fb9692bisc73a2bhd0', # 1e6 shots for each of the 4 circuits on q48, perfect reconstruction with 0 margin_tol
    'd4fcer6lo8as739ohimg', # 1e5 shots for each of the 4 circuits on q48, perfect reconstruction with 0 margin_tol (trial 1)
    'd4fchbulo8as739ohkug', # 1e5 shots for each of the 4 circuits on q48, margin_tol = 0.00014 (trial 2)

    'd4fk7kccdebc73f08fqg', # 1e5 shots for each of the 4 circuits, on qubit #4, perfect reconstruction with 0 margin_tol
    'd4ii9iiv0j9c73e2cd3g', # 1e5 shots for each of the 4 circuits, on qubit #4, perfect reconstruction with 0 margin_tol (trial 2
]

i_wanna_retrieve_job_with_id = "d4j7bhh0i6jc73ddgcf0" # replace with some random below-5-digit string to skip retrieval

if len(i_wanna_retrieve_job_with_id) > 5:
    service = QiskitRuntimeService(
        channel='ibm_quantum_platform',
        instance='crn:v1:bluemix:public:quantum-computing:us-east:a/ed5d7d2fb3b249c6baea6864058116cb:e46b4031-d4ee-48a6-88cd-1ae20a49ec21::'
    )
    job = service.job(i_wanna_retrieve_job_with_id)
    job_result = job.result()
    print(f"Retrieved job result for job ID: {i_wanna_retrieve_job_with_id}\n")
    display_job_details(job)
else:
    print("No job ID specified for retrieval., proceeding with existing job_result (if exists).")

Retrieved job result for job ID: d4j7bhh0i6jc73ddgcf0

Job ID: d4j7bhh0i6jc73ddgcf0
Backend: <IBMBackend('ibm_kingston')>
Status: DONE
Creation Time: 2025-11-25 21:36:38.673134-06:00
Tags: ['1e6 shots', 'accuracy test', 'non-zero-tol for q105', 'non-zero-tol for q117', 'q117 tol=0.00028', 'q105 tol=0.000012', 'q105q117', "q117's result is accurate"]
Metadata: {'execution_spans': ExecutionSpans([DoubleSliceSpan(<start='2025-11-26 03:36:41', stop='2025-11-26 03:41:07', size=1000000>), DoubleSliceSpan(<start='2025-11-26 03:41:07', stop='2025-11-26 03:45:33', size=1000000>), DoubleSliceSpan(<start='2025-11-26 03:45:33', stop='2025-11-26 03:50:00', size=1000000>), DoubleSliceSpan(<start='2025-11-26 03:50:00', stop='2025-11-26 03:54:26', size=1000000>), DoubleSliceSpan(<start='2025-11-26 03:54:26', stop='2025-11-26 03:58:52', size=1000000>), DoubleSliceSpan(<start='2025-11-26 03:58:52', stop='2025-11-26 04:03:18', size=1000000>), DoubleSliceSpan(<start='2025-11-26 04:03:57', stop='2025-11-26 04:08:26', size=1000000>), DoubleSliceSpan(<start='2025-11-26 04:08:59', stop='2025-11-26 04:13:27', size=1000000>), DoubleSliceSpan(<start='2025-11-26 04:13:27', stop='2025-11-26 04:17:54', size=1000000>), DoubleSliceSpan(<start='2025-11-26 04:17:54', stop='2025-11-26 04:22:20', size=1000000>)])}

Reconstruction

display_job_details(job)

# plot histogram from the job result
counting_lists = plot_job_results(job_result[:-num_qubits*1], backend_name=backend.name, num_qubits=num_qubits)

Job ID: d4j7bhh0i6jc73ddgcf0
Backend: <IBMBackend('ibm_kingston')>
Status: DONE
Creation Time: 2025-11-25 21:36:38.673134-06:00
Tags: ['1e6 shots', 'accuracy test', 'non-zero-tol for q105', 'non-zero-tol for q117', 'q117 tol=0.00028', 'q105 tol=0.000012', 'q105q117', "q117's result is accurate"]
Metadata: {'execution_spans': ExecutionSpans([DoubleSliceSpan(<start='2025-11-26 03:36:41', stop='2025-11-26 03:41:07', size=1000000>), DoubleSliceSpan(<start='2025-11-26 03:41:07', stop='2025-11-26 03:45:33', size=1000000>), DoubleSliceSpan(<start='2025-11-26 03:45:33', stop='2025-11-26 03:50:00', size=1000000>), DoubleSliceSpan(<start='2025-11-26 03:50:00', stop='2025-11-26 03:54:26', size=1000000>), DoubleSliceSpan(<start='2025-11-26 03:54:26', stop='2025-11-26 03:58:52', size=1000000>), DoubleSliceSpan(<start='2025-11-26 03:58:52', stop='2025-11-26 04:03:18', size=1000000>), DoubleSliceSpan(<start='2025-11-26 04:03:57', stop='2025-11-26 04:08:26', size=1000000>), DoubleSliceSpan(<start='2025-11-26 04:08:59', stop='2025-11-26 04:13:27', size=1000000>), DoubleSliceSpan(<start='2025-11-26 04:13:27', stop='2025-11-26 04:17:54', size=1000000>), DoubleSliceSpan(<start='2025-11-26 04:17:54', stop='2025-11-26 04:22:20', size=1000000>)])}

Counts for classical register index 0 from PUB 0:
{'000000': 994687, '000001': 831, '111111': 380, '010000': 671, '100000': 733, '000100': 701, '000010': 685, '001111': 130, '001000': 742, '011000': 4, '000011': 118, '011111': 122, '000111': 135, '110000': 6, '101110': 2, '101111': 2, '001010': 3, '111101': 3, '111110': 2, '111011': 1, '011110': 4, '100111': 1, '000110': 7, '111100': 4, '010100': 4, '000101': 3, '111000': 5, '001110': 1, '001100': 4, '100100': 1, '100010': 1, '101000': 3, '010001': 1, '110111': 1, '010111': 1, '010010': 1}

Counts for classical register index 0 from PUB 1:
{'111111': 937186, '111000': 11458, '111100': 11431, '000000': 5101, '100000': 11756, '110000': 11585, '111110': 11211, '010000': 7, '110010': 9, '110111': 27, '111001': 10, '110100': 10, '101111': 26, '110001': 11, '011110': 1, '111011': 31, '110011': 6, '111101': 39, '100111': 5, '100010': 9, '000001': 7, '100100': 8, '100001': 9, '011111': 25, '001000': 3, '010111': 1, '101000': 6, '000010': 4, '000110': 1, '110110': 1, '111010': 6, '100011': 4, '000100': 5, '101100': 1}

Counts for classical register index 0 from PUB 2:
{'111111': 937675, '100000': 11558, '111110': 11023, '110000': 11404, '111100': 11246, '000000': 5385, '111000': 11408, '100001': 19, '011111': 32, '100010': 7, '110111': 38, '111101': 46, '111001': 8, '101111': 41, '110001': 7, '100100': 13, '111011': 25, '110011': 6, '001111': 1, '101000': 7, '111010': 10, '000010': 3, '100011': 6, '101010': 1, '110010': 4, '000100': 7, '100111': 2, '110100': 7, '010000': 4, '000001': 2, '000011': 1, '001000': 2, '110110': 1, '000111': 1}

Counts for classical register index 0 from PUB 3:
{'000000': 991604, '001000': 685, '000010': 727, '111111': 3426, '000111': 144, '100000': 683, '000001': 807, '000100': 714, '010000': 681, '010110': 1, '001111': 139, '011111': 121, '000011': 180, '011101': 1, '101000': 2, '000110': 8, '010001': 2, '111101': 4, '000101': 7, '110000': 3, '001100': 5, '111000': 8, '111110': 7, '011110': 5, '001110': 2, '011000': 4, '100010': 3, '100100': 2, '111100': 8, '111011': 6, '001010': 1, '011011': 3, '001011': 2, '011100': 1, '010100': 1, '110111': 2, '100111': 1}

Counts for classical register index 1 from PUB 4:
{'000000': 998941, '111110': 8, '000100': 95, '000001': 158, '001000': 110, '011111': 35, '000010': 102, '010000': 115, '111111': 108, '000111': 51, '011100': 3, '011110': 2, '100000': 92, '000011': 55, '001111': 38, '000110': 10, '001101': 6, '010111': 2, '110110': 2, '000101': 8, '101111': 6, '111101': 5, '111011': 3, '001100': 6, '110111': 3, '011101': 3, '101011': 6, '011001': 1, '011000': 4, '110010': 1, '101000': 1, '001110': 4, '111001': 1, '001010': 2, '001011': 3, '011011': 1, '010110': 2, '111000': 1, '101110': 1, '101100': 1, '110000': 1, '010101': 1, '010010': 1, '111100': 1}

Counts for classical register index 1 from PUB 5:
{'111100': 11335, '111111': 939592, '100000': 11674, '111110': 11589, '111000': 11151, '111101': 303, '110000': 11256, '000000': 1940, '101111': 182, '110011': 29, '011111': 141, '110111': 203, '111011': 257, '110001': 23, '100011': 9, '101000': 16, '110010': 16, '111010': 45, '110101': 24, '111001': 57, '110110': 22, '100010': 7, '100111': 19, '100001': 5, '101100': 7, '000111': 2, '101101': 8, '001011': 3, '101001': 5, '110100': 16, '011100': 1, '011000': 2, '010010': 1, '010100': 2, '000011': 2, '010110': 1, '010101': 1, '001000': 4, '001111': 7, '010000': 3, '010001': 1, '101110': 7, '001110': 4, '011011': 2, '011101': 1, '000110': 2, '101011': 6, '010111': 1, '100101': 4, '011110': 2, '101010': 5, '001100': 1, '100100': 2, '100110': 1, '011010': 1}

Counts for classical register index 1 from PUB 6:
{'111111': 938647, '111110': 11456, '111100': 10973, '110000': 11432, '100000': 11399, '111000': 11246, '111011': 429, '000000': 2077, '101111': 377, '011111': 399, '110011': 40, '110111': 442, '111101': 487, '100101': 6, '100010': 8, '100001': 7, '101000': 13, '100111': 31, '101110': 28, '111010': 79, '110010': 19, '100011': 12, '110110': 36, '110101': 32, '110100': 23, '011010': 3, '011011': 24, '101100': 12, '110001': 26, '011110': 21, '101011': 23, '000011': 3, '111001': 63, '001111': 18, '011101': 13, '101101': 20, '000101': 1, '010111': 15, '100110': 7, '000111': 5, '001011': 3, '100100': 7, '010110': 3, '001001': 2, '010101': 1, '101010': 1, '001110': 3, '101001': 7, '010011': 2, '000001': 1, '011001': 3, '011000': 3, '010010': 2, '011100': 4, '001000': 1, '001101': 1, '010000': 2, '010001': 1, '001100': 1}

Counts for classical register index 1 from PUB 7:
{'000000': 994271, '111111': 3211, '111110': 210, '111011': 195, '111100': 22, '111101': 193, '000011': 80, '010111': 18, '000100': 94, '111000': 12, '001000': 90, '110111': 226, '000001': 144, '110011': 14, '110110': 11, '100011': 3, '111010': 16, '011111': 320, '101111': 220, '000111': 69, '000010': 96, '100000': 102, '011110': 17, '110101': 10, '001111': 79, '010000': 106, '001101': 1, '101101': 9, '011011': 16, '100111': 21, '110000': 11, '101011': 13, '110100': 2, '001110': 12, '010110': 2, '100101': 3, '101100': 4, '001011': 6, '101110': 9, '011101': 10, '111001': 15, '000110': 3, '000101': 5, '101000': 3, '100110': 2, '011100': 4, '110010': 2, '010011': 2, '010100': 1, '010101': 3, '100001': 2, '001010': 1, '001001': 1, '001100': 4, '011000': 1, '110001': 1, '011001': 1, '011010': 1}

Real hardware run on ibm_kingston

which_qubit_to_analyze = 1  # change this index to analyze different qubits from the chosen_qubits list

print(f"Total number of chosen qubits: {len(chosen_qubits)}")
print(f"\nAnalyzing qubit index {which_qubit_to_analyze}, physical qubit {chosen_qubits[which_qubit_to_analyze]}:\n")

Total number of chosen qubits: 2

Analyzing qubit index 1, physical qubit 117:

counting_list = counting_lists[which_qubit_to_analyze] 

chosen_sol_index = 0  # Choose which solution to analyze (no real effect because 2 solutions belongs to same gauge family)
print(f"\nAnalyzing results for qubit index {chosen_qubits[which_qubit_to_analyze]}:\n") 

ibm_prob_list = counts_to_probabilities(counting_list)
# display(ibm_prob_list)

emp_inv = derived_constraints_from_empirical_probs(ibm_prob_list)
# print(emp_inv)

sols = reconstruct_instrument_from_invariants_mixed_det(
    trM0 = emp_inv[0],
    detM0 = emp_inv[1],
    trM1 = emp_inv[2],
    detM1 = emp_inv[3],
    S0 = emp_inv[6] * 2, # S0 = sum of all elements of P0 = 2 * probability of observing 0 from v0=(0.5,0.5)
    gauge_p00 = 1.0
)

print(f"\nthere are {len(sols)} solutions from empirical data,\nwe choose solution with index {chosen_sol_index}.\n")

reconstructed_MCM, t_bound = plot_gauge_transformation_effects(
    Instrument2x2(sols[chosen_sol_index][0], sols[chosen_sol_index][1]), 
    t_width_factor=1.5,
    verbose=False,
    MCM_reference=[],
    resolution=1e5,
    margin_tol=0.00028
    )

print("\nReconstructed Instrument from invariants (empirical data):")
if reconstructed_MCM is not None:
    reconstructed_MCM.reveal()

Analyzing results for qubit index 117:


there are 2 solutions from empirical data,
we choose solution with index 0.

================================================================================
Valid Intervals for Derived Quantities (Readout & Back-action Errors)
================================================================================

	Quantity	Min	Max	Width	Value
0	prep 0 meas 1	0.00020175	0.00048543	2.83679e-04	(3.436 ± 1.418)e-4
1	prep 1 meas 0	0.00017875	0.00046243	2.83679e-04	(3.206 ± 1.418)e-4
2	prep 0 excite to 1	0.00010132	0.00010460	3.27872e-06	(1.030 ± 0.016)e-4
3	prep 1 decay to 0	0.01165149	0.01165477	3.27872e-06	(1.165 ± 0.000)e-2

================================================================================
Valid Intervals for Matrix Entries
================================================================================

	Entry	Min	Max	Width	Value
0	M^0[0,0] = p_0^(0,0)	0.99979455	0.99979457	2.01313e-08	(9.998 ± 0.000)e-1
1	M^0[0,1] = p_1^(0,0)	-0.00028000	0.00000368	2.83679e-04	(-1.382 ± 1.418)e-4
2	M^0[1,0] = p_0^(0,1)	-0.00028000	0.00000368	2.83679e-04	(-1.382 ± 1.418)e-4
3	M^0[1,1] = p_1^(0,1)	0.00045875	0.00045877	2.01313e-08	(4.588 ± 0.000)e-4
4	M^1[0,0] = p_0^(1,0)	0.00010083	0.00010411	3.27872e-06	(1.025 ± 0.016)e-4
5	M^1[0,1] = p_1^(1,0)	0.01165109	0.01193149	2.80400e-04	(1.179 ± 0.014)e-2
6	M^1[1,0] = p_0^(1,1)	0.00010092	0.00038132	2.80400e-04	(2.411 ± 1.402)e-4
7	M^1[1,1] = p_1^(1,1)	0.98788648	0.98788976	3.27872e-06	(9.879 ± 0.000)e-1

================================================================================



/var/folders/dk/_dd6ng8n3yv_crq6vrdt8qzc0000gn/T/ipykernel_62248/724936655.py:540: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()

Reconstructed Instrument from invariants (empirical data):

MCM Instrument:

M0 matrix (outcome 0):
  [[p_0^(0,0), p_1^(0,0)]] = [[0.99979456, -0.00013816]]
  [[p_0^(0,1), p_1^(0,1)]]   [[-0.00013816, 0.00045876]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00010247, 0.01179129]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00024112, 0.98788812]]

1. Obtain the framework consistency/integrity (%) -> done.

2. use learned MCM to bound SP error

3. apply QEM using the learned SP (includes initial SP and back action of MCM) and M error models, then compare with IBM state-of-the-art QEM performance

4. Applications:

5. Superdense Coding: https://quantum.cloud.ibm.com/learning/en/courses/basics-of-quantum-information/entanglement-in-action/qiskit-implementation#superdense-coding

6. Long-range entanglement with dynamic circuits: https://quantum.cloud.ibm.com/docs/en/tutorials/long-range-entanglement?utm_source=chatgpt.com

---

References:

• Blog post of MCM availability on Nov 19: https://www.ibm.com/quantum/blog/utility-scale-dynamic-circuits

Additional Tasks

1. Obtain the framework consistency/integrity (%)

• Actually we can use the MCM reconstructed from len 3 probs -> to verify the len 4 and larger ones' probs -> get the framework consistency (%)

• RMSE for prob strings of length =3 should be same magnitude as shot noise, and RMSE for prob strings of length >3 should be slightly increasing but still close to shot noise level.

• RMSE will also be affected by leakage errors and single-qubit gate errors.

# Use the reconstructed MCM to compare with IBM data
if reconstructed_MCM is not None:
    theoretical_probs, rmse, df_comparison = compare_instrument_to_empirical(
        instrument=reconstructed_MCM,
        empirical_probs=ibm_prob_list,
        v0=np.array([0.5, 0.5]),
        max_len=max_word_length
    )

    print(f"Overall RMSE between reconstructed MCM and IBM data: {rmse:.6e}")
    print(f"\nNumber of strings compared: {len(df_comparison)}")
    print("\nDetailed comparison:")
    # sort df_comparison's row by first using String's length ascending, then by string lexicographically
    df_comparison = df_comparison.sort_values(by=['Length', 'String']).reset_index(drop=True)

    display(df_comparison)

    # Summary statistics by string length
    print("\n" + "="*80)
    print("RMSE by string length:")
    print("="*80)
    for length in range(1, max_word_length + 1):
        mask = df_comparison['Length'] == length
        if mask.any():
            length_diffs = [float(df_comparison[mask]['Difference'].iloc[i].replace('+', ''))
                            for i in range(mask.sum())]
            length_rmse = np.sqrt(np.mean([d**2 for d in length_diffs]))
            print(f"Length {length}: RMSE = {length_rmse:.6e} ({mask.sum()} strings)")

Overall RMSE between reconstructed MCM and IBM data: 2.199958e-05

Number of strings compared: 126

Detailed comparison:

	String	Length	Theoretical	Empirical (ibm)	Difference	Abs Diff	Rel Error
0	0	1	0.49998850	0.49998850	+5.551115e-17	5.551115e-17	0.0000%
1	1	1	0.50001150	0.50001150	+0.000000e+00	0.000000e+00	0.0000%
2	00	2	0.49965651	0.49965650	+1.005945e-08	1.005945e-08	0.0000%
3	01	2	0.00033198	0.00033200	-1.508894e-08	1.508894e-08	0.0045%
4	10	2	0.00610324	0.00610325	-1.474015e-08	1.474015e-08	0.0002%
...	...	...	...	...	...	...	...
121	111011	6	0.00021904	0.00022100	-1.960858e-06	1.960858e-06	0.8873%
122	111100	6	0.00561718	0.00558275	+3.442735e-05	3.442735e-05	0.6167%
123	111101	6	0.00021960	0.00024700	-2.740168e-05	2.740168e-05	11.0938%
124	111110	6	0.00576610	0.00581575	-4.965179e-05	4.965179e-05	0.8537%
125	111111	6	0.47041635	0.47038950	+2.684994e-05	2.684994e-05	0.0057%

126 rows × 7 columns

================================================================================
RMSE by string length:
================================================================================
Length 1: RMSE = 3.925231e-17 (2 strings)
Length 2: RMSE = 1.526758e-08 (4 strings)
Length 3: RMSE = 1.370968e-05 (8 strings)
Length 4: RMSE = 2.284385e-05 (16 strings)
Length 5: RMSE = 2.358547e-05 (32 strings)
Length 6: RMSE = 2.281988e-05 (64 strings)

2. Obtain SP error bounds from learned MCM

# pauli_labels must match the order used in create_list_of_circuits
PAULI_LABELS = ("I", "X", "Y", "Z")


def extract_counts_by_qubit_and_pauli(
    job_result,
    num_qubits: int,
    pauli_labels: Tuple[str, ...] = PAULI_LABELS,
) -> Dict[Tuple[int, str], Dict[str, int]]:
    """
    From a SamplerV2 PrimitiveResult, reconstruct:
        (qubit_index, Pauli_label) -> counts over repeated-MCM bitstrings.

    Assumes:
        - Circuits were created with create_list_of_circuits(max_word_length, num_qubits)
        - Circuits are ordered as:
            for q in range(num_qubits):
                for gate_name in pauli_labels:
                    ...
        - For each circuit, the classical register used for qubit q is the
          q-th entry in pub_result.data.keys() (as in your plot_job_results).

    Returns
    -------
    result_counts : dict[(int,str), dict[str,int]]
        Example key: (0, "X") -> {'000000': ..., '111000': ..., ...}
    """
    num_pubs = len(job_result)
    num_paulis = len(pauli_labels)

    if num_pubs != num_qubits * num_paulis:
        raise ValueError(
            f"Expected {num_qubits * num_paulis} pubs, got {num_pubs}. "
            "Check pauli_labels and circuit ordering."
        )

    result_counts: Dict[Tuple[int, str], Dict[str, int]] = {}

    for i, pub in enumerate(job_result):
        q = i // num_paulis          # which qubit this circuit targets
        p_idx = i % num_paulis       # which Pauli in PAULI_LABELS
        P = pauli_labels[p_idx]

        # pub.data keys are in the order of classical registers (c0, c1, ...)
        keys = list(pub.data.keys())
        if q >= len(keys):
            raise RuntimeError(
                f"For PUB {i}, expected at least {q+1} classical subsets, got {len(keys)}"
            )
        creg_key = keys[q]
        counts = pub.data[creg_key].get_counts()

        result_counts[(q, P)] = dict(counts)

    return result_counts

def get_prob_meas_0_from_counts(
        result_counts: Dict[Tuple[int, str], Dict[str, int]],
        qubit_index: int,
        filter_paulis: Tuple[str, ...] = ("I", "Z"),
        ) -> float:
    """
    Extracts the probability of measuring '0' for a specific qubit from the result counts.
    It aggregates counts from 'I' and 'Z' Pauli basis preparations.
    """
    # filter the result_counts to only show the specified qubit's I and Z data
    counts_list = []
    for key in result_counts:
        if key[0] == qubit_index and key[1] in filter_paulis:
            # print(f"\nCounts for qubit {key[0]} with Pauli {key[1]}:")
            # print(result_counts[key])
            counts_list.append(result_counts[key])

    # convert the filtered counting dict into prob dict by using `counts_to_probabilities()`
    filtered_ibm_prob_list = counts_to_probabilities(counts_list)

    # print(filtered_ibm_prob_list)
    prob_0 = filtered_ibm_prob_list.get('0', 0.0)

    print(f"Qubit {qubit_index}: prob of reading out 0 (noisy state prep): {prob_0:.5e}")
    print(f"Qubit {qubit_index}: prob of reading out 1 (noisy state prep): {filtered_ibm_prob_list.get('1', 0.0):.5e}")

    return prob_0

result_counts = extract_counts_by_qubit_and_pauli(job_result[:-num_qubits*1], num_qubits=num_qubits) # have to take subset of job_result if needed to exclude accuracy verification circuits

prob_meas_0_val = get_prob_meas_0_from_counts(result_counts, qubit_index=which_qubit_to_analyze) 

print(f"\nAnalyzing results for qubit index {chosen_qubits[which_qubit_to_analyze]}:\n") 
print(t_bound, prob_meas_0_val)

Qubit 1: prob of reading out 0 (noisy state prep): 9.97608e-01
Qubit 1: prob of reading out 1 (noisy state prep): 2.39150e-03

Analyzing results for qubit index 117:

(np.float64(-0.00014195384834125857), np.float64(0.00014191355798992555)) 0.9976084999999999

refined_t_bound, epsilon_bound = refine_bounds_with_sp_error(
    reconstructed_MCM,
    t_bound,
    prob_meas_0_val,
    margin_tol=0.00,
    verbose=True
)

print(f"\nRefined t bound: {refined_t_bound}")
print(f"😍 Estimated SP error bound for noisy |0⟩ preparation (epsilon):\n{epsilon_bound[0]:.4e}, {epsilon_bound[1]:.4e}") 

--- Refine Bounds Debug ---
P(outcome=1): 0.002392
Baseline P(1|0): 0.000344, P(1|1): 0.999679
Input t_bound: [-0.000142, 0.000142]
Epsilon range on input t: [0.001908, 0.002191]
Target epsilon range: [0.000000, 1.000000]
Refined t_bound: [-0.000142, 0.000142]
Refined epsilon: [0.001908, 0.002191]
---------------------------

Refined t bound: (np.float64(-0.00014195384834125835), np.float64(0.00014191355798992545))
😍 Estimated SP error bound for noisy |0⟩ preparation (epsilon):
1.9079e-03, 2.1906e-03

def calculate_and_display_refined_intervals(
    reconstructed_MCM: Instrument2x2 | Any,
    refined_t_bound: Tuple[float, float],
    resolution: int = 1000
) -> Tuple[pd.DataFrame, pd.DataFrame]:
    """
    Calculates and displays the min/max intervals for matrix entries and derived quantities
    of a reconstructed MCM instrument over a refined gauge parameter range.

    Args:
        reconstructed_MCM: The Instrument2x2 object at the center of the gauge (t=0).
        refined_t_bound: The valid interval (t_min, t_max) for the gauge parameter.
        resolution: Number of points to sample within the t_bound for min/max calculation.

    Returns:
        Tuple containing:
            - DataFrame for derived quantities intervals.
            - DataFrame for matrix entries intervals.
    """
    t_min, t_max = refined_t_bound

    # Handle the case where t_min == t_max (single point)
    if abs(t_max - t_min) < 1e-15:
        t_values = np.array([t_min])
    else:
        t_values = np.linspace(t_min, t_max, resolution)

    # Store the 8 entries of the transformed instrument for each value of t
    transformed_entries = []

    for t_val in t_values:
        M0_prime, M1_prime = gauge_transform_instrument_numerically(
            reconstructed_MCM.M0, reconstructed_MCM.M1, t_val
        )
        entries = np.concatenate((M0_prime.flatten(), M1_prime.flatten()))
        transformed_entries.append(entries)

    transformed_entries = np.array(transformed_entries)

    # Helper for formatting value column
    def fmt_val_err(min_v, max_v):
        c = (min_v + max_v) / 2
        h = (max_v - min_v) / 2

        ref = abs(c) if abs(c) > 0 else abs(h)
        if ref == 0:
            return "(0.000 ± 0.000)e+0"

        exponent = int(np.floor(np.log10(ref)))
        scale = 10.0 ** (-exponent)

        c_s = c * scale
        h_s = h * scale

        return f"({c_s:.3f} ± {h_s:.3f})e{exponent:+d}"

    # ===== 1. Matrix Entries =====
    entry_labels = [
        "M^0[0,0] = p_0^(0,0)",
        "M^0[0,1] = p_1^(0,0)",
        "M^0[1,0] = p_0^(0,1)",
        "M^0[1,1] = p_1^(0,1)",
        "M^1[0,0] = p_0^(1,0)",
        "M^1[0,1] = p_1^(1,0)",
        "M^1[1,0] = p_0^(1,1)",
        "M^1[1,1] = p_1^(1,1)"
    ]

    entry_ranges_data = []

    for entry_idx, entry_label in enumerate(entry_labels):
        block_values = transformed_entries[:, entry_idx]
        min_val = np.min(block_values)
        max_val = np.max(block_values)
        width = max_val - min_val

        entry_ranges_data.append({
            'Entry': entry_label,
            'Min': f"{min_val:.8f}",
            'Max': f"{max_val:.8f}",
            'Width': f"{width:.5e}",
            'Value': fmt_val_err(min_val, max_val)
        })

    # ===== 2. Derived Quantities =====
    prep0_meas1 = transformed_entries[:, 4] + transformed_entries[:, 6]
    prep1_meas0 = transformed_entries[:, 1] + transformed_entries[:, 3]
    prep0_excite = transformed_entries[:, 2] + transformed_entries[:, 6]
    prep1_decay = transformed_entries[:, 1] + transformed_entries[:, 5]

    derived_quantities = {
        "prep 0 meas 1": prep0_meas1,
        "prep 1 meas 0": prep1_meas0,
        "prep 0 excite to 1": prep0_excite,
        "prep 1 decay to 0": prep1_decay
    }

    quantity_ranges_data = []
    for quantity_name, block_values in derived_quantities.items():
        min_val = np.min(block_values)
        max_val = np.max(block_values)
        width = max_val - min_val

        quantity_ranges_data.append({
            'Quantity': quantity_name,
            'Min': f"{min_val:.8f}",
            'Max': f"{max_val:.8f}",
            'Width': f"{width:.5e}",
            'Value': fmt_val_err(min_val, max_val)
        })

    # ===== Display DataFrames =====
    print("="*80)
    print(f"Refined Intervals over t ∈ [{t_min:.6e}, {t_max:.6e}]")
    print("="*80)

    print("\nDerived Quantities (Readout & Back-action Errors):")
    df_quantities = pd.DataFrame(quantity_ranges_data)
    display(df_quantities)

    print("\nMatrix Entries:")
    df_entries = pd.DataFrame(entry_ranges_data)
    display(df_entries)
    print("="*80 + "\n")

    return df_quantities, df_entries

bounded_quantities, bounded_entries = calculate_and_display_refined_intervals(reconstructed_MCM, refined_t_bound)

================================================================================
Refined Intervals over t ∈ [-1.419538e-04, 1.419136e-04]
================================================================================

Derived Quantities (Readout & Back-action Errors):

	Quantity	Min	Max	Width	Value
0	prep 0 meas 1	0.00020173	0.00048541	2.83679e-04	(3.436 ± 1.418)e-4
1	prep 1 meas 0	0.00017874	0.00046242	2.83679e-04	(3.206 ± 1.418)e-4
2	prep 0 excite to 1	0.00010132	0.00010459	3.27871e-06	(1.030 ± 0.016)e-4
3	prep 1 decay to 0	0.01165149	0.01165477	3.27871e-06	(1.165 ± 0.000)e-2

Matrix Entries:

	Entry	Min	Max	Width	Value
0	M^0[0,0] = p_0^(0,0)	0.99979456	0.99979458	2.01318e-08	(9.998 ± 0.000)e-1
1	M^0[0,1] = p_1^(0,0)	-0.00028000	0.00000368	2.83679e-04	(-1.382 ± 1.418)e-4
2	M^0[1,0] = p_0^(0,1)	-0.00028000	0.00000368	2.83679e-04	(-1.382 ± 1.418)e-4
3	M^0[1,1] = p_1^(0,1)	0.00045874	0.00045876	2.01318e-08	(4.588 ± 0.000)e-4
4	M^1[0,0] = p_0^(1,0)	0.00010081	0.00010409	3.27872e-06	(1.025 ± 0.016)e-4
5	M^1[0,1] = p_1^(1,0)	0.01165109	0.01193149	2.80400e-04	(1.179 ± 0.014)e-2
6	M^1[1,0] = p_0^(1,1)	0.00010092	0.00038132	2.80400e-04	(2.411 ± 1.402)e-4
7	M^1[1,1] = p_1^(1,1)	0.98788650	0.98788978	3.27872e-06	(9.879 ± 0.000)e-1

================================================================================

3. Learn Pauli-X error rate from learned MCM

4. Accuracy test

accuracy_verification_lists = plot_job_results(job_result[-num_qubits*1:], backend_name=backend.name, num_qubits=num_qubits)
# accuracy_verification_lists

Counts for classical register index 0 from PUB 0:
{'101010': 939744, '010101': 4897, '110101': 12665, '101011': 14431, '101101': 13262, '101001': 5500, '100101': 5462, '111111': 706, '110010': 93, '111010': 855, '100100': 30, '101110': 872, '011010': 81, '110111': 19, '101111': 586, '110100': 85, '111101': 35, '010100': 39, '101000': 42, '010010': 23, '110110': 162, '010110': 68, '001010': 45, '100010': 67, '101100': 82, '111011': 15, '010111': 11, '011111': 6, '100110': 74, '110001': 1, '111110': 14, '011101': 5, '111001': 6, '001111': 1, '100111': 6, '111100': 2, '011011': 2, '001011': 1, '000101': 2, '100011': 1, '110011': 2}

Counts for classical register index 1 from PUB 1:
{'101010': 960684, '101101': 9947, '101011': 10772, '110101': 9861, '100101': 2042, '101000': 178, '010101': 1859, '101111': 571, '101001': 1989, '101110': 308, '111111': 462, '110111': 54, '100010': 143, '111010': 265, '001010': 145, '011111': 23, '111101': 104, '111110': 64, '110110': 129, '010111': 5, '111011': 51, '110010': 31, '100100': 10, '111000': 5, '111100': 9, '110100': 28, '100110': 28, '011010': 13, '100111': 35, '010110': 22, '110000': 1, '100001': 10, '010010': 5, '010100': 10, '011011': 4, '000101': 4, '101100': 50, '110001': 7, '001101': 4, '011101': 7, '111001': 13, '110011': 10, '100011': 11, '010011': 3, '001001': 1, '010000': 2, '001011': 2, '100000': 2, '001111': 4, '000010': 1, '011110': 4, '001000': 2, '011001': 2, '001110': 3, '000111': 1}

Real hardware run on ibm_kingston

ibm_verify_prob_list = counts_to_probabilities(accuracy_verification_lists[which_qubit_to_analyze]) 

def predict_and_compare_mcm_readouts(
    epsilon_bound: Tuple[float, float],
    bounded_quantities: pd.DataFrame,
    bounded_entries: pd.DataFrame,
    measure_2_error: float,
    pauli_x_error: float,
    ibm_verify_prob_list: Dict[str, float],
    max_word_length: int = 6
):
    """
    Predicts the probability of measuring '0' at each step of a repeated X-MCM circuit
    using three models (Learned, IBM Simplified, Ideal) and compares them with
    empirical QPU data.

    The circuit consists of repeated cycles of: X gate -> Mid-Circuit Measurement.

    Args:
        epsilon_bound: Tuple (min, max) for the state preparation error epsilon.
        bounded_quantities: DataFrame containing min/max values for MCM properties
                          (prep 0 meas 1, prep 1 meas 0, prep 0 excite to 1, prep 1 decay to 0).
        bounded_entries: DataFrame containing min/max values for MCM entries (unused).
        measure_2_error: The single scalar error rate provided by IBM for MCM readout.
        pauli_x_error: The single scalar error rate for X gate (depolarizing-like).
        ibm_verify_prob_list: Dictionary of empirical probabilities from the QPU.
        max_word_length: Number of MCMs in the sequence.
    """

    # --- Helper to extract interval from DataFrame ---
    def get_interval(name):
        row = bounded_quantities[bounded_quantities['Quantity'] == name].iloc[0]
        return float(row['Min']), float(row['Max'])

    # --- 1. Extract Learned Model Parameters (Intervals) ---
    # Readout errors
    p0_m1_int = get_interval("prep 0 meas 1")
    p1_m0_int = get_interval("prep 1 meas 0")

    # Back-action (state update) probabilities
    p0_excite_int = get_interval("prep 0 excite to 1")
    p1_decay_int = get_interval("prep 1 decay to 0")

    # State prep error
    eps_int = epsilon_bound

    # --- 2. Initialize States for all 3 models ---
    # We track p(1) (probability of being in state |1>)

    # A. Learned Model (Interval arithmetic)
    # Initial state |0> noisy: v0 = [1-eps, eps] -> p(1) is eps
    current_p1_interval = eps_int

    # B. IBM Simplified Model (Scalar)
    # Initial state |0> perfect: p(1) = 0
    current_p1_ibm = 0.0

    # C. Ideal Model (Scalar)
    # Initial state |0> perfect: p(1) = 0
    current_p1_ideal = 0.0

    # Storage for predictions
    predictions = []

    # --- 3. Simulate the Circuit Step-by-Step ---
    # The circuit is: Repeat(X -> Measure)

    for step in range(1, max_word_length + 1):
        # --- Step 3a: Apply X Gate ---
        # Rule: v_after = [(1 - px) X + px I] v_before
        # p1_new = (1 - px)(1 - p1) + px * p1 = 1 - p1 - px + 2*px*p1 = (1 - px) + p1 * (2*px - 1)

        # Learned (Interval)
        # Since (2*px - 1) is negative (assuming error < 0.5), max maps to min and min maps to max.
        px = pauli_x_error
        p1_min, p1_max = current_p1_interval

        next_p1_min = (1 - px) + p1_max * (2 * px - 1)
        next_p1_max = (1 - px) + p1_min * (2 * px - 1)
        current_p1_interval = (next_p1_min, next_p1_max)

        # IBM (Scalar)
        current_p1_ibm = (1 - px) + current_p1_ibm * (2 * px - 1)

        # Ideal (Scalar) - Perfect X
        current_p1_ideal = 1.0 - current_p1_ideal

        # --- Step 3b: Calculate Readout Probability for '0' ---
        # P(r=0) = (1 - p0m1)(1 - p1) + p1m0 * p1
        #        = (1 - p0m1) + p1 * (p1m0 + p0m1 - 1)
        # Coefficient (p1m0 + p0m1 - 1) is negative.

        # Learned Model Prediction (Interval)
        # Maximize P(r=0): use min p1, min p0m1, max p1m0 (since partials: dp1<0, dp0m1<0, dp1m0>0)
        # Minimize P(r=0): use max p1, max p0m1, min p1m0

        p1_min, p1_max = current_p1_interval
        e01_min, e01_max = p0_m1_int
        e10_min, e10_max = p1_m0_int

        # Max P(r=0)
        pred_learned_max = (1 - e01_min) * (1 - p1_min) + e10_max * p1_min
        # Min P(r=0)
        pred_learned_min = (1 - e01_max) * (1 - p1_max) + e10_min * p1_max

        pred_learned = (pred_learned_min, pred_learned_max)

        # IBM Model Prediction
        # P(r=0) = (1 - m2_err)(1 - p1) + m2_err * p1
        pred_ibm = (1 - measure_2_error) * (1 - current_p1_ibm) + measure_2_error * current_p1_ibm

        # Ideal Model Prediction
        pred_ideal = 1.0 - current_p1_ideal

        # --- Step 3c: Apply Back-action (State Update) ---
        # p1_new = u01(1 - p1) + (1 - u10)p1
        #        = u01 + p1 * (1 - u10 - u01)
        # Coefficient (1 - u10 - u01) is positive.

        # Learned Model Update
        # Max p1_new: max p1, max u01, min u10
        # Min p1_new: min p1, min u01, max u10

        u01_min, u01_max = p0_excite_int
        u10_min, u10_max = p1_decay_int

        next_state_p1_max = u01_max * (1 - p1_max) + (1 - u10_min) * p1_max
        next_state_p1_min = u01_min * (1 - p1_min) + (1 - u10_max) * p1_min

        current_p1_interval = (next_state_p1_min, next_state_p1_max)

        # IBM Model Update (No back-action modeled -> Identity)
        # current_p1_ibm remains unchanged

        # Ideal Model Update (Identity)
        # current_p1_ideal remains unchanged

        # --- Step 3d: Extract Empirical Data ---
        # Sum probabilities of all strings of length `step` ending in '0'
        emp_prob_0 = 0.0
        found_keys = False
        for s, prob in ibm_verify_prob_list.items():
            if len(s) == step and s.endswith('0'):
                emp_prob_0 += prob
                found_keys = True

        if not found_keys and step == 1:
             emp_prob_0 = ibm_verify_prob_list.get('0', 0.0)

        predictions.append({
            "MCM Index": step,
            "Ideal": pred_ideal,
            "IBM Simple": pred_ibm,
            "Learned Min": pred_learned[0],
            "Learned Max": pred_learned[1],
            "Empirical": emp_prob_0
        })

    # --- 4. Visualization ---
    df_res = pd.DataFrame(predictions)

    # Calculate differences relative to Empirical
    df_res["Diff Ideal"] = df_res["Ideal"] - df_res["Empirical"]
    df_res["Diff IBM"] = df_res["IBM Simple"] - df_res["Empirical"]
    df_res["Diff Learned Min"] = df_res["Learned Min"] - df_res["Empirical"]
    df_res["Diff Learned Max"] = df_res["Learned Max"] - df_res["Empirical"]

    # --- Plot 1: Absolute Probabilities ---
    plt.figure(figsize=(12, 6))

    # Plot Learned Interval (Shaded Region)
    plt.fill_between(df_res["MCM Index"], df_res["Learned Min"], df_res["Learned Max"], 
                     color='green', alpha=0.3, label='Learned Model (Interval)')

    # Plot IBM Simple Model
    plt.plot(df_res["MCM Index"], df_res["IBM Simple"], 'b--o', label=f'IBM Simple Model (err={measure_2_error:.1e})')

    # Plot Ideal Model
    plt.plot(df_res["MCM Index"], df_res["Ideal"], 'k:', alpha=0.5, label='Ideal')

    # Plot Empirical Data
    plt.plot(df_res["MCM Index"], df_res["Empirical"], 'r-x', linewidth=1.5, label='Empirical Data')

    plt.xlabel('MCM Index (Step in Circuit)')
    plt.ylabel('Probability of Measuring 0')
    plt.title(f'Repeated X-MCM: Absolute Probabilities (P(0))\n(State Prep Error $\epsilon \in [{epsilon_bound[0]:.1e}, {epsilon_bound[1]:.1e}]$)')
    plt.xticks(df_res["MCM Index"])
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.show()

    # --- Plot 2: Deviations from Empirical ---
    plt.figure(figsize=(12, 6))

    # Plot Learned Interval (Shaded Region) relative to Empirical (which is at 0)
    plt.fill_between(df_res["MCM Index"], df_res["Diff Learned Min"], df_res["Diff Learned Max"], 
                     color='green', alpha=0.3, label='Learned Model Deviation (Interval)')

    # Plot IBM Simple Model Deviation
    plt.plot(df_res["MCM Index"], df_res["Diff IBM"], 'b--o', label=f'IBM Simple Model Deviation (err={measure_2_error:.1e})')

    # Plot Ideal Model Deviation
    plt.plot(df_res["MCM Index"], df_res["Diff Ideal"], 'k:', alpha=0.5, label='Ideal Deviation')

    # Plot Empirical Data Baseline (at y=0)
    plt.axhline(y=0, color='r', linestyle='-', linewidth=1.5, label='Empirical Baseline (0)')

    plt.xlabel('MCM Index (Step in Circuit)')
    plt.ylabel('Difference from Empirical Probability (P(0))')
    plt.title(f'Repeated X-MCM: Model Deviations from Empirical Data\n(State Prep Error $\epsilon \in [{epsilon_bound[0]:.1e}, {epsilon_bound[1]:.1e}]$)')
    plt.xticks(df_res["MCM Index"])
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.show()

    # Display Table
    print("\nDetailed Predictions vs Empirical:")
    # Format columns for display
    df_disp = df_res.copy()

    # Calculate deviations for table
    df_disp['Learned Mid'] = (df_disp['Learned Min'] + df_disp['Learned Max']) / 2
    df_disp['Dev IBM'] = df_disp['IBM Simple'] - df_disp['Empirical']
    df_disp['Dev Learned Mid'] = df_disp['Learned Mid'] - df_disp['Empirical']

    # Calculate ratio of absolute errors (IBM Error / Learned Error)
    # A value > 1 means IBM error is larger (Learned is better)
    df_disp['Error Ratio (IBM/Learned)'] = df_disp['Dev IBM'].abs() / df_disp['Dev Learned Mid'].abs()

    df_disp['Learned Interval'] = df_disp.apply(lambda x: f"[{x['Learned Min']:.4f}, {x['Learned Max']:.4f}]", axis=1)
    df_disp['In Interval?'] = df_disp.apply(lambda x: "✅" if x['Learned Min'] <= x['Empirical'] <= x['Learned Max'] else "❌", axis=1)

    cols = ['MCM Index', 'Ideal', 'IBM Simple', 'Learned Interval', 'Empirical', 'In Interval?', 'Dev IBM', 'Dev Learned Mid', 'Error Ratio (IBM/Learned)']
    display(df_disp[cols])

    return df_res

predict_and_compare_mcm_readouts(
    epsilon_bound,
    bounded_quantities,
    bounded_entries,
    pauli_x_error=1.617e-4, # look up from IBM Compute resources page for the backend used
    measure_2_error=3.418e-3, # look up from IBM Compute resources page for the backend used
    ibm_verify_prob_list=ibm_verify_prob_list,
    max_word_length=max_word_length
)

Detailed Predictions vs Empirical:

	MCM Index	Ideal	IBM Simple	Learned Interval	Empirical	In Interval?	Dev IBM	Dev Learned Mid	Error Ratio (IBM/Learned)
0	1	0.0	0.003579	[0.0022, 0.0028]	0.002126	❌	0.001453	0.000403	3.600879
1	2	1.0	0.996261	[0.9854, 0.9860]	0.986947	❌	0.009314	-0.001276	7.299643
2	3	0.0	0.003900	[0.0141, 0.0147]	0.014314	✅	-0.010414	0.000088	119.011577
3	4	1.0	0.995940	[0.9737, 0.9742]	0.974352	❌	0.021588	-0.000410	52.694949
4	5	0.0	0.004220	[0.0257, 0.0263]	0.026147	✅	-0.021927	-0.000158	138.440903
5	6	1.0	0.995619	[0.9622, 0.9628]	0.962142	❌	0.033477	0.000353	94.800299

	MCM Index	Ideal	IBM Simple	Learned Min	Learned Max	Empirical	Diff Ideal	Diff IBM	Diff Learned Min	Diff Learned Max
0	1	0.0	0.003579	0.002246	0.002812	0.002126	-0.002126	0.001453	0.000120	0.000686
1	2	1.0	0.996261	0.985388	0.985954	0.986947	0.013053	0.009314	-0.001559	-0.000993
2	3	0.0	0.003900	0.014119	0.014684	0.014314	-0.014314	-0.010414	-0.000195	0.000370
3	4	1.0	0.995940	0.973659	0.974225	0.974352	0.025648	0.021588	-0.000693	-0.000127
4	5	0.0	0.004220	0.025706	0.026271	0.026147	-0.026147	-0.021927	-0.000441	0.000124
5	6	1.0	0.995619	0.962212	0.962778	0.962142	0.037858	0.033477	0.000070	0.000636

see all functions defined here

import inspect, textwrap

def user_functions(namespace=None):
    ns = globals() if namespace is None else namespace
    funcs = []
    for name, obj in ns.items():
        if inspect.isfunction(obj):
            # keep only functions defined in this notebook/kernel
            mod = getattr(obj, "__module__", None)
            if mod in (None, "__main__"):
                funcs.append((name, obj))
    return sorted(funcs, key=lambda x: x[0].lower())

def show_functions_with_docs(namespace=None, width=88):
    for name, fn in user_functions(namespace):
        sig = str(inspect.signature(fn))
        doc = inspect.getdoc(fn) or "(no docstring)"
        doc_wrapped = textwrap.fill(doc, width=width)
        print(f"🟢{name}{sig}\n{doc_wrapped}\n" + "-"*width)

show_functions_with_docs()

🟢_abgd_from_M(M)
(no docstring)
----------------------------------------------------------------------------------------
🟢_D_intervals_for_M(M, D_positive=True, tol=1e-24, margin_tol=0.0)
Return feasible D-intervals for a single M. For D>0: enforce -margin_tol <= a',b',c',d'
<= 1+margin_tol. For D<0: same but inequality directions flip (because 2D < 0).
----------------------------------------------------------------------------------------
🟢_gauge_transform_pair(M0, M1, t)
Return (M0', M1') = (R^-1 M0 R, R^-1 M1 R).
----------------------------------------------------------------------------------------
🟢_intersect_interval_lists(A, B, tol=1e-24)
(no docstring)
----------------------------------------------------------------------------------------
🟢_intersect_many(list_of_interval_lists, tol=1e-24)
(no docstring)
----------------------------------------------------------------------------------------
🟢_intervals_quad(a, b, c, rel, tol=1e-24)
Solve { x : a x^2 + b x + c (rel) 0 }, rel in {'ge','le'}. Returns a list of (lo, hi)
with lo/hi allowed to be ±inf.
----------------------------------------------------------------------------------------
🟢_map_D_to_t(D_intervals, tol=1e-24)
t = (1 - D)/2. Monotone decreasing map. Map each [D_lo, D_hi] to [t_lo, t_hi] with
t_lo=(1-D_hi)/2, t_hi=(1-D_lo)/2. Excludes t=1/2 automatically because D-intervals never
include D=0 (domain split).
----------------------------------------------------------------------------------------
🟢_merge_intervals(intervals, tol=1e-24)
(no docstring)
----------------------------------------------------------------------------------------
🟢_R(t)
R(t) = [[1-t, t],[t, 1-t]], valid for t != 1/2.
----------------------------------------------------------------------------------------
🟢_R_inv(t)
Closed-form inverse of R(t).
----------------------------------------------------------------------------------------
🟢allowed_t_regions_for_list(M_list, tol=1e-24, margin_tol=0.0)
Intersect allowed regions across several matrices.
----------------------------------------------------------------------------------------
🟢allowed_t_regions_for_M(M, tol=1e-24, margin_tol=0.0)
Input: M (2x2 numpy array, real). Output: list of (t_min, t_max) intervals such that
R(t)^{-1} @ M @ R(t) has all entries in [-margin_tol, 1+margin_tol].         Endpoints
may be ±np.inf. t=1/2 is excluded by construction.
----------------------------------------------------------------------------------------
🟢calculate_and_display_refined_intervals(reconstructed_MCM: __main__.Instrument2x2, refined_t_bound: Tuple[float, float], resolution: int = 1000) -> Tuple[pandas.core.frame.DataFrame, pandas.core.frame.DataFrame]
Calculates and displays the min/max intervals for matrix entries and derived quantities
of a reconstructed MCM instrument over a refined gauge parameter range.  Args:
reconstructed_MCM: The Instrument2x2 object at the center of the gauge (t=0).
refined_t_bound: The valid interval (t_min, t_max) for the gauge parameter.
resolution: Number of points to sample within the t_bound for min/max calculation.
Returns:     Tuple containing:         - DataFrame for derived quantities intervals.
- DataFrame for matrix entries intervals.
----------------------------------------------------------------------------------------
🟢calculate_average_probs_np(empirical_probs_list) -> Dict[str, float]
Calculates the average probabilities for each key across a list of dictionaries.  Args:
empirical_probs_list (list[dict]):          A list of dictionaries, where each
dictionary maps string keys to         probability floats.  Returns:     dict:
A single dictionary mapping each string key to its         average probability across
all input dictionaries.
----------------------------------------------------------------------------------------
🟢calculate_exact_all_string_probabilities_from_v0_and_instrument(inst: __main__.Instrument2x2, v0: numpy.ndarray, max_len: int) -> Dict[str, float]
Calculate exact probabilities for all binary strings up to max_len (L) given an initial
state distribution v0 and a MCM instrument. Exact Prob(w) = [1,1]^T M^{(w_L)} ...
M^{(w_1)} v0 for all |w|<=L. Args:     inst: Instrument2x2 instance representing the MCM
instrument.     v0: Initial state distribution as a numpy array of shape (2,).
max_len: Maximum length of binary strings to consider. Returns:     Dictionary mapping
binary strings to their exact probabilities.
----------------------------------------------------------------------------------------
🟢compare_instrument_to_empirical(instrument: __main__.Instrument2x2, empirical_probs: Dict[str, float], v0: numpy.ndarray = array([0.5, 0.5]), max_len: int = 4) -> Tuple[Dict[str, float], float, pandas.core.frame.DataFrame]
Calculate theoretical probabilities from an instrument and compare to empirical data.
Args:     instrument: Instrument2x2 object to evaluate     empirical_probs: Dictionary
of empirical probabilities from hardware     v0: Initial state distribution (default:
maximally mixed [0.5, 0.5])     max_len: Maximum length of binary strings to consider
Returns:     Tuple of:         - Dictionary of theoretical probabilities         -
Overall RMSE between theoretical and empirical         - DataFrame with detailed
comparison
----------------------------------------------------------------------------------------
🟢construct_instrument_symbolically(a: sympy.core.symbol.Symbol, b: sympy.core.symbol.Symbol, c: sympy.core.symbol.Symbol, e: sympy.core.symbol.Symbol, f: sympy.core.symbol.Symbol, g: sympy.core.symbol.Symbol) -> Tuple[sympy.matrices.dense.MutableDenseMatrix, sympy.matrices.dense.MutableDenseMatrix]
Constructs symbolic instrument matrices M0 and M1 from 6 independent parameters. The
parameters are probabilities for each column, which sum to 1. p_s^{(o,s')} are the
probabilities. col0 = [p_0^{(0,0)}, p_0^{(0,1)}, p_0^{(1,0)}, p_0^{(1,1)}] col1 =
[p_1^{(0,0)}, p_1^{(0,1)}, p_1^{(1,0)}, p_1^{(1,1)}]  Args:     a, b, c: Parameters for
column 0. The fourth, p1_10, is 1 - a - b - c.     e, f, g: Parameters for column 1. The
fourth, p1_11, is 1 - e - f - g.  Returns:     A tuple of symbolic matrices (M0, M1).
----------------------------------------------------------------------------------------
🟢counts_to_probabilities(counts_input: Union[Dict[str, int], List[Dict[str, int]]], max_len: int = None) -> Dict[str, float]
Convert count dictionary/dictionaries to a single probability dictionary.  This function
takes either a single dictionary or a list of dictionaries containing bit-string counts,
aggregates all counts, and normalizes them to probabilities. Probabilities for shorter
strings are computed by marginalizing (summing) over the longer strings that extend
them.  Args:     counts_input: Either a single dict mapping bit-strings to counts, or a
list of such dicts.     max_len: Maximum length of binary strings to include. If None,
inferred from input data.  Returns:     Dictionary mapping all binary strings (up to
max_len) to their probabilities.     Shorter strings' probabilities are computed by
summing over extensions.  Example:     >>> counts1 = {'00': 80, '01': 20}     >>>
counts2 = {'00': 70, '10': 30}     >>> counts_to_probabilities([counts1, counts2])
{'0': 0.75, '1': 0.25, '00': 0.75, '01': 0.1, '10': 0.15, '11': 0.0}
----------------------------------------------------------------------------------------
🟢create_list_of_circuits(max_word_length, num_qubits=2, visualize_first=4)
(no docstring)
----------------------------------------------------------------------------------------
🟢create_list_of_circuits_v2(max_word_length)
(no docstring)
----------------------------------------------------------------------------------------
🟢derived_constraints_from_empirical_probs(emp_probs: Dict[str, float]) -> Tuple[float, float, float, float, None | float, None | float, float, float]
Derive the constraint values from empirical probabilities. Args:     emp_probs:
Dictionary of empirical probabilities for binary strings. Returns:     Tuple of derived
constraint values:         - trace_M0         - det_M0         - trace_M1         -
det_M1         - trace_M0M1 and trace_M1M0 (they are the same but we can use both to
check consistency)         - prob_string_0 (need assumption of v0 being in maximally-
mixed state)         - prob_string_1 (need assumption of v0 being in maximally-mixed
state)
----------------------------------------------------------------------------------------
🟢display_job_details(job)
Display detailed information about a Qiskit Runtime job.  Args:     job: Qiskit Runtime
Job object
----------------------------------------------------------------------------------------
🟢extract_counts_by_qubit_and_pauli(job_result, num_qubits: int, pauli_labels: Tuple[str, ...] = ('I', 'X', 'Y', 'Z')) -> Dict[Tuple[int, str], Dict[str, int]]
From a SamplerV2 PrimitiveResult, reconstruct:     (qubit_index, Pauli_label) -> counts
over repeated-MCM bitstrings.  Assumes:     - Circuits were created with
create_list_of_circuits(max_word_length, num_qubits)     - Circuits are ordered as:
for q in range(num_qubits):             for gate_name in pauli_labels:
...     - For each circuit, the classical register used for qubit q is the       q-th
entry in pub_result.data.keys() (as in your plot_job_results).  Returns -------
result_counts : dict[(int,str), dict[str,int]]     Example key: (0, "X") -> {'000000':
..., '111000': ..., ...}
----------------------------------------------------------------------------------------
🟢gauge_transform_instrument_numerically(M0: numpy.ndarray, M1: numpy.ndarray, t: float) -> Tuple[numpy.ndarray, numpy.ndarray]
Apply a numerical gauge transformation to the instrument matrices M0 and M1. M' = R^{-1}
@ M @ R where R is a 2x2 matrix [[1 - t, t], [t, 1 - t]]  Args:     M0: Numpy array for
M0.     M1: Numpy array for M1.     t: Float value for the gauge parameter in [0,1].
Returns:     Transformed matrices M0' and M1'.
----------------------------------------------------------------------------------------
🟢gauge_transform_instrument_symbolically(M0: sympy.matrices.dense.MutableDenseMatrix, M1: sympy.matrices.dense.MutableDenseMatrix, t: sympy.core.symbol.Symbol) -> Tuple[sympy.matrices.dense.MutableDenseMatrix, sympy.matrices.dense.MutableDenseMatrix]
Apply a symbolic gauge transformation to the instrument matrices M0 and M1. M' = R^{-1}
@ M @ R where R is a 2x2 matrix [[1 - t, t], [t, 1 - t]]  Args:     M0: Symbolic matrix
for M0.     M1: Symbolic matrix for M1.     t: Symbolic variable for the gauge
parameter. Returns:     Transformed matrices M0' and M1'.
----------------------------------------------------------------------------------------
🟢generate_and_cache_simulated_data(GT: __main__.Instrument2x2, num_seeds: int, shots_per_seed: int = 100000, L: int = 4, use_same_v0_for_all_seeds: None | numpy.ndarray = array([0.5, 0.5]), seed_for_reproduce: int = 123, chosen_dir: str = 'simulated_data') -> Tuple[Dict[str, Any], List[numpy.ndarray], List[Dict[str, float]]]
Generate and cache simulated data for multiple seeds. Args:     GT: Ground truth
Instrument2x2 instance.     num_seeds: Number of different random seeds to simulate.
shots_per_seed: Number of Monte Carlo shots per seed.     L: Maximum length of binary
strings to simulate.     use_same_v0_for_all_seeds: If provided, use this initial state
distribution for all seeds.     seed_for_reproduce: Seed for reproducibility of the
entire simulation.     chosen_dir: Directory to save the simulated data. Returns:
Tuple of:         - Experiment info dictionary.         - List of initial state
distributions (v0) for each seed.         - List of empirical probability dictionaries
for each seed.
----------------------------------------------------------------------------------------
🟢get_all_binary_strings(max_len: int) -> List[str]
Generates all binary strings up to a given maximum length. Args:     max_len: Maximum
length of binary strings to generate. Returns:     List of binary strings. Example:
get_all_binary_strings(2) returns ['0', '1', '00', '01', '10', '11']
----------------------------------------------------------------------------------------
🟢get_prob_meas_0_from_counts(result_counts: Dict[Tuple[int, str], Dict[str, int]], qubit_index: int, filter_paulis: Tuple[str, ...] = ('I', 'Z')) -> float
Extracts the probability of measuring '0' for a specific qubit from the result counts.
It aggregates counts from 'I' and 'Z' Pauli basis preparations.
----------------------------------------------------------------------------------------
🟢make_instrument_from_columns(col0: numpy.ndarray, col1: numpy.ndarray) -> __main__.Instrument2x2
Build Instrument2x2 from two columns of 4 probabilities each: col_s = [p_s^{(0,0)},
p_s^{(0,1)}, p_s^{(1,0)}, p_s^{(1,1)}], sums to 1. Args:     col0: Column for starting
state |0>     col1: Column for starting state |1> Returns:     Instrument2x2 instance
----------------------------------------------------------------------------------------
🟢monte_carlo_empirical_probabilities(inst: __main__.Instrument2x2, v0: numpy.ndarray, L: int, shots: int, seed: int = 1234) -> Dict[str, float]
Simulate the MCM process to obtain empirical probabilities for all binary strings up to
length L. Args:     inst: Instrument2x2 instance representing the MCM instrument.
v0: Initial state distribution as a numpy array of shape (2,).     L: Maximum length of
binary strings to observe.     shots: Number of Monte Carlo shots to simulate.     seed:
Random seed for reproducibility. Returns:     Dictionary mapping binary strings to their
empirical probabilities.
----------------------------------------------------------------------------------------
🟢plot_gauge_transformation_effects(MCM_to_transform, t_width_factor: float = 1.2, verbose: bool = False, MCM_reference: List[__main__.Instrument2x2] = None, resolution: float = 100000.0, p00_min: float = 0.5, margin_tol: float = 0.0)
Analyzes and plots the effect of a gauge transformation on an instrument.  This function
takes an instrument, applies a gauge transformation over a range of the gauge parameter
't', and plots how each of the 8 matrix entries evolves. It uses analytical methods to
determine valid gauge parameter regions.  Note: t=0.5 is excluded from the analysis as
the gauge transformation matrix is non-invertible at that point.  Args:
MCM_to_transform: An Instrument2x2 object to be transformed.     t_width_factor: Factor
to scale the plotting range around valid regions. Default 1.20.                     1.00
means plot exactly the valid regions, >1.00 adds padding.     verbose: If True, print
detailed information about valid intervals. Default False.     MCM_reference: List of
reference Instrument2x2 objects to compare RMSE against. Default [].     resolution:
Number of points to sample per unit t-range.     p00_min: Minimum value for M0[0,0] when
identifying focus regions. Default 0.5.     margin_tol: Margin tolerance for allowed
regions. Default 0.0.  Returns:     Tuple containing:         - Instrument2x2 object
constructed from center points of valid entry ranges, or None if no valid ranges.
- List of valid t-regions for the center_instrument (relative to itself).
----------------------------------------------------------------------------------------
🟢plot_job_results(job_result, num_qubits, backend_name=None)
Plot histogram comparing results across multiple circuits (PUBs), grouped by classical
register index.  Args:     job_result: PrimitiveResult object from Qiskit Runtime
num_qubits: Number of qubits (groups) to split the results into.     backend_name:
Optional name of the backend used for execution      Returns:     list: A list of length
num_qubits, where each element is a list of count dictionaries            for that
classical register index.
----------------------------------------------------------------------------------------
🟢random_instrument(corr_strength: float = 0.2, seed: int = 1, fidelity: float = 0.85) -> __main__.Instrument2x2
Sample a random but 'physical' instrument with a tunable correlation flavor. We bias
columns so that outcome=0 is more likely when s=0 and vice versa, and inject some back-
action to create realistic correlations. For high-fidelity cases, p_0^(0,0) and
p_1^(1,1) should be close to 1.  Args:     corr_strength: Tunes how much the Dirichlet
distribution is biased.     seed: Random seed for reproducibility.     fidelity: The
base probability for the 'correct' outcome (e.g., p_0(0,0)). Must be between 0.25 and
1.0. Returns:     Instrument2x2 instance
----------------------------------------------------------------------------------------
🟢rebase_and_anchor_instrument(M0, M1, t_regions, p00_min=0.5, tol=1e-24)
For each allowed t-interval [t_lo, t_hi], anchor at t_lo, transform (M0,M1), and return
the anchored pair with its adjusted local-gauge interval [0, s_max], where s_max = (t_hi
- t_lo)/(1 - 2*t_lo).  Parameters ---------- M0, M1 : (2,2) np.ndarray     Trial
instrument matrices. t_regions : list of (t_lo, t_hi)     Allowed t-intervals from
allowed_t_regions_for_list([M0, M1]). p00_min : float     Required lower bound for
anchored [M0']_{00}. tol : float     Numerical tolerance for boundary checks.  Returns
------- results : list of dict     Each dict contains:         - 'anchor_t' : float
- 'anchored_M0' : (2,2) np.ndarray         - 'anchored_M1' : (2,2) np.ndarray         -
'adjusted_region' : (0.0, s_max)  # can have s_max < 0 when anchor_t > 1/2         -
'ok_p00' : bool  # whether [anchored_M0]_{00} >= p00_min (within tol)         - 'p00' :
float
----------------------------------------------------------------------------------------
🟢reconstruct_instrument_from_invariants_mixed_det(trM0: float, detM0: float, trM1: float, detM1: float, S0: float, gauge_p00: float = 1.0, noise_tol: float = 1e-12) -> List[Tuple[numpy.ndarray, numpy.ndarray]]
Inputs: trM0, detM0, trM1, detM1, S0, and gauge M0[0,0]=a. the parametrization is:
M0 = [[a, b],           [c, d]]     M1 = [[e, f],           [g, h]] Solve M0 from:     d
= trM0 - a     b + c = S0 - trM0     b * c = a*d - detM0 Generally we will have two set
of solutions for (b,c). Then solve M1 from:     e+g = 1 - a - c   (= u)     f+h = 1 - b
- d   (= v)     e+h = trM1        (= T1)     eh - fg = detM1 The last equation gives h
linearly: h = [detM1 + v(u - T1)]/(u - v) unless u=v. (this is unlikely for general
MCMs) If u≈v, fall back to tr(M0M1) computed from detM1 by identity.
----------------------------------------------------------------------------------------
🟢refine_bounds_with_sp_error(reconstructed_MCM: Union[__main__.Instrument2x2, Any], t_bound: Union[Tuple[float, float], Any], prob_meas_0_val: Union[float, Any], margin_tol: float = 0.0001, verbose: bool = False) -> Tuple[Tuple[float, float], Tuple[float, float]]
Refines the valid gauge parameter interval (t_bound) by incorporating the constraint
that the state preparation error (epsilon) must be physical (within [0, 1]).  Logic:
P(outcome=1) = P(1|0_noisy)                  = P(1|0_ideal) * (1 - epsilon) +
P(1|1_ideal) * epsilon          Where P(1|s_ideal) are functions of the gauge parameter
t derived from reconstructed_MCM.     We solve for epsilon(t) and find the subset of
t_bound where 0 <= epsilon(t) <= 1.  Args:     reconstructed_MCM: The Instrument2x2
object at the center of the gauge (t=0).     t_bound: The current valid interval (t_min,
t_max) derived from MCM positivity.     prob_meas_0_val: The empirical probability of
measuring '0' (from prob dictionary).     margin_tol: Tolerance for the physical bounds
of epsilon (allows -tol to 1+tol).     verbose: If True, prints detailed intermediate
values.  Returns:     Tuple containing:         - refined_t_bound: (new_t_min,
new_t_max)         - epsilon_bound: (eps_min, eps_max) over the refined t interval.
----------------------------------------------------------------------------------------
🟢show_functions_with_docs(namespace=None, width=88)
(no docstring)
----------------------------------------------------------------------------------------
🟢summarize_instrument(M0: numpy.ndarray, M1: numpy.ndarray) -> Dict[str, Any]
(no docstring)
----------------------------------------------------------------------------------------
🟢user_functions(namespace=None)
(no docstring)
----------------------------------------------------------------------------------------