WIP: general leakage modeling

[rileyjmurray]

This PR promotes pygsti.tools.leakage from a monolithic 616-line module into the
top-level subpackage pygsti.leakage, and simultaneously generalizes leakage analysis
so that the computational subspace is inferred from the basis rather than supplied as
an explicit n_leak: int parameter.

This PR is still WIP. I plan on revisiting the following before merge.

• The contents of the Leakage-automagic-2Q.ipynb notebook.
• The name of the augment_for_leakage_modeling function, and the labels
  for elements of its returned Basis. Latest unimplemented thought is that if the i-th element of the new basis is (numerically) equal to the i-th element of the old basis then its label should stay the same.
• How n_leak > 0 is currently used as a flag in some functions to say
  "activate leakage modeling behavior that you can infer from the basis."
• The change to EmbeddedOp._iter_matrix_elements.

Subpackage layout

pygsti/tools/leakage.py is deleted. Its contents are reorganized into five modules:

Module	Contents
leakage/core.py	computational_effect, computational_superkets, computational_projector, augment_for_leakage_modeling
leakage/metrics.py	Choi-induced metrics, projected metrics, transport profiles
leakage/models.py	leaky_qubit_model_from_pspec, promote_bb_to_bt
leakage/gaugeopt.py	lagoified_gopparams_dicts, std_lago_gopsuite, add_lago_models
leakage/reports.py	construct_leakage_report

pygsti/__init__.py now imports leakage, making pygsti.leakage a first-class
subpackage alongside pygsti.algorithms, pygsti.models, etc.

pygsti/tools/__init__.py replaces from .leakage import * with explicit re-exports
from pygsti.leakage, so all pygsti.tools.* leakage names continue to resolve without
any changes to their callers.

Basis-aware leakage detection

The central design change is a new Basis.implies_leakage_modeling property. A basis
implies leakage modeling when it has a unique element whose label consists entirely of
the character 'I' (e.g., 'I', 'II', 'III') and that element is proportional to
an orthogonal projector of rank strictly less than the full Hilbert space dimension.

This property drives three new functions in leakage/core.py:

• computational_effect(basis) — extracts the Hermitian operator that orthogonally
  projects onto the computational subspace.
• computational_superkets(basis) — returns a column-orthonormal matrix whose columns
  span the computational subspace of superoperator space.
• computational_projector(basis) — the Hermicity-preserving superoperator that
  projects M[H] onto M[C].

The n_leak: int keyword argument is removed from gaugeopt_to_target,
_create_objective_fn, and _legacy_create_scalar_objective. Those functions now call
mx_basis.implies_leakage_modeling directly, and look up the projector via
computational_projector(mxBasis) when needed. The signatures of
subspace_jtracedist, diamond_distance_projection_model, leaky_maximum_trace_dist,
and diamonddist_to_leakfree_cptp are similarly updated.

Basis.is_hermitian() is a new method used internally by computational_superkets to
validate its input.

Basis construction: augment_for_leakage_modeling

augment_for_leakage_modeling(basis, E) takes any Hermitian basis and a real orthogonal
projector E, and returns a new ExplicitBasis that implies leakage modeling with
computational subspace C = range(E). Calling this function on the Gell-Mann basis
with E = np.diag([1, 1, 0]) basically gives you the l2p1 built-in basis.

The returned basis has these properties (k = rank(E)):

• Element 0 is proportional to E, labeled with all 'I's (longest available).
• Elements 0 ... k²−1 span M[C].
• Elements k² ... dim−1 span M[H] \\ M[C]. All are labeled 'L[orig_label]', except the
  final element, which is proportional to I − E and labeled 'L'.
• Every element has unit Frobenius norm.

Element selection uses pivot_indices_after_deflation (pivoted QR after orthogonal
deflation of a fixed vector) to prefer the original basis elements that are "most
supported" on each subspace.

Input validation: E must be real (‖imag(E)‖ ≤ 1e-10), Hermitian, and (proportional
to) an orthogonal projector; otherwise ValueError is raised.

Leakage gauge groups and new model functions

leaky_qubit_model_from_pspec now assigns a DirectSumUnitaryGroup as its default
gauge group. This group is a block-diagonal subgroup of U(3), decomposed as U(2) ⨁ U(1):
U(2) acts on the 2-dimensional computational subspace and U(1) acts on the 1-dimensional
leakage level. The two summands are specified to _direct_sum_unitary_group by their
bases in Hilbert-Schmidt space — 'pp' of dimension 4 for the U(2) factor, 'pp' of
dimension 1 for the U(1) factor — rather than by the dimension of the underlying Hilbert
space itself.

Supporting U(1) as a gauge-group factor requires two new classes in gaugegroup.py:

• U1Group — a one-parameter GaugeGroup whose elements are complex phases.
• U1GroupElement — implements transform_matrix, transform_matrix_inverse,
  from_vector, to_vector, and inverse.

DirectSumUnitaryGroup and DirectSumUnitaryGroupElement are extended to accept
U1Group / U1GroupElement components alongside the existing UnitaryGaugeGroup and
TrivialGaugeGroup types.

promote_bb_to_bt (qubit ⨂ qutrit) is a new function in leakage/models.py. Given a
2-qubit ExplicitOpModel, it produces a 6-dimensional model in a tensor-product basis
pp ⨂ l2p1 where the second qubit factor has an accessible leakage level. The function
handles the lifting of 2-qubit unitaries to 2-qudit unitaries and POVM construction with
four effects ('00', '01', '10', '11'). Its default gauge group is a plain
UnitaryGaugeGroup over the full 6-dimensional space.

The Leakage-automagic-2Q.ipynb notebook is updated to use promote_bb_to_bt and
demonstrates a complete 2-qubit leakage GST workflow.

Other changes

• matrixtools.is_projector — new utility function (used in computational_superkets
  and pop_transport_profile).
• jamiolkowski_iso — adds a fallback when basis.create_simple_equivalent() raises,
  which happens for leakage tensor-product bases.
• find_closest_unitary_opmx — gains an op_basis parameter (defaults to 'pp').
• EmbeddedOp — early return in _iter_matrix_elements when the embedded op fills
  the full state space, avoiding unnecessary iteration.
• UnitaryGaugeGroup.__init__ — raises ValueError when state_space.udim == 1;
  callers should use U1Group instead.
• gaugeopt_to_target — docstring consolidation: the Notes section that was split
  across two separate """...""" blocks is merged into one.
• Type annotation typos in gaugegroup.py corrected (_np.array → _np.ndarray,
  _np.ndarry → _np.ndarray).
• algorithms/core.py import cleanup: BasisLike imported from baseobjs.basis
  directly rather than via the star-import chain.

Tests

A new test/unit/leakage/ package adds unit tests for the subpackage's public API:

• test_core.py — Basis.implies_leakage_modeling, Basis.is_hermitian(),
  matrixtools.is_projector, computational_effect, computational_superkets, and
  computational_projector (both the 1-arg basis-aware form and the 3-arg explicit form,
  plus consistency between them); augment_for_leakage_modeling (structural properties,
  spanning of M[C] and M[C]⊥, label format, non-diagonal projector round-trip
  through computational_effect / computational_superkets, and input validation).
• test_metrics.py — choi_state (Hermitian and positive-semidefinite output),
  self-distance = 0 and self-fidelity = 1 for subspace_jtracedist,
  subspace_superop_fro_dist, and subspace_entanglement_fidelity; transport profiles
  (gate_leakage_profile, gate_seepage_profile, pop_transport_profile) verified
  against identity and explicitly-constructed leaky gates; subspace_diamonddist
  guarded with @needs_cvxpy.
• test_models.py — leaky_qubit_model_from_pspec (state space dimension,
  DirectSumUnitaryGroup gauge group, trace-1 prep, POVM completeness) and
  promote_bb_to_bt (6-dimensional state space, UnitaryGaugeGroup gauge group,
  four-effect POVM completeness, gate label preservation).

test_gaugeopt_correctness.py adds FindPerfectGauge_DirectSumGaugeGroupTester, which:

1. Builds a leaky-qubit model via leaky_qubit_model_from_pspec.
2. Applies a random gauge transformation using a block-diagonal unitary (2×2 block on
   the computational subspace, a U(1) phase on the leakage level).
3. Runs gaugeopt_to_target with method='L-BFGS-B' over a DirectSumUnitaryGroup.
4. Verifies that all gate metrics return to near zero after optimization.

The test runs with both a U1Group and a TrivialGaugeGroup for the leakage factor,
confirming that both configurations recover the gauge.

test_gaugegroup.py gains U1GroupTester, covering construction, the identity transform
at angle zero, unitarity for arbitrary angles, from_vector/to_vector round-trip, 2π
angle wrapping, and the inverse element.

Migration notes

Code that imported directly from pygsti.tools.leakage will need to update to
pygsti.leakage. Anything that accessed leakage names through pygsti.tools.* requires
no change.

Callers that passed n_leak to gaugeopt_to_target will get a ValueError at runtime
if the model's basis does not imply leakage modeling. The intended replacement is to use
a leakage basis (e.g., 'l2p1') and let Basis.implies_leakage_modeling drive the
behavior automatically.

[rileyjmurray]

Todo: add the following function. Maybe not strictly necessary, but I wrote it in another PR where it's distinctly less necessary.

def trace_effect(op: _np.ndarray, op_basis: BasisLike, on_space: SpaceT = 'HilbertSchmidt'):
    """
    Let `op` be the array representation of a superoperator G in `op_basis`,
    where G maps from and to the space of order-d Hermitian operators.
    
    The trace effect of G is the Heritian operator E that satifies

        trace(G(ρ)) = trace(E ρ) for all order-d Hermitian matrices ρ.

    If on_space='HilbertSchmidt', then this function returns a superket representation
    of E in `op_basis`. If on_space='Hilbert', then we return E itself.
    """
    d = int(_np.round(op.size ** 0.25))
    assert op.shape == (d**2, d**2)
    basis = op_basis if isinstance(op_basis, _Basis) else _Basis.cast(op_basis, dim=d**2)
    vecI = _bt.stdmx_to_vec(_np.eye(d), basis)
    vecE = op.T.conj() @ vecI
    if on_space == 'HilbertSchmidt':
        return vecE
    else:
        E = _bt.vec_to_stdmx(vecE, op_basis)
        return E

[rileyjmurray]

The test in test/integration/test_leakage_gst.py is failing due to numerical issues. The body of the test includes the following info in a long-form comment.

        Original results are shown below. We don't rely on the exact numbers here. What matters is
        qualitative aspects of how Gxpi2 and Gypi2 deviate from their respective targets. Since our
        data generating model only applied leakage to Gxpi2, a "good" result reports much more error
        in Gxpi2 than Gypi2. (It's not clear to me why stdgaugeopt lacks wildcard error.)

        Leakage-aware guage optimization.

            | Gate    | ent. infidelity | 1/2 trace dist | 1/2 diamond dist | Max TOP  | Unmodeled error |
            |---------|-----------------|----------------|------------------|----------|-----------------|
            | []      | 0.000009        | 0.002639   	 | 0.003724	        | 0.000014 | 0.00007         |
            | Gxpi2:0 | 0.000965        | 0.030946       | 0.040751	        | 0.001629 | 0.000764        |
            | Gypi2:0 | 0.000375        | 0.019367       | 0.025308         | 0.000531 | 0.000475        |

        
        Standard gauge optimization

            | Gate    | ent. infidelity | 1/2 trace dist | 1/2 diamond dist | Max TOP  |
            |---------|-----------------|----------------|------------------|----------|
            | []      | 0.000016        | 0.003644       | 0.005146         | 0.000026 |
            | Gxpi2:0 | 0.000620        | 0.024758       | 0.032616         | 0.001000 |
            | Gypi2:0 | 0.000611        | 0.024714       | 0.033323         | 0.001001 |

        We'll run tests with subspace entanglement infidelity.

            * For LAGO, infidelity of Gxpi2 is 2.57x larger than that of Gypi2;
              we'll test for a 2x difference.
            
            * For standard gauge optimization, Gxpi2 and Gypi2 have almost identical infidelities;
              we'll test for a factor 1.1x there.

The issue is that the factor 1.1x is, apparently, not enough to cover the change after merging in develop, using the latest numpy. The tests pass locally when I run with older numpy. I don't know exactly what changes in numpy / develop would cause the changes here, but that's the nature of integration tests. I'll probably regenerate the table and pick a new tolerance.

Tests pass as of March 24 with the attached conda environment. Unclear how this is different from what's on CI.
leakage_integration_ok_env.txt

[rileyjmurray]

Must address before merge.

[rileyjmurray]

Should open an issue and remove this comment from the PR changes.

[rileyjmurray]

Not quite the case. More accurate to say that there is at least one element whose label is a string consisting of one or more copies of I (and no other characters), and that the element with the longest such label is proportional to a real orthogonal projector on H.

[rileyjmurray]

Resolve before merge.

[rileyjmurray]

Potentially dubious change worth discussing: use the letter 'C' for the first label, then strings of the form 'C[lbl]' with existing labels lbl for the next k^2-1 labels, then strings of the form 'L[lbl]' with exsiting labels lbl, and finally just the letter 'L'.

[rileyjmurray]

Maybe add a TODO that we should revisit the question of where this lives the next time we're making changes to pygsti.report.factory.py. Could open a GitHub issue about that and point to this line so we don't lose track of it.

[rileyjmurray]

TODO: open a GitHub issue about np.NaN vs np.nan.

[rileyjmurray]

Need close review

[rileyjmurray]

Need close review

[rileyjmurray]

Need close review