WIP: Improvements to leakage modeling and gauge optimization

jupyter_notebooks/Examples/Leakage-automagic.ipynb

@@ -2,14 +2,21 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
-    "from pygsti.modelpacks import smq1Q_XY, smq1Q_ZN\n",
-    "from pygsti.tools.leakage import leaky_qubit_model_from_pspec, construct_leakage_report\n",
+    "%load_ext autoreload\n",
+    "%autoreload 2\n",
+    "import pygsti\n",
+    "from pygsti.modelpacks import smq1Q_XYI as mp\n",
+    "from pygsti.tools.leakage import leaky_qubit_model_from_pspec, construct_leakage_report, add_lago_models, gate_leakage_profile\n",
+    "from pygsti.report import construct_standard_report\n",
+    "from pygsti.tools.optools import unitary_to_superop\n",
     "from pygsti.data import simulate_data\n",
-    "from pygsti.protocols import StandardGST, ProtocolData"
+    "from pygsti.protocols import StandardGST, ProtocolData, GateSetTomography, GSTGaugeOptSuite, GSTBadFitOptions\n",
+    "import numpy as np\n",
+    "import scipy.linalg as la"
    ]
   },
   {
@@ -23,30 +30,222 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
-    "mp = smq1Q_XY\n",
-    "ed = mp.create_gst_experiment_design(max_max_length=32)\n",
-    "tm3 = leaky_qubit_model_from_pspec(mp.processor_spec(), mx_basis='l2p1')\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": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/modelmembers/povms/__init__.py:615: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
+      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
+      "  errgen_vec = _np.linalg.lstsq(phys_directions, soln.x)[0]\n"
+     ]
+    }
+   ],
+   "source": [
+    "ed = mp.create_gst_experiment_design(max_max_length=16)\n",
+    "ps = mp.processor_spec()\n",
+    "tm3 = leaky_qubit_model_from_pspec(ps, mx_basis='l2p1')\n",
+    "tm3.convert_members_inplace('CPTPLND')\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",
+    "se = 'multinomial'\n",
+    "leaky_gate = ('Gxpi2', 0)\n",
+    "dm3, leaking_state = with_leaky_gate(tm3, leaky_gate, strength=0.125)\n",
+    "if se == 'none':\n",
+    "    import pygsti.objectivefns.objectivefns as pg_ofns\n",
+    "    # Forward simulation can raise errors if there are positive\n",
+    "    # probabilities < pg_ofns.DEFAULT_MIN_PROB_CLIP.\n",
+    "    # We change that default here so forward simulation works with\n",
+    "    # the ideal model.\n",
+    "    pg_ofns.DEFAULT_MIN_PROB_CLIP = 1e-12\n",
+    "ds = simulate_data(dm3, ed.all_circuits_needing_data, num_samples=10_000, sample_error=se, seed=1997)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "--- Iterative GST: [--------------------------------------------------] 0.0%  92 circuits ---\r"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/objectivefns/objectivefns.py:4494: RuntimeWarning: divide by zero encountered in divide\n",
+      "  p5over_lsvec = 0.5/lsvec\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "--- Iterative GST: [##########----------------------------------------] 20.0%  168 circuits ---\r"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "--- Iterative GST: [####################------------------------------] 40.0%  285 circuits ---\r"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "--- Iterative GST: [##############################--------------------] 60.0%  448 circuits ---\r"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "--- Iterative GST: [########################################----------] 80.0%  616 circuits ---\r"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n",
+      "\n",
+      "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n",
+      "\n",
+      "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "--- Iterative GST: [##################################################] 100.0%  616 circuits ---\n"
+     ]
+    }
+   ],
+   "source": [
+    "gopsuite = GSTGaugeOptSuite.cast('stdgaugeopt')\n",
+    "gopsuite.gaugeopt_target = tm3\n",
+    "bfo = {'actions': ['wildcard1d'], 'wildcard1d_reference': 'diamond distance', 'threshold': 0.0}\n",
+    "bfo = GSTBadFitOptions.cast(bfo)\n",
+    "gst = GateSetTomography(initial_model=dm3, gaugeopt_suite=gopsuite, verbosity=1, name='CPTPLND', badfit_options=bfo)\n",
     "pd = ProtocolData(ed, ds)\n",
     "res = gst.run(pd)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 5,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/models/model.py:145: UserWarning: Model.num_modeltest_params could not obtain number of *non-gauge* parameters - using total instead\n",
+      "  _warnings.warn((\"Model.num_modeltest_params could not obtain number of *non-gauge* parameters\"\n"
+     ]
+    }
+   ],
+   "source": [
+    "import copy\n",
+    "resc = copy.deepcopy(res)\n",
+    "add_lago_models(resc, 'CPTPLND', verbosity=1)\n",
+    "est = resc.estimates['CPTPLND']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Running idle tomography\n",
+      "Computing switchable properties\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/report/plothelpers.py:367: UserWarning:\n",
+      "\n",
+      "Max-model params (92) <= model params (158)!  Using k == 1.\n",
+      "\n",
+      "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/forwardsims/mapforwardsim.py:732: UserWarning:\n",
+      "\n",
+      "Generating dense process matrix representations of circuits or gates \n",
+      "can be inefficient and should be avoided for the purposes of forward \n",
+      "simulation/calculation of circuit outcome probability distributions \n",
+      "when using the MapForwardSimulator.\n",
+      "\n",
+      "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/report/plothelpers.py:367: UserWarning:\n",
+      "\n",
+      "Max-model params (92) <= model params (158)!  Using k == 1.\n",
+      "\n"
+     ]
+    }
+   ],
    "source": [
-    "report_dir = 'example_files/leakage-report-automagic'\n",
-    "report_object, updated_res = construct_leakage_report(res, title='easy leakage analysis!')\n",
+    "report_dir = f'example_files/leakage-report-automagic-251016-NleakNone-{se}'\n",
+    "report_object = construct_standard_report(resc, advanced_options={'n_leak': None}, title=f'fidelity-gop sample_error={se}')\n",
+    "# report_object, updated_res = construct_leakage_report(res, title='easy leakage analysis!')\n",
     "# ^ Each estimate in updated_res has a new gauge-optimized model.\n",
     "#   The gauge optimization was done to reflect how our target gates\n",
     "#   are only _really_ defined on the first two levels of our\n",