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(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:
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)
----------------------------------------
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)

    # 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 len(valid_t_regions) > 0:
            t_lo, t_hi = valid_t_regions[0]
            mask = (t_values >= t_lo) & (t_values <= t_hi)
            if np.any(mask):
                block_values = transformed_entries[mask, 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
    }

    quantity_valid_ranges = {}
    quantity_ranges_data = []
    for quantity_name, quantity_values in derived_quantities.items():
        if len(valid_t_regions) > 0:
            # Check first valid region
            t_lo, t_hi = valid_t_regions[0]
            mask = (t_values >= t_lo) & (t_values <= t_hi)
            if np.any(mask):
                block_values = quantity_values[mask]

                # We accept whatever values are in the valid region defined by matrix entries
                # (removed the check for derived quantity validity to ensure all are shown)

                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
        )

    return center_instrument, center_t_regions

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
)

if 1==2:
    for i, sol in enumerate(sols):
        regions = allowed_t_regions_for_list([sol[0], sol[1]], tol=1e-24)
        print(f"\nallowed t-regions for the {i+1}th solution pair:")
        display(regions)
        print("width of allowed t-regions:")
        display([reg[1] - reg[0] for reg in regions])

        # Suppose sols[k] = (M0_trial, M1_trial) and regions = allowed_t_regions_for_list([M0_trial, M1_trial])
        res = rebase_and_anchor_instrument(sol[0], sol[1], regions, p00_min=0.5, tol=1e-24)

        for item in res:
            if item['ok_p00']:
                print("------------------------------------------------------------")
                print("anchor t =", item['anchor_t'])
                print("adjusted region (local) =", item['adjusted_region'])
                print("p00 anchored =", item['p00'], "ok?", item['ok_p00'])
                # item['anchored_M0'], item['anchored_M1'] are the rebased instrument matrices
                display(item['anchored_M0'], item['anchored_M1'])

        print("============================================================")

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.00003190	0.00127620	1.24430e-03	(6.541 ± 6.222)e-4	IN (+6.08221e-04)
1	prep 1 meas 0	0.00027571	0.00152002	1.24430e-03	(8.979 ± 6.222)e-4	IN (+6.08221e-04)
2	prep 0 excite to 1	0.00124866	0.00124972	1.06778e-06	(1.249 ± 0.001)e-3	IN (-5.21907e-07)
3	prep 1 decay to 0	0.00210546	0.00210653	1.06778e-06	(2.106 ± 0.001)e-3	IN (+5.21907e-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.45140e-07	(9.987 ± 0.000)e-1	IN (+2.04138e-07)
1	M^0[0,1] = p_1^(0,0)	0.00009245	0.00133686	1.24441e-03	(7.147 ± 6.222)e-4	IN (+6.08258e-04)
2	M^0[1,0] = p_0^(0,1)	0.00000346	0.00124787	1.24441e-03	(6.257 ± 6.222)e-4	IN (-6.08258e-04)
3	M^0[1,1] = p_1^(0,1)	0.00018315	0.00018360	4.45140e-07	(1.834 ± 0.002)e-4	IN (-2.04138e-07)
4	M^1[0,0] = p_0^(1,0)	0.00003005	0.00003106	1.01354e-06	(3.056 ± 0.051)e-5	IN (+4.56558e-07)
5	M^1[0,1] = p_1^(1,0)	0.00076966	0.00201301	1.24335e-03	(1.391 ± 0.622)e-3	IN (-6.07736e-04)
6	M^1[1,0] = p_0^(1,1)	0.00000185	0.00124520	1.24335e-03	(6.235 ± 6.217)e-4	IN (+6.07736e-04)
7	M^1[1,1] = p_1^(1,1)	0.99771026	0.99771127	1.01354e-06	(9.977 ± 0.000)e-1	IN (-4.56558e-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_4263/1689714025.py:509: 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.99872011, 0.00071466]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00062567, 0.00018338]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00003056, 0.00139134]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00062352, 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.00003190	0.00127620	1.24430e-03	(6.541 ± 6.222)e-4	IN (+6.08221e-04)
1	prep 1 meas 0	0.00027571	0.00152002	1.24430e-03	(8.979 ± 6.222)e-4	IN (+6.08221e-04)
2	prep 0 excite to 1	0.00124866	0.00124972	1.06778e-06	(1.249 ± 0.001)e-3	IN (-5.21907e-07)
3	prep 1 decay to 0	0.00210546	0.00210653	1.06778e-06	(2.106 ± 0.001)e-3	IN (+5.21907e-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.45140e-07	(9.987 ± 0.000)e-1	IN (+2.04138e-07)
1	M^0[0,1] = p_1^(0,0)	0.00009245	0.00133686	1.24441e-03	(7.147 ± 6.222)e-4	IN (+6.08258e-04)
2	M^0[1,0] = p_0^(0,1)	0.00000346	0.00124787	1.24441e-03	(6.257 ± 6.222)e-4	IN (-6.08258e-04)
3	M^0[1,1] = p_1^(0,1)	0.00018315	0.00018360	4.45140e-07	(1.834 ± 0.002)e-4	IN (-2.04138e-07)
4	M^1[0,0] = p_0^(1,0)	0.00003005	0.00003106	1.01354e-06	(3.056 ± 0.051)e-5	IN (+4.56558e-07)
5	M^1[0,1] = p_1^(1,0)	0.00076966	0.00201301	1.24335e-03	(1.391 ± 0.622)e-3	IN (-6.07736e-04)
6	M^1[1,0] = p_0^(1,1)	0.00000185	0.00124520	1.24335e-03	(6.235 ± 6.217)e-4	IN (+6.07736e-04)
7	M^1[1,1] = p_1^(1,1)	0.99771026	0.99771127	1.01354e-06	(9.977 ± 0.000)e-1	IN (-4.56558e-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_4263/1689714025.py:509: 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.99872011, 0.00071466]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00062567, 0.00018338]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00003056, 0.00139134]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00062352, 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.3, 
    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.00830149	0.01103481	2.73332e-03	(9.668 ± 1.367)e-3	IN (-1.19783e-03)
1	prep 1 meas 0	0.00686036	0.00959368	2.73332e-03	(8.227 ± 1.367)e-3	IN (-1.25409e-03)
2	prep 0 excite to 1	0.02658732	0.02665272	6.54041e-05	(2.662 ± 0.003)e-2	OUT below (-3.89935e-05)
3	prep 1 decay to 0	0.00308704	0.00315245	6.54041e-05	(3.120 ± 0.033)e-3	IN (+3.00934e-05)

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

	Entry	Min	Max	Width	Value	Ref1
0	M^0[0,0] = p_0^(0,0)	0.96505624	0.96512162	6.53810e-05	(9.651 ± 0.000)e-1	IN (+1.16419e-05)
1	M^0[0,1] = p_1^(0,0)	0.00041697	0.00308491	2.66794e-03	(1.751 ± 1.334)e-3	IN (-1.20505e-03)
2	M^0[1,0] = p_0^(0,1)	0.02390895	0.02657689	2.66794e-03	(2.524 ± 0.133)e-2	IN (+1.18618e-03)
3	M^0[1,1] = p_1^(0,1)	0.00644339	0.00650877	6.53810e-05	(6.476 ± 0.033)e-3	OUT below (-4.90429e-05)
4	M^1[0,0] = p_0^(1,0)	0.00829104	0.00829296	1.91334e-06	(8.292 ± 0.001)e-3	OUT above (+2.64065e-05)
5	M^1[0,1] = p_1^(1,0)	0.00000214	0.00273548	2.73334e-03	(1.369 ± 1.367)e-3	IN (+1.23515e-03)
6	M^1[1,0] = p_0^(1,1)	0.00001043	0.00274377	2.73334e-03	(1.377 ± 1.367)e-3	IN (-1.22518e-03)
7	M^1[1,1] = p_1^(1,1)	0.99040227	0.99040418	1.91334e-06	(9.904 ± 0.000)e-1	OUT above (+1.98946e-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_4263/1689714025.py:509: 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.96508893, 0.00175094]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.02524292, 0.00647608]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00829200, 0.00136881]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00137710, 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.00832328	0.01107024	2.74696e-03	(9.697 ± 1.373)e-3	IN (-1.22644e-03)
1	prep 1 meas 0	0.00682588	0.00957284	2.74696e-03	(8.199 ± 1.373)e-3	IN (-1.22644e-03)
2	prep 0 excite to 1	0.02657753	0.02664323	6.57016e-05	(2.661 ± 0.003)e-2	IN (-2.93524e-05)
3	prep 1 decay to 0	0.00308764	0.00315334	6.57016e-05	(3.120 ± 0.033)e-3	IN (+2.93524e-05)

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

	Entry	Min	Max	Width	Value	Ref1
0	M^0[0,0] = p_0^(0,0)	0.96503876	0.96510446	6.56970e-05	(9.651 ± 0.000)e-1	IN (+2.89621e-05)
1	M^0[0,1] = p_1^(0,0)	0.00040273	0.00308399	2.68126e-03	(1.743 ± 1.341)e-3	IN (-1.19747e-03)
2	M^0[1,0] = p_0^(0,1)	0.02389100	0.02657226	2.68126e-03	(2.523 ± 0.134)e-2	IN (+1.19747e-03)
3	M^0[1,1] = p_1^(0,1)	0.00642315	0.00648885	6.56970e-05	(6.456 ± 0.033)e-3	IN (-2.89621e-05)
4	M^1[0,0] = p_0^(1,0)	0.00831801	0.00831994	1.92308e-06	(8.319 ± 0.001)e-3	IN (-5.68947e-07)
5	M^1[0,1] = p_1^(1,0)	0.00000364	0.00275061	2.74696e-03	(1.377 ± 1.373)e-3	IN (+1.22683e-03)
6	M^1[1,0] = p_0^(1,1)	0.00000527	0.00275223	2.74696e-03	(1.379 ± 1.373)e-3	IN (-1.22683e-03)
7	M^1[1,1] = p_1^(1,1)	0.99042159	0.99042351	1.92308e-06	(9.904 ± 0.000)e-1	IN (+5.68947e-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_4263/1689714025.py:509: 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.96507161, 0.00174336]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.02523163, 0.00645600]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00831897, 0.00137713]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00137875, 0.99042255]]

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)
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()

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):

    circuits_for_local_test = []
    circuits_for_QPU = []
    qreg_q = QuantumRegister(1, 'q')
    creg_c = ClassicalRegister(max_word_length, 'c')

    circuit_0 = QuantumCircuit(qreg_q, creg_c)

    circuit_1 = QuantumCircuit(qreg_q, creg_c)
    circuit_1.x(0)

    circuit_2 = QuantumCircuit(qreg_q, creg_c)
    circuit_2.y(0)

    circuit_3 = QuantumCircuit(qreg_q, creg_c)
    circuit_3.z(0)


    circuits_for_local_test.extend([circuit_0, circuit_1, circuit_2, circuit_3])
    circuits_for_QPU.extend([circuit_0, circuit_1, circuit_2, circuit_3])

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

    for i in range(max_word_length-1,-1,-1):
        for circ in circuits_for_local_test:
            circ.append(Measure(f"measure_c{i}"), [0], [i])
        for circ in circuits_for_QPU:
            circ.append(MidCircuitMeasure(), [0], [i]) # it seems like we cannot name the MidCircuitMeasure


    # === Visualize the circuit ===
    print("Visualizing circuit for local test:")
    for circ in circuits_for_local_test:
        circ.draw('mpl', style="iqp")

    print("Visualizing circuit for QPU:")
    for circ in circuits_for_QPU:
        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 plot_job_results(job_result, backend_name=None):
    """
    Plot histogram comparing results across multiple circuits (PUBs).

    Args:
        job_result: PrimitiveResult object from Qiskit Runtime
        backend_name: Optional name of the backend used for execution
    """
    # Create a single histogram with grouped bars for all PUB results
    all_bitstrings = set()
    pub_counts_list = []

    # Collect all bitstrings and counts from each PUB
    for i, pub_result in enumerate(job_result):
        counts = pub_result.data.c.get_counts()
        pub_counts_list.append(counts)
        all_bitstrings.update(counts.keys())
        print(f"\nPUB {i} counts:")
        print(counts)

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

    # Prepare data for grouped bar chart
    num_pubs = len(job_result)
    colors = ['blue', 'orange', 'green', 'red', 'purple', 'brown']
    bar_width = 0.8 / num_pubs
    x = np.arange(len(sorted_bitstrings))

    # Create single plot with narrower width
    fig, ax = plt.subplots(figsize=(10, 6))

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

        # Add count labels above each bar
        for j, (bar, val) in enumerate(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)

    ax.set_xlabel('Bitstring')
    ax.set_ylabel('Counts')
    ax.set_title(f'Comparison of Results Across {num_pubs} Circuits')
    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 pub_counts_list

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 string length 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())

Available to run on:



[<IBMBackend('ibm_pittsburgh')>,
 <IBMBackend('ibm_fez')>,
 <IBMBackend('ibm_torino')>,
 <IBMBackend('ibm_kingston')>,
 <IBMBackend('ibm_marrakesh')>]



✅ But the ones equipped with mid-circuit measurement are:



[<IBMBackend('ibm_pittsburgh')>, <IBMBackend('ibm_kingston')>]


Service Usage:



{'instance_id': 'crn:v1:bluemix:public:quantum-computing:us-east:a/ed5d7d2fb3b249c6baea6864058116cb:e46b4031-d4ee-48a6-88cd-1ae20a49ec21::',
 'plan_id': '7f666d17-7893-47d8-bf9d-2b2389fc4dfc',
 'usage_consumed_seconds': 2646,
 'usage_period': {'start_time': '2025-10-28T03:56:09.980Z',
  'end_time': '2025-11-25T03:56:09.980Z'},
 'usage_allocation_seconds': 48000,
 'usage_limit_reached': False,
 'usage_remaining_seconds': 45354}

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

circuits_for_local, circuits_for_QPU = create_list_of_circuits(max_word_length)

Visualizing circuit for local test:
Visualizing circuit for QPU:
Visualizing circuit for QPU:


/Users/trainerblade/Documents/02_myDocs/venv/lib/python3.10/site-packages/qiskit/circuit/quantumcircuit.py:3943: RuntimeWarning: Cregbundle set to False since an instruction needs to refer to individual classical wire
  return circuit_drawer(

# 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:")

noisy simulation run on backend 'fake_manila' results:

Run Jobs

chosen_qubits = [
    4
    ]

# 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
    # "ibm_marrakesh", # Heron r2, many jobs pending
    # "ibm_pittsburgh", # Heron r3 & MCM capability, many jobs pending
    )

# 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]

# # 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

print(f"🖥️ Running on {backend.name},\non qubit(s) with index: {chosen_qubits},\na total of {len(circuits_for_QPU)} 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', # margin_tol = 0.0, perfect reconstruction on ibm_kingston
]

i_wanna_retrieve_job_with_id = "d4ii9iiv0j9c73e2cd3g" # 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}")
else:
    print("No job ID specified for retrieval., proceeding with existing job_result (if exists).")

Retrieved job result for job ID: d4ii9iiv0j9c73e2cd3g

Reconstruction

display_job_details(job)

# plot histogram from the job result
counting_list = plot_job_results(job_result, backend_name=backend.name)

Job ID: d4ii9iiv0j9c73e2cd3g
Backend: <IBMBackend('ibm_kingston')>
Status: DONE
Creation Time: 2025-11-24 21:38:50.606320-06:00
Tags: ['q4']
Metadata: {'execution_spans': ExecutionSpans([DoubleSliceSpan(<start='2025-11-25 03:38:53', stop='2025-11-25 03:40:44', size=400000>)])}

PUB 0 counts:
{'000000': 93348, '100000': 904, '000100': 921, '000001': 951, '111111': 765, '000010': 923, '010000': 887, '001000': 920, '001001': 7, '001111': 34, '110000': 27, '001010': 13, '100100': 12, '100001': 13, '000111': 52, '000110': 12, '000011': 45, '010001': 14, '111100': 20, '010100': 9, '100010': 6, '011111': 26, '111000': 15, '101000': 11, '000101': 7, '001011': 1, '111110': 16, '001100': 9, '111001': 1, '011000': 15, '101111': 1, '111101': 1, '110011': 1, '001110': 2, '100111': 2, '010010': 3, '111011': 2, '110111': 1, '010011': 1, '011100': 1, '110010': 1}

PUB 1 counts:
{'111111': 87548, '111110': 2165, '111100': 1976, '111000': 2062, '000000': 943, '100000': 2040, '110000': 1984, '111010': 32, '111011': 203, '011111': 104, '101111': 147, '111101': 230, '001111': 1, '110100': 22, '101011': 6, '110111': 163, '100111': 31, '110010': 18, '000001': 13, '100001': 31, '100011': 20, '100010': 20, '110110': 6, '110001': 37, '110011': 25, '101000': 25, '101101': 5, '000010': 8, '000100': 12, '010000': 13, '100100': 20, '111001': 54, '011000': 4, '101110': 6, '010111': 2, '011110': 3, '001000': 14, '110101': 4, '000111': 1, '100101': 1, '011100': 1}

PUB 2 counts:
{'111111': 87941, '000000': 975, '111000': 1895, '100000': 1926, '101111': 118, '111100': 1914, '111010': 30, '111110': 2101, '110000': 1975, '111101': 223, '110111': 184, '000010': 20, '111011': 208, '000001': 15, '110001': 30, '111001': 49, '110010': 21, '011111': 126, '011000': 5, '100111': 26, '101000': 23, '110100': 17, '100001': 25, '100100': 16, '110011': 27, '100011': 21, '010111': 4, '100010': 18, '110101': 4, '110110': 10, '001000': 9, '000100': 4, '101011': 3, '010000': 11, '000110': 1, '101110': 8, '101100': 3, '000111': 2, '010100': 1, '011110': 2, '011100': 2, '101101': 2, '101010': 2, '010010': 1, '000011': 1, '001111': 1}

PUB 3 counts:
{'000000': 92668, '000001': 1027, '001000': 899, '111111': 1063, '000100': 865, '100000': 929, '010000': 955, '000010': 908, '110000': 23, '000011': 125, '010010': 14, '000111': 117, '001111': 112, '011111': 115, '111000': 14, '111100': 17, '100001': 11, '011000': 15, '000101': 6, '111101': 4, '100100': 9, '001010': 9, '001100': 8, '111011': 3, '010100': 7, '010001': 7, '101000': 8, '001001': 8, '100010': 6, '100111': 3, '011110': 2, '011011': 2, '111110': 14, '110010': 2, '010111': 2, '000110': 11, '001110': 4, '001011': 1, '110111': 3, '011100': 2, '101111': 1, '010011': 1}
Metadata: {'execution_spans': ExecutionSpans([DoubleSliceSpan(<start='2025-11-25 03:38:53', stop='2025-11-25 03:40:44', size=400000>)])}

PUB 0 counts:
{'000000': 93348, '100000': 904, '000100': 921, '000001': 951, '111111': 765, '000010': 923, '010000': 887, '001000': 920, '001001': 7, '001111': 34, '110000': 27, '001010': 13, '100100': 12, '100001': 13, '000111': 52, '000110': 12, '000011': 45, '010001': 14, '111100': 20, '010100': 9, '100010': 6, '011111': 26, '111000': 15, '101000': 11, '000101': 7, '001011': 1, '111110': 16, '001100': 9, '111001': 1, '011000': 15, '101111': 1, '111101': 1, '110011': 1, '001110': 2, '100111': 2, '010010': 3, '111011': 2, '110111': 1, '010011': 1, '011100': 1, '110010': 1}

PUB 1 counts:
{'111111': 87548, '111110': 2165, '111100': 1976, '111000': 2062, '000000': 943, '100000': 2040, '110000': 1984, '111010': 32, '111011': 203, '011111': 104, '101111': 147, '111101': 230, '001111': 1, '110100': 22, '101011': 6, '110111': 163, '100111': 31, '110010': 18, '000001': 13, '100001': 31, '100011': 20, '100010': 20, '110110': 6, '110001': 37, '110011': 25, '101000': 25, '101101': 5, '000010': 8, '000100': 12, '010000': 13, '100100': 20, '111001': 54, '011000': 4, '101110': 6, '010111': 2, '011110': 3, '001000': 14, '110101': 4, '000111': 1, '100101': 1, '011100': 1}

PUB 2 counts:
{'111111': 87941, '000000': 975, '111000': 1895, '100000': 1926, '101111': 118, '111100': 1914, '111010': 30, '111110': 2101, '110000': 1975, '111101': 223, '110111': 184, '000010': 20, '111011': 208, '000001': 15, '110001': 30, '111001': 49, '110010': 21, '011111': 126, '011000': 5, '100111': 26, '101000': 23, '110100': 17, '100001': 25, '100100': 16, '110011': 27, '100011': 21, '010111': 4, '100010': 18, '110101': 4, '110110': 10, '001000': 9, '000100': 4, '101011': 3, '010000': 11, '000110': 1, '101110': 8, '101100': 3, '000111': 2, '010100': 1, '011110': 2, '011100': 2, '101101': 2, '101010': 2, '010010': 1, '000011': 1, '001111': 1}

PUB 3 counts:
{'000000': 92668, '000001': 1027, '001000': 899, '111111': 1063, '000100': 865, '100000': 929, '010000': 955, '000010': 908, '110000': 23, '000011': 125, '010010': 14, '000111': 117, '001111': 112, '011111': 115, '111000': 14, '111100': 17, '100001': 11, '011000': 15, '000101': 6, '111101': 4, '100100': 9, '001010': 9, '001100': 8, '111011': 3, '010100': 7, '010001': 7, '101000': 8, '001001': 8, '100010': 6, '100111': 3, '011110': 2, '011011': 2, '111110': 14, '110010': 2, '010111': 2, '000110': 11, '001110': 4, '001011': 1, '110111': 3, '011100': 2, '101111': 1, '010011': 1}

Real hardware run on ibm_kingston

""" 
1e5 shots on ibm_kingston, on qubit #48
    Quantity            Min         Max         Width
0   prep 0 meas 1       0.00080893  0.00083010  2.11681e-05
1   prep 1 meas 0       0.00265893  0.00268010  2.11681e-05
2   prep 0 excite to 1  0.00003727  0.00003749  2.24144e-07
3   prep 1 decay to 0   0.01058908  0.01058931  2.24144e-07


after ~5 minutes, run learning again with all same settings:  on qubit #48
    Quantity            Min         Max         Width
0   prep 0 meas 1       0.00102677  0.00114363  1.16860e-04
1   prep 1 meas 0       0.00283677  0.00295363  1.16860e-04
2   prep 0 excite to 1  -0.00026578 -0.00026391 1.86942e-06
3   prep 1 decay to 0   0.01566764  0.01566951  1.86942e-06


run using 1e6 shots, on qubit #48
    Quantity            Min         Max         Width
0   prep 0 meas 1       0.00176879  0.00196512  1.96331e-04
1   prep 1 meas 0       0.00466479  0.00486112  1.96331e-04
2   prep 0 excite to 1  0.00019532  0.00020243  7.11004e-06
3   prep 1 decay to 0   0.03616983  0.03617694  7.11004e-06
"""

' \n1e5 shots on ibm_kingston, on qubit #48\n\tQuantity\t        Min\t        Max\t        Width\n0\tprep 0 meas 1\t    0.00080893\t0.00083010\t2.11681e-05\n1\tprep 1 meas 0\t    0.00265893\t0.00268010\t2.11681e-05\n2\tprep 0 excite to 1\t0.00003727\t0.00003749\t2.24144e-07\n3\tprep 1 decay to 0\t0.01058908\t0.01058931\t2.24144e-07\n\n\nafter ~5 minutes, run learning again with all same settings:  on qubit #48\n\tQuantity\t        Min\t        Max\t        Width\n0\tprep 0 meas 1\t    0.00102677\t0.00114363\t1.16860e-04\n1\tprep 1 meas 0\t    0.00283677\t0.00295363\t1.16860e-04\n2\tprep 0 excite to 1\t-0.00026578\t-0.00026391\t1.86942e-06\n3\tprep 1 decay to 0\t0.01566764\t0.01566951\t1.86942e-06\n\n\nrun using 1e6 shots, on qubit #48\n\tQuantity\t        Min\t        Max\t        Width\n0\tprep 0 meas 1\t    0.00176879\t0.00196512\t1.96331e-04\n1\tprep 1 meas 0\t    0.00466479\t0.00486112\t1.96331e-04\n2\tprep 0 excite to 1\t0.00019532\t0.00020243\t7.11004e-06\n3\tprep 1 decay to 0\t0.03616983\t0.03617694\t7.11004e-06\n'

9.82962e-04/2

0.000491481

chosen_sol_index = 0  # Choose which solution to analyze

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, _ = plot_gauge_transformation_effects(
    Instrument2x2(sols[chosen_sol_index][0], sols[chosen_sol_index][1]), 
    t_width_factor=1.3,
    verbose=False,
    MCM_reference=[],
    resolution=1e5,
    margin_tol=0.0000
    )

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

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

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

	Quantity	Min	Max	Width	Value
0	prep 0 meas 1	0.00947263	0.01045559	9.82962e-04	(9.964 ± 0.491)e-3
1	prep 1 meas 0	0.00142263	0.00240559	9.82962e-04	(1.914 ± 0.491)e-3
2	prep 0 excite to 1	0.00124021	0.00126086	2.06482e-05	(1.251 ± 0.010)e-3
3	prep 1 decay to 0	0.02199678	0.02201743	2.06482e-05	(2.201 ± 0.001)e-2

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

	Entry	Min	Max	Width	Value
0	M^0[0,0] = p_0^(0,0)	0.98953416	0.98953441	2.49129e-07	(9.895 ± 0.000)e-1
1	M^0[0,1] = p_1^(0,0)	0.00000065	0.00098361	9.82953e-04	(4.921 ± 4.915)e-4
2	M^0[1,0] = p_0^(0,1)	0.00001001	0.00099296	9.82953e-04	(5.015 ± 4.915)e-4
3	M^0[1,1] = p_1^(0,1)	0.00142198	0.00142223	2.49129e-07	(1.422 ± 0.000)e-3
4	M^1[0,0] = p_0^(1,0)	0.00920473	0.00922539	2.06575e-05	(9.215 ± 0.010)e-3
5	M^1[0,1] = p_1^(1,0)	0.02103383	0.02199613	9.62304e-04	(2.151 ± 0.048)e-2
6	M^1[1,0] = p_0^(1,1)	0.00026790	0.00123020	9.62304e-04	(7.491 ± 4.812)e-4
7	M^1[1,1] = p_1^(1,1)	0.97656058	0.97658124	2.06575e-05	(9.766 ± 0.000)e-1

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



/var/folders/dk/_dd6ng8n3yv_crq6vrdt8qzc0000gn/T/ipykernel_4263/1689714025.py:509: 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.98953429, 0.00049213]]
  [[p_0^(0,1), p_1^(0,1)]]   [[0.00050148, 0.00142210]]

M1 matrix (outcome 1):
  [[p_0^(1,0), p_1^(1,0)]] = [[0.00921506, 0.02151498]]
  [[p_0^(1,1), p_1^(1,1)]]   [[0.00074905, 0.97657091]]

To-Dos

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

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: 1.605224e-04

Number of strings compared: 126

Detailed comparison:

	String	Length	Theoretical	Empirical (ibm)	Difference	Abs Diff	Rel Error
0	0	1	0.49597500	0.49597500	+5.551115e-17	5.551115e-17	0.0000%
1	1	1	0.50402500	0.50402500	+0.000000e+00	0.000000e+00	0.0000%
2	00	2	0.49008262	0.49008250	+1.230631e-07	1.230631e-07	0.0000%
3	01	2	0.00589232	0.00589250	-1.823054e-07	1.823054e-07	0.0031%
4	10	2	0.01614733	0.01614750	-1.746473e-07	1.746473e-07	0.0011%
...	...	...	...	...	...	...	...
121	111011	6	0.00066002	0.00104000	-3.799846e-04	3.799846e-04	36.5370%
122	111100	6	0.01013983	0.00981750	+3.223254e-04	3.223254e-04	3.2832%
123	111101	6	0.00075312	0.00114500	-3.918810e-04	3.918810e-04	34.2254%
124	111110	6	0.01063791	0.01074000	-1.020882e-04	1.020882e-04	0.9505%
125	111111	6	0.44373018	0.44329250	+4.376845e-04	4.376845e-04	0.0987%

126 rows × 7 columns

================================================================================
RMSE by string length:
================================================================================
Length 1: RMSE = 3.925231e-17 (2 strings)
Length 2: RMSE = 1.819531e-07 (4 strings)
Length 3: RMSE = 1.492315e-04 (8 strings)
Length 4: RMSE = 1.959335e-04 (16 strings)
Length 5: RMSE = 1.708640e-04 (32 strings)
Length 6: RMSE = 1.541141e-04 (64 strings)

Nov 20 Discussion

0. Related works: bounding SP and M (individually / separately)

• https://arxiv.org/abs/2506.09131, from Senrui: learning SPAM noise on n-qudit assuming perfect single-qudit control. One message is there are exactly 2^n-1 independent gauge degrees of freedom, corresponding to depolarizing gauge acting on each non-empty subset of qudit.

---

https://arxiv.org/abs/2506.09131	This work
Requires perfect single-quDit control.	Requires perfect single-quBit control.
Bound/Learn all learnable Pauli fidelity exactly (up to shot noise) -> reconstruct error rates up to First Order.	Bound all error rates (non-unital) exactly (up to shot noise).
Gauge: depolarizing gauge	Gauge: continuum gauge under our parametrization (same 1 DoF)
Model the terminating error channel's M errors	Model the non-terminating error channel's M error (MCM): POVM is non-diagonal
Model the initial state's SP errors	Model the initial state's SP errors and MCM's back-action (similar to its SP errors)

one-way property (irreversible):

• if given the non-unital noise rates, we can calculate the Pauli error rates after twirling.

• if given the Pauli error rates, we cannot calculate the non-unital noise rates before twirling.

Additional for this work:

• learning circuits used are not cycle benchmarking circuits, just simple noisy state preparation + noisy consecutive MCMs.

• of all the unlocked devices, IBM Kingston and IBM Pittsburgh have MCM capabilities.

• IBM Pittsburgh uses Heron r3 processor, which is more advanced than Kingston's Heron r2 but also too popular and busy. (queue time is too long ~ 6 hours)

• so I choose IBM Kingston for experiments.

1. Choose one qubit from IBM Kingston with low MEASURE_2 error

• At Nov 20, 00:07 Chicago time, the qubit 48 has MEASURE_2 error 4.883e-4, which is the lowest among all qubits in IBM Kingston.

• At Nov 20, 00:07 Chicago time, the qubit 31 has MEASURE_2 error 7.983e-2, which is the higher among all qubits in IBM Kingston.

2. Choose qubit #48 and #4 and run learning protocol circuits

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

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

• for job d4fb9692bisc73a2bhd0, we used qubit #48 from IBM Kingston and the circuits with 1e6 shots each for all 4 circuits.

• check framework consistency/integrity (%), if our model is good enough to explain the data within shot noise, the RMSE for prob strings of length >=3 should increase at the rate of shot noise + single-qubit gate errors. (can see from below snapshot that Pauli-X error rate on qubit #48 is 2.082e-4)

  ==============================================
  RMSE by string length:
  ==============================================
  Length 1: RMSE = 7.850462e-17 (2 strings)
  Length 2: RMSE = 9.803064e-09 (4 strings)
  Length 3: RMSE = 1.418754e-04 (8 strings)
  Length 4: RMSE = 1.614584e-04 (16 strings)
  Length 5: RMSE = 1.722368e-04 (32 strings)
  Length 6: RMSE = 1.782331e-04 (64 strings)
  

# compare to 4.883e-4 of MEASURE_2 error from ibm's calibration data for qubit 48 


# 1e6 shots for each of the 4 circuits, perfect reconstruction with 0 margin_tol
"""
    Quantity            Min         Max         Width
0   prep 0 meas 1       0.00176795  0.00196864  2.00690e-04
1   prep 1 meas 0       0.00466395  0.00486464  2.00690e-04
2   prep 0 excite to 1  0.00019519  0.00020246  7.26793e-06
3   prep 1 decay to 0   0.03616980  0.03617706  7.26793e-06
"""


# 1e5 shots for each of the 4 circuits, perfect reconstruction with 0 margin_tol (first trial)
"""
    Quantity            Min         Max         Width
0   prep 0 meas 1       0.00080893  0.00083010  2.11681e-05
1   prep 1 meas 0       0.00265893  0.00268010  2.11681e-05
2   prep 0 excite to 1  0.00003727  0.00003749  2.24144e-07
3   prep 1 decay to 0   0.01058908  0.01058931  2.24144e-07
"""


# 1e5 shots for each of the 4 circuits, margin_tol = 0.00014 (second trial)
"""
    Quantity            Min         Max         Width
0   prep 0 meas 1       0.00108596  0.00108596  0.00000e+00
1   prep 1 meas 0       0.00289596  0.00289596  0.00000e+00
2   prep 0 excite to 1  -0.00026486 -0.00026486 0.00000e+00
3   prep 1 decay to 0   0.01566859  0.01566859  0.00000e+00
"""

'\n\tQuantity\t        Min\t        Max\t        Width\n0\tprep 0 meas 1\t    0.00108596\t0.00108596\t0.00000e+00\n1\tprep 1 meas 0\t    0.00289596\t0.00289596\t0.00000e+00\n2\tprep 0 excite to 1\t-0.00026486\t-0.00026486\t0.00000e+00\n3\tprep 1 decay to 0\t0.01566859\t0.01566859\t0.00000e+00\n'

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_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 string length 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)
(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
----------------------------------------------------------------------------------------
🟢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']
----------------------------------------------------------------------------------------
🟢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, backend_name=None)
Plot histogram comparing results across multiple circuits (PUBs).  Args:     job_result:
PrimitiveResult object from Qiskit Runtime     backend_name: Optional name of the
backend used for execution
----------------------------------------------------------------------------------------
🟢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.
----------------------------------------------------------------------------------------
🟢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)
----------------------------------------------------------------------------------------