Leakage modeling: bugfixes for wildcard error and gauge optimization

jupyter_notebooks/Examples/Leakage-automagic.ipynb

@@ -6,10 +6,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from pygsti.modelpacks import smq1Q_XY, smq1Q_ZN\n",
+    "from pygsti.modelpacks import smq1Q_XYI as mp\n",
     "from pygsti.tools.leakage import leaky_qubit_model_from_pspec, construct_leakage_report\n",
     "from pygsti.data import simulate_data\n",
-    "from pygsti.protocols import StandardGST, ProtocolData"
+    "from pygsti.protocols import StandardGST, ProtocolData\n",
+    "import numpy as np\n",
+    "import scipy.linalg as la"
    ]
   },
   {
@@ -27,14 +29,51 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "mp = smq1Q_XY\n",
-    "ed = mp.create_gst_experiment_design(max_max_length=32)\n",
+    "def with_leaky_gate(m, gate_label, strength):\n",
+    "    rng = np.random.default_rng(0)\n",
+    "    v = np.concatenate([[0.0], rng.standard_normal(size=(2,))])\n",
+    "    v /= la.norm(v)\n",
+    "    H = v.reshape((-1, 1)) @ v.reshape((1, -1))\n",
+    "    H *= strength\n",
+    "    U = la.expm(1j*H)\n",
+    "    m_copy = m.copy()\n",
+    "    G_ideal = m_copy.operations[gate_label]\n",
+    "    from pygsti.modelmembers.operations import ComposedOp, StaticUnitaryOp\n",
+    "    m_copy.operations[gate_label] = ComposedOp([G_ideal, StaticUnitaryOp(U, basis=m.basis)])\n",
+    "    return m_copy, v\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ed = mp.create_gst_experiment_design(max_max_length=8)\n",
+    "# ^ The default max length is small so we don't have to wait as long \n",
+    "#   for the GST fit (just for purposes of this notebook).\n",
     "tm3 = leaky_qubit_model_from_pspec(mp.processor_spec(), mx_basis='l2p1')\n",
-    "# ^ We could use basis = 'gm' instead of 'l2p1'. We prefer 'l2p1'\n",
-    "#   because it makes process matrices easier to interpret in leakage\n",
-    "#   modeling.\n",
-    "ds = simulate_data(tm3, ed.all_circuits_needing_data, num_samples=1000, seed=1997)\n",
-    "gst = StandardGST( modes=('CPTPLND',), target_model=tm3, verbosity=2)\n",
+    "# ^ Target model. \"Leaky\" is a bit of a misnomer here. The returned model\n",
+    "#   is simply a qutrit lift of the qubit model; leakage erorrs in the\n",
+    "#   qubit model can manifest as CPTP Markovian errors in the qutrit model.\n",
+    "dgm3, leaking_state = with_leaky_gate(tm3, ('Gxpi2', 0), strength=0.125)\n",
+    "# ^ Data generating model. \n",
+    "num_samples = 100_000\n",
+    "# ^ The number of samples is large to compensate for short circuit length.\n",
+    "#   Feel free to change the number of samples to something more \"realistic\"\n",
+    "#   if you'd like.\n",
+    "if num_samples > 10_000:\n",
+    "    from pygsti.objectivefns import objectivefns\n",
+    "    objectivefns.DEFAULT_MIN_PROB_CLIP = objectivefns.DEFAULT_RADIUS = 1e-12\n",
+    "    # ^ There are numerical thresholding rules in objective function evaluation\n",
+    "    #   that lead to errors when the number of samples is extremely large.\n",
+    "    #   The lines above change those thresholding rules to be appropriate in\n",
+    "    #   the unusual setting that is this notebook.\n",
+    "ds = simulate_data(dgm3, ed.all_circuits_needing_data, num_samples=num_samples, seed=1997)\n",
+    "gst = StandardGST(\n",
+    "    modes=('CPTPLND',), target_model=tm3, verbosity=2,\n",
+    "    badfit_options={'actions': ['wildcard1d'], 'threshold': 0.0}\n",
+    ")\n",
     "pd = ProtocolData(ed, ds)\n",
     "res = gst.run(pd)"
    ]

pygsti/modelmembers/operations/linearop.py

@@ -15,6 +15,7 @@
 from pygsti.baseobjs.opcalc import bulk_eval_compact_polynomials_complex as _bulk_eval_compact_polynomials_complex
 from pygsti.modelmembers import modelmember as _modelmember
 from pygsti.tools import optools as _ot
+from pygsti.tools import matrixtools as _mt
 from pygsti import SpaceT
 
 from typing import Any
@@ -416,11 +417,14 @@ def frobeniusdist_squared(self, other_op, transform=None, inv_transform=None) ->
         float
         """
         self_mx = self.to_dense("minimal")
+        other_mx = other_op.to_dense("minimal")
         if transform is not None:
             self_mx = self_mx @ transform
+            if isinstance(inv_transform, _mt.IdentityOperator):
+                other_mx = other_mx @ transform
         if inv_transform is not None:
             self_mx = inv_transform @ self_mx
-        return _ot.frobeniusdist_squared(self_mx, other_op.to_dense("minimal"))
+        return _ot.frobeniusdist_squared(self_mx, other_mx)
 
 
     def frobeniusdist(self, other_op, transform=None, inv_transform=None):

pygsti/modelmembers/povms/effect.py

@@ -14,6 +14,7 @@
 
 from pygsti.modelmembers import modelmember as _modelmember
 from pygsti.tools import optools as _ot
+from pygsti.tools import matrixtools as _mt
 from pygsti.baseobjs.opcalc import bulk_eval_compact_polynomials_complex as _bulk_eval_compact_polynomials_complex
 
 from typing import Any
@@ -142,9 +143,12 @@ def frobeniusdist_squared(self, other_spam_vec, transform=None, inv_transform=No
         float
         """
         vec = self.to_dense()
+        other_vec = other_spam_vec.to_dense()
         if transform is not None:
             vec = transform.T @ vec
-        return _ot.frobeniusdist_squared(vec, other_spam_vec.to_dense())
+            if isinstance(inv_transform, _mt.IdentityOperator):
+                other_vec = transform.T @ other_vec
+        return _ot.frobeniusdist_squared(vec, other_vec)
 
     def residuals(self, other_spam_vec, transform=None, inv_transform=None):
         """

pygsti/models/gaugegroup.py

@@ -817,7 +817,7 @@ class UnitaryGaugeGroup(OpGaugeGroupWithBasis):
         to specifying the value of `pygsti.evotypes.Evotype.default_evotype`.
     """
 
-    def __init__(self, state_space: _StateSpace, basis: Optional[Union[_Basis, str]],
+    def __init__(self, state_space: _StateSpace, basis: Union[_Basis, str],
                  evotype: Optional[Union[_Evotype, str]] = 'default'):
         state_space = _StateSpace.cast(state_space)
         evotype = _Evotype.cast(str(evotype), default_prefer_dense_reps=True)  # since we use deriv_wrt_params
@@ -1129,25 +1129,36 @@ def _from_nice_serialization(cls, state: dict):  # memo holds already de-seriali
 
 class DirectSumUnitaryGroup(GaugeGroup):
     """
-    A subgroup of the unitary group, where the unitary operators in the group all have a
-    shared block-diagonal structure.
+    A subgroup of the unitary group on a Hilbert space H, where H has a direct sum structure.
 
-    Example setting where this is useful:
-        The system's Hilbert space is naturally expressed as a direct sum, H = U ⨁ V,
-        and we want gauge optimization to preserve the natural separation between U and V.
+    Unitaries in this group are block diagonal and preserve the direct sum structure of H.
     """
 
     def __init__(self, subgroups: Tuple[Union[UnitaryGaugeGroup, TrivialGaugeGroup], ...],
                  basis, name="Direct sum gauge group"):
+
+        # Step 1. Get the size of the direct sum Hilbert space
+        udim = 0
+        for sg in subgroups:
+            udim += sg.state_space.udim
+
+        # Step 2. Infer and validate the provided basis.
+        if isinstance(basis, str):
+            basis = _Basis.cast(basis, udim**2)
+            # NOTE: Strings get cast to BuiltinBasis objects, which subclass LazyBasis.
+            # If a LazyBasis is instantiated with invalid dimensions (e.g., dim=9 when
+            # basis='pp'), then an exception will be raised when we try to actually use
+            # the basis, such as by accessing basis.labels.
+
+        if basis.dim != udim**2:
+            msg = """
+            Inconsistency between `basis.dim` and the dimension of the direct sum Hilbert
+            space implied by the `state_space` members of the provided subgroups.
+            """
+            raise ValueError(msg)
+
+        self.basis = basis
         self.subgroups = subgroups
-        if isinstance(basis, _Basis):
-            self.basis = basis
-        elif isinstance(basis, str):
-            udim = 0
-            for sg in subgroups:
-                udim += sg.state_space.udim
-            ss = _ExplicitStateSpace(['all'], [udim])
-            self.basis = _BuiltinBasis("std", ss)
         self._param_dims = _np.array([sg.num_params for sg in subgroups])
         self._num_params = _np.sum(self._param_dims)
         super().__init__(name)