diff --git a/jupyter_notebooks/Examples/Leakage-automagic.ipynb b/jupyter_notebooks/Examples/Leakage-automagic.ipynb index c586ff789..07f93c469 100644 --- a/jupyter_notebooks/Examples/Leakage-automagic.ipynb +++ b/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", diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index c332553fa..0dfaecf5d 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -20,34 +20,79 @@ from pygsti import tools as _tools from pygsti.tools import mpitools as _mpit from pygsti.tools import slicetools as _slct +from pygsti.models import ( + ExplicitOpModel as _ExplicitOpModel +) +from pygsti.models.model import OpModel as _OpModel from pygsti.models.gaugegroup import ( TrivialGaugeGroupElement as _TrivialGaugeGroupElement, GaugeGroupElement as _GaugeGroupElement ) -from typing import Callable, Union, Optional +from typing import Callable, Union, Optional, Any -def gaugeopt_to_target(model, target_model, item_weights=None, - cptp_penalty_factor=0, spam_penalty_factor=0, - gates_metric="frobenius", spam_metric="frobenius", - gauge_group=None, method='auto', maxiter=100000, - maxfev=None, tol=1e-8, oob_check_interval=0, - convert_model_to=None, return_all=False, comm=None, - verbosity=0, check_jac=False, n_leak=0): +class GaugeoptToTargetArgs: + """ + This class is basically a namespace. It was introduced to strip out tons of complexity + in gaugeopt_to_target(...) without breaking old code that might call gaugeopt_to_target. """ - Optimize the gauge degrees of freedom of a model to that of a target. - - Parameters - ---------- - model : Model - The model to gauge-optimize - - target_model : Model - The model to optimize to. The metric used for comparing models - is given by `gates_metric` and `spam_metric`. - item_weights : dict, optional + old_trailing_positional_args = ( + 'item_weights', 'cptp_penalty_factor', 'spam_penalty_factor', + 'gates_metric', 'spam_metric', 'gauge_group', 'method', + 'maxiter', 'maxfev', 'tol', 'oob_check_interval', + 'convert_model_to', 'return_all', 'comm', 'verbosity', + 'check_jac' + ) + + @staticmethod + def parsed_model(model: Union[_OpModel, _ExplicitOpModel], convert_model_to: Optional[Any]) -> _OpModel: + if convert_model_to is None: + return model + + if not isinstance(model, _ExplicitOpModel): + raise ValueError('Gauge optimization only supports model conversion for ExplicitOpModels.') + conversion_args = convert_model_to if isinstance(convert_model_to, (list, tuple)) else (convert_model_to,) + model_out = model.copy() # don't alter the original model's parameterization (this would be unexpected) + for args in conversion_args: + if isinstance(args, str): + model_out.convert_members_inplace(args, set_default_gauge_group=True) + elif isinstance(args, dict): + model_out.convert_members_inplace(**args) + else: + raise ValueError("Invalid `convert_model_to` arguments: %s" % str(args)) + return model_out + + @staticmethod + def parsed_method(target_model: Optional[_OpModel], method: str, gates_metric: str, spam_metric: str, n_leak: Optional[int]) -> str: + ls_mode_allowed = bool(target_model is not None + and gates_metric.startswith("frobenius") + and spam_metric.startswith("frobenius") + and isinstance(n_leak, int) and n_leak == -1) + # and model.dim < 64: # least squares optimization seems uneffective if more than 3 qubits + # -- observed by Lucas - should try to debug why 3 qubits seemed to cause trouble... + if method == "ls" and not ls_mode_allowed: + raise ValueError("Least-squares method is not allowed! Target" + " model must be non-None and frobenius metrics" + " must be used.") + elif method == "auto": + return 'ls' if ls_mode_allowed else 'L-BFGS-B' + else: + return method + + # The static dicts of default values are substituted in gaugeopt_to_target's kwargs. + # This is a safe thing to do because no invocation of "gaugeopt_to_target" within pyGSTi + # used positional arguments past target_model. + + create_objective_passthrough_kwargs : dict[str,Any] = dict( + item_weights=None, # gets cast to a dict. + cptp_penalty_factor=0, + spam_penalty_factor=0, + check_jac=False, + ) + """ + item_weights : dict Dictionary of weighting factors for gates and spam operators. Keys can be gate, state preparation, or POVM effect, as well as the special values "spam" or "gates" which apply the given weighting to *all* spam operators @@ -56,16 +101,56 @@ def gaugeopt_to_target(model, target_model, item_weights=None, "spam" values. The precise use of these weights depends on the model metric(s) being used. - cptp_penalty_factor : float, optional + cptp_penalty_factor : float If greater than zero, the objective function also contains CPTP penalty terms which penalize non-CPTP-ness of the gates being optimized. This factor multiplies these CPTP penalty terms. - spam_penalty_factor : float, optional + spam_penalty_factor : float If greater than zero, the objective function also contains SPAM penalty terms which penalize non-positive-ness of the state preps being optimized. This factor multiplies these SPAM penalty terms. + check_jac : bool + When True, check least squares analytic jacobian against finite differences. + """ + + gaugeopt_custom_passthrough_kwargs : dict[str, Any] = dict( + maxiter=100000, + maxfev=None, + tol=1e-8, + oob_check_interval=0, + verbosity=0, + ) + """ + maxiter : int + Maximum number of iterations for the gauge optimization. + + maxfev : int + Maximum number of function evaluations for the gauge optimization. + Defaults to maxiter. + + tol : float + The tolerance for the gauge optimization. + + oob_check_interval : int + If greater than zero, gauge transformations are allowed to fail (by raising + any exception) to indicate an out-of-bounds condition that the gauge optimizer + will avoid. If zero, then any gauge-transform failures just terminate the + optimization. + """ + + other_kwargs : dict[str, Any] = dict( + gates_metric="frobenius", # defines ls_mode_allowed, then passed to _create_objective_fn, + spam_metric="frobenius", # defines ls_mode_allowed, then passed to _create_objective_fn, + gauge_group=None, # used iff convert_model_to is not None, then passed to _create_objective_fn and gaugeopt_custom + method='auto', # validated, then passed to _create_objective_fn and gaugeopt_custom + return_all=False, # used in the function body (only for branching after passing to gaugeopt_custom) + convert_model_to=None, # passthrough to parsed_model. + comm=None, # passthrough to _create_objective_fn and gaugeopt_custom + n_leak=-1 # passthrough to parsed_objective and _create_objective_fn; default value is a flag. + ) + """ gates_metric : {"frobenius", "fidelity", "tracedist"}, optional The metric used to compare gates within models. "frobenius" computes the normalized sqrt(sum-of-squared-differences), with weights @@ -98,22 +183,6 @@ def gaugeopt_to_target(model, target_model, item_weights=None, - 'evolve' -- evolutionary global optimization algorithm using DEAP - 'brute' -- Experimental: scipy.optimize.brute using 4 points along each dimensions - maxiter : int, optional - Maximum number of iterations for the gauge optimization. - - maxfev : int, optional - Maximum number of function evaluations for the gauge optimization. - Defaults to maxiter. - - tol : float, optional - The tolerance for the gauge optimization. - - oob_check_interval : int, optional - If greater than zero, gauge transformations are allowed to fail (by raising - any exception) to indicate an out-of-bounds condition that the gauge optimizer - will avoid. If zero, then any gauge-transform failures just terminate the - optimization. - convert_model_to : str, dict, list, optional For use when `model` is an `ExplicitOpModel`. When not `None`, calls `model.convert_members_inplace(convert_model_to, set_default_gauge_group=False)` if @@ -130,58 +199,103 @@ def gaugeopt_to_target(model, target_model, item_weights=None, When not None, an MPI communicator for distributing the computation across multiple processors. - verbosity : int, optional - How much detail to send to stdout. + n_leak : int + Used in leakage modeling. If positive, this specifies how modify defintiions + of gate and SPAM metrics to reflect the fact that target gates do not have + well-defined actions outside the computational subspace. + """ - check_jac : bool - When True, check least squares analytic jacobian against finite differences + @staticmethod + def create_full_kwargs(args: tuple[Any,...], kwargs: dict[str, Any]): + full_kwargs = GaugeoptToTargetArgs.create_objective_passthrough_kwargs.copy() + full_kwargs.update(GaugeoptToTargetArgs.gaugeopt_custom_passthrough_kwargs) + full_kwargs.update(GaugeoptToTargetArgs.other_kwargs) + full_kwargs.update(kwargs) + + if extra_posargs := len(args) > 0: + msg = \ + f""" + Recieved {extra_posargs} positional arguments past `model` and `target_model`. + These should be passed as appropriate keyword arguments instead. This version + of pyGSTi will infer intended keyword arguments based on the legacy argument + positions. Future versions of pyGSTi will raise an error. + """ + _warnings.warn(msg) + for k,v in zip(GaugeoptToTargetArgs.old_trailing_positional_args, args): + full_kwargs[k] = v + + if full_kwargs['item_weights'] is None: + full_kwargs['item_weights'] = dict() + + return full_kwargs + + +def gaugeopt_to_target(model, target_model, *args, **kwargs): + """ + Legacy function to optimize the gauge degrees of freedom of a model to that of a target. + + Use of more than two positional arguments is deprecated; keyword arguments should be used + instead. See the GaugeoptToTargetArgs class for available keyword arguments and their + default values. Returns ------- - model : if return_all == False + model : if kwargs.get('return_all', False) == False + + (goodnessMin, gaugeMx, model) : if kwargs.get('return_all', False) == True - (goodnessMin, gaugeMx, model) : if return_all == True Where goodnessMin is the minimum value of the goodness function (the best 'goodness') found, gaugeMx is the gauge matrix used to transform the model, and model is the final gauge-transformed model. """ - - if item_weights is None: item_weights = {} - - ls_mode_allowed = bool(target_model is not None - and gates_metric.startswith("frobenius") - and spam_metric.startswith("frobenius")) - #and model.dim < 64: # least squares optimization seems uneffective if more than 3 qubits - # -- observed by Lucas - should try to debug why 3 qubits seemed to cause trouble... - - if method == "ls" and not ls_mode_allowed: - raise ValueError("Least-squares method is not allowed! Target" - " model must be non-None and frobenius metrics" - " must be used.") - if method == "auto": - method = 'ls' if ls_mode_allowed else 'L-BFGS-B' - - if convert_model_to is not None: - conversion_args = convert_model_to if isinstance(convert_model_to, (list, tuple)) else (convert_model_to,) - model = model.copy() # don't alter the original model's parameterization (this would be unexpected) - for args in conversion_args: - if isinstance(args, str): - model.convert_members_inplace(args, set_default_gauge_group=True) - elif isinstance(args, dict): - model.convert_members_inplace(**args) - else: - raise ValueError("Invalid `convert_model_to` arguments: %s" % str(args)) - objective_fn, jacobian_fn = _create_objective_fn( - model, target_model, item_weights, - cptp_penalty_factor, spam_penalty_factor, - gates_metric, spam_metric, method, comm, check_jac, n_leak) + full_kwargs = GaugeoptToTargetArgs.create_full_kwargs(args, kwargs) + """ + This function handles a strange situation where `target_model` can be None. - result = gaugeopt_custom(model, objective_fn, gauge_group, method, - maxiter, maxfev, tol, oob_check_interval, - return_all, jacobian_fn, comm, verbosity) + In this case, the objective function will only depend on `cptp_penality_factor` + and `spam_penalty_factor`; it's forbidden to use method == 'ls'; and the + returned OpModel might not have its basis set. + """ - #If we've gauge optimized to a target model, declare that the + # arg parsing: validating `method` + gates_metric = full_kwargs['gates_metric'] + spam_metric = full_kwargs['spam_metric'] + n_leak = full_kwargs['n_leak'] + method = full_kwargs['method'] + method = GaugeoptToTargetArgs.parsed_method( + target_model, method, gates_metric, spam_metric, n_leak + ) + + # arg parsing: (possibly) converting `model` + convert_model_to = full_kwargs['convert_model_to'] + model = GaugeoptToTargetArgs.parsed_model(model, convert_model_to) + + # actual work: constructing objective_fn and jacobian_fn + item_weights = full_kwargs['item_weights'] + cptp_penalty = full_kwargs['cptp_penalty_factor'] + spam_penalty = full_kwargs['spam_penalty_factor'] + comm = full_kwargs['comm'] + check_jac = full_kwargs['check_jac'] + objective_fn, jacobian_fn = _create_objective_fn( + model, target_model, item_weights, cptp_penalty, spam_penalty, + gates_metric, spam_metric, method, comm, check_jac, n_leak + ) + + # actual work: calling the (wrapper of the wrapper of the ...) optimizer + gauge_group = full_kwargs['gauge_group'] + maxiter = full_kwargs['maxiter'] + maxfev = full_kwargs['maxfev'] + tol = full_kwargs['tol'] + oob_check = full_kwargs['oob_check_interval'] + return_all = full_kwargs['return_all'] + verbosity = full_kwargs['verbosity'] + result = gaugeopt_custom( + model, objective_fn, gauge_group, method, maxiter, maxfev, + tol, oob_check, return_all, jacobian_fn, comm, verbosity + ) + + # If we've gauge optimized to a target model, declare that the # resulting model is now in the same basis as the target. if target_model is not None: newModel = result[-1] if return_all else result @@ -190,9 +304,14 @@ def gaugeopt_to_target(model, target_model, item_weights=None, return result -def gaugeopt_custom(model, objective_fn, gauge_group=None, +GGElObjective = Callable[[_GaugeGroupElement, bool], Union[float, _np.ndarray]] + +GGElJacobian = Union[None, Callable[[_GaugeGroupElement], _np.ndarray]] + + +def gaugeopt_custom(model, objective_fn: GGElObjective, gauge_group=None, method='L-BFGS-B', maxiter=100000, maxfev=None, tol=1e-8, - oob_check_interval=0, return_all=False, jacobian_fn=None, + oob_check_interval=0, return_all=False, jacobian_fn: Optional[GGElJacobian]=None, comm=None, verbosity=0): """ Optimize the gauge of a model using a custom objective function. @@ -203,8 +322,9 @@ def gaugeopt_custom(model, objective_fn, gauge_group=None, The model to gauge-optimize objective_fn : function - The function to be minimized. The function must take a single `Model` - argument and return a float. + The function to be minimized. The function must take a GaugeGroupElement + and a bool. If method == 'ls' then objective_fn must return an ndarray; if + method != 'ls' then objective_fn must return a float. gauge_group : GaugeGroup, optional The gauge group which defines which gauge trasformations are optimized @@ -356,20 +476,239 @@ def _call_jacobian_fn(gauge_group_el_vec): return newModel -GGElObjective = Callable[[_GaugeGroupElement,bool], Union[float, _np.ndarray]] +def _transform_with_oob_check(mdl, gauge_group_el, oob_check): + """ Helper function that sometimes checks if mdl.transform_inplace(gauge_group_el) fails. """ + mdl = mdl.copy() + if oob_check: + try: + mdl.transform_inplace(gauge_group_el) + except Exception as e: + raise ValueError("Out of bounds: %s" % str(e)) # signals OOB condition + else: + mdl.transform_inplace(gauge_group_el) + return mdl + + +def _gate_fidelity_targets(model, target_model): + from pygsti.report.reportables import eigenvalue_entanglement_infidelity + gate_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + for lbl in target_model.operations: + G_target = target_model.operations[lbl].to_dense() + G_curest = model.operations[lbl].to_dense() + t = 1 - eigenvalue_entanglement_infidelity(G_curest, G_target, model.basis) + t = _np.clip(t, a_min=0.0, a_max=1.0) + gate_fidelity_targets[lbl] = (G_target, t) + return gate_fidelity_targets + + +def _prep_fidelity_targets(model, target_model): + prep_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + for preplbl in target_model.preps: + rho_target = _tools.vec_to_stdmx( target_model.preps[preplbl].to_dense(), model.basis ) + rho_curest = _tools.vec_to_stdmx( model.preps[preplbl].to_dense(), model.basis ) + t = _tools.eigenvalue_fidelity(rho_curest, rho_target) + t = _np.clip(t, a_min=0.0, a_max=1.0) + prep_fidelity_targets[preplbl] = (rho_target, t) + return prep_fidelity_targets + + +def _povm_effect_fidelity_targets(model, target_model): + povm_fidelity_targets : dict[ _baseobjs.Label, tuple[dict, dict] ] = dict() + for povmlbl in target_model.povms: + M_target = target_model.povms[povmlbl] + M_curest = model.povms[povmlbl] + Es_target = {elbl : _tools.vec_to_stdmx(e_target.to_dense(), model.basis) for (elbl, e_target) in M_target.items() } + Es_curest = {elbl : _tools.vec_to_stdmx(e_target.to_dense(), model.basis) for (elbl, e_target) in M_curest.items() } + ts = {elbl : _tools.eigenvalue_fidelity(Es_target[elbl], Es_curest[elbl]) for elbl in M_target } + ts = {elbl : _np.clip(t, a_min=0.0, a_max=1.0) for (elbl, t) in ts.items()} + povm_fidelity_targets[povmlbl] = (Es_target, ts) + return povm_fidelity_targets + + +def _legacy_create_scalar_objective(model, target_model, + item_weights: dict[str,float], cptp_penalty_factor: float, spam_penalty_factor: float, + gates_metric: str, spam_metric: str, n_leak: Optional[int] + ) -> tuple[GGElObjective, GGElJacobian]: + opWeight = item_weights.get('gates', 1.0) + spamWeight = item_weights.get('spam', 1.0) + prepWeight = item_weights.get('prep', spamWeight) + mxBasis = model.basis -GGElJacobian = Union[None, Callable[[_GaugeGroupElement], _np.ndarray]] + #Use the target model's basis if model's is unknown + # (as it can often be if it's just come from an logl opt, + # since from_vector will clear any basis info) + if mxBasis.name == "unknown" and target_model is not None: + mxBasis = target_model.basis -def _create_objective_fn(model, target_model, item_weights: Optional[dict[str,float]]=None, - cptp_penalty_factor: float=0.0, spam_penalty_factor: float=0.0, - gates_metric="frobenius", spam_metric="frobenius", - method=None, comm=None, check_jac=False, n_leak=0) -> tuple[GGElObjective, GGElJacobian]: - """ - Creates the objective function and jacobian (if available) - for gaugeopt_to_target - """ - if item_weights is None: item_weights = {} + udim = int(_np.sqrt(mxBasis.dim)) + I = _tools.matrixtools.IdentityOperator() + if n_leak is None and mxBasis.first_element_is_identity: + P = I + elif n_leak is None: + P = _tools.superop_subspace_projector(mxBasis) + elif n_leak > 0: + P = _tools.superop_subspace_projector(udim - n_leak, udim, mxBasis) + else: + # this also handles the case n_leak==-1, which is a flag. + P = I + n_leak = 0 + transform_mx_arg = (P, I) + # ^ The semantics of this tuple are defined by the frobeniusdist function + # in the ExplicitOpModelCalc class. + + gate_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + if gates_metric == 'fidelity': + gate_fidelity_targets.update(_gate_fidelity_targets(model, target_model)) + + prep_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + povm_fidelity_targets : dict[ _baseobjs.Label, tuple[dict, dict] ] = dict() + if spam_metric == 'fidelity': + prep_fidelity_targets.update(_prep_fidelity_targets(model, target_model)) + povm_fidelity_targets.update(_povm_effect_fidelity_targets(model, target_model)) + + def _objective_fn(gauge_group_el: _GaugeGroupElement, oob_check: bool) -> float: + mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) + ret = 0 + + if cptp_penalty_factor > 0: + mdl.basis = mxBasis # set basis for jamiolkowski iso + cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) + ret += _np.sum(cpPenaltyVec) + + if spam_penalty_factor > 0: + mdl.basis = mxBasis + spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) + ret += _np.sum(spamPenaltyVec) + + if target_model is None: + return ret + + if "frobenius" in gates_metric: + if spam_metric == gates_metric: + val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) + else: + wts = item_weights.copy() + wts['spam'] = 0.0 + for k in wts: + if k in mdl.preps or k in mdl.povms: + wts[k] = 0.0 + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts, n_leak) + if "squared" in gates_metric: + val = val ** 2 + ret += val + + elif gates_metric == 'commutators': + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * _np.linalg.norm((mop - top) @ P)**2 + + elif gates_metric == "fidelity": + # Leakage-aware metrics NOT available + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + mop = mdl.operations[opLbl].to_dense() + top, t = gate_fidelity_targets[opLbl] + v = _tools.entanglement_fidelity(mop, top, mxBasis) + z = _np.abs(t - v) + ret += wt * z + + elif gates_metric == "tracedist": + # If n_leak==0, then subspace_jtracedist is just jtracedist. + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * _tools.subspace_jtracedist(top, mop, mxBasis, n_leak) + + else: + raise ValueError("Invalid gates_metric: %s" % gates_metric) + + if "frobenius" in spam_metric and gates_metric == spam_metric: + # We already handled SPAM error in this case. Just return. + return ret + + if "frobenius" in spam_metric: + # SPAM and gates can have different choices for squared vs non-squared. + wts = item_weights.copy(); wts['gates'] = 0.0 + for k in wts: + if k in mdl.operations or k in mdl.instruments: + wts[k] = 0.0 + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts) + if "squared" in spam_metric: + val = val ** 2 + ret += val + + elif spam_metric == 'commutators': + for preplabel, m_prep in mdl.preps.items(): + wt_prep = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + com = rhoMx1 @ rhoMx2 - rhoMx2 - rhoMx1 + ret += wt_prep * _np.linalg.norm(com)**2 + + for povmlbl in target_model.povms: + wt_povm = item_weights.get(povmlbl, spamWeight) + M_curest = model.povms[povmlbl] + M_target = target_model.povms[povmlbl] + for elbl in M_target: + E_curest = _tools.vec_to_stdmx(M_curest[elbl].to_dense(), mxBasis) + E_target = _tools.vec_to_stdmx(M_target[elbl].to_dense(), mxBasis) + com_mat = E_curest @ E_target - E_target @ E_curest + com_vec = P @ _tools.stdmx_to_vec(com_mat, mxBasis) + val = _np.linalg.norm(com_vec) + ret += wt_povm * val**2 + + elif spam_metric == "fidelity": + # Leakage-aware metrics NOT available + val_prep = 0.0 + for preplbl in target_model.preps: + wt_prep = item_weights.get(preplbl, prepWeight) + rho_curest = _tools.vec_to_stdmx(model.preps[preplbl].to_dense(), mxBasis) + rho_target, t = prep_fidelity_targets[preplbl] + v = _tools.fidelity(rho_curest, rho_target) + val_prep += wt * abs(t - v) + val_povm = 0.0 + for povmlbl in target_model.povms: + wt_povm = item_weights.get(povmlbl, spamWeight) + M_target = target_model.povms[povmlbl] + M_curest = model.povms[povmlbl] + Es_target, ts = povm_fidelity_targets[povmlbl] + for elbl in M_target: + t_e = ts[elbl] + E = _tools.vec_to_stdmx(M_curest[elbl].to_dense(), mxBasis) + v_e = _tools.fidelity(E, Es_target[elbl]) + val_povm += wt * abs(t_e - v_e) + ret += (val_prep + val_povm) + + elif spam_metric == "tracedist": + # Leakage-aware metrics NOT available. + for preplabel, m_prep in mdl.preps.items(): + wt = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + ret += wt * _tools.tracedist(rhoMx1, rhoMx2) + + for povmlabel in mdl.povms.keys(): + wt = item_weights.get(povmlabel, spamWeight) + ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) + else: + raise ValueError("Invalid spam_metric: %s" % spam_metric) + + return ret + + return _objective_fn, None + + +def _legacy_create_least_squares_objective(model, target_model, + item_weights: dict[str,float], cptp_penalty_factor: float, spam_penalty_factor: float, + gates_metric: str, spam_metric: str, comm: Optional[Any], check_jac: bool + ) -> tuple[GGElObjective, GGElJacobian]: + opWeight = item_weights.get('gates', 1.0) spamWeight = item_weights.get('spam', 1.0) mxBasis = model.basis @@ -380,350 +719,228 @@ def _create_objective_fn(model, target_model, item_weights: Optional[dict[str,fl if mxBasis.name == "unknown" and target_model is not None: mxBasis = target_model.basis - def _transform_with_oob_check(mdl, gauge_group_el, oob_check): - """ Helper function that sometimes checks if mdl.transform_inplace(gauge_group_el) fails. """ - mdl = mdl.copy() - if oob_check: - try: mdl.transform_inplace(gauge_group_el) - except Exception as e: - raise ValueError("Out of bounds: %s" % str(e)) # signals OOB condition - else: - mdl.transform_inplace(gauge_group_el) - return mdl + assert(gates_metric.startswith("frobenius") and spam_metric.startswith("frobenius")), \ + "Only 'frobenius' and 'frobeniustt' metrics can be used when `method='ls'`!" + assert(gates_metric == spam_metric) + frobenius_transform_target = bool(gates_metric == 'frobeniustt') # tt = "transform target" - if method == "ls": - # least-squares case where objective function returns an array of - # the before-they're-squared difference terms and there's an analytic jacobian + if frobenius_transform_target: + full_target_model = target_model.copy() + full_target_model.convert_members_inplace("full") # so we can gauge-transform the target model. + else: + full_target_model = None # in case it get's referenced by mistake - assert(gates_metric.startswith("frobenius") and spam_metric.startswith("frobenius")), \ - "Only 'frobenius' and 'frobeniustt' metrics can be used when `method='ls'`!" - assert(gates_metric == spam_metric) - frobenius_transform_target = bool(gates_metric == 'frobeniustt') # tt = "transform target" + def _objective_fn(gauge_group_el: _GaugeGroupElement, oob_check: bool) -> _np.ndarray: if frobenius_transform_target: - full_target_model = target_model.copy() - full_target_model.convert_members_inplace("full") # so we can gauge-transform the target model. + transformed = _transform_with_oob_check(full_target_model, gauge_group_el.inverse(), oob_check) + other = model else: - full_target_model = None # in case it get's referenced by mistake + transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) + other = target_model - def _objective_fn(gauge_group_el, oob_check): + residuals, _ = transformed.residuals(other, None, item_weights) + # We still the non-target model to be transformed and checked for these penalties + if cptp_penalty_factor > 0 or spam_penalty_factor > 0: if frobenius_transform_target: - transformed = _transform_with_oob_check(full_target_model, gauge_group_el.inverse(), oob_check) - other = model - else: transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) - other = target_model - - residuals, _ = transformed.residuals(other, None, item_weights) - # We still the non-target model to be transformed and checked for these penalties - if cptp_penalty_factor > 0 or spam_penalty_factor > 0: - if frobenius_transform_target: - transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) - - if cptp_penalty_factor > 0: - transformed.basis = mxBasis - cpPenaltyVec = _cptp_penalty(transformed, cptp_penalty_factor, transformed.basis) - else: cpPenaltyVec = [] # so concatenate ignores + if cptp_penalty_factor > 0: + transformed.basis = mxBasis + cpPenaltyVec = _cptp_penalty(transformed, cptp_penalty_factor, transformed.basis) + else: cpPenaltyVec = [] # so concatenate ignores - if spam_penalty_factor > 0: - transformed.basis = mxBasis - spamPenaltyVec = _spam_penalty(transformed, spam_penalty_factor, transformed.basis) - else: spamPenaltyVec = [] # so concatenate ignores + if spam_penalty_factor > 0: + transformed.basis = mxBasis + spamPenaltyVec = _spam_penalty(transformed, spam_penalty_factor, transformed.basis) + else: spamPenaltyVec = [] # so concatenate ignores - return _np.concatenate((residuals, cpPenaltyVec, spamPenaltyVec)) - else: - return residuals + return _np.concatenate((residuals, cpPenaltyVec, spamPenaltyVec)) + else: + return residuals - def _jacobian_fn(gauge_group_el): + def _jacobian_fn(gauge_group_el: _GaugeGroupElement) -> _np.ndarray: - #Penalty terms below always act on the transformed non-target model. - original_gauge_group_el = gauge_group_el + #Penalty terms below always act on the transformed non-target model. + original_gauge_group_el = gauge_group_el - if frobenius_transform_target: - gauge_group_el = gauge_group_el.inverse() - mdl_pre = full_target_model.copy() - mdl_post = mdl_pre.copy() - else: - mdl_pre = model.copy() - mdl_post = mdl_pre.copy() - mdl_post.transform_inplace(gauge_group_el) - - # Indices: Jacobian output matrix has shape (L, N) - start = 0 - d = mdl_pre.dim - N = gauge_group_el.num_params - L = mdl_pre.num_elements - - #Compute "extra" (i.e. beyond the model-element) rows of jacobian - if cptp_penalty_factor != 0: L += _cptp_penalty_size(mdl_pre) - if spam_penalty_factor != 0: L += _spam_penalty_size(mdl_pre) - - #Set basis for pentaly term calculation - if cptp_penalty_factor != 0 or spam_penalty_factor != 0: - mdl_pre.basis = mxBasis - mdl_post.basis = mxBasis - - jacMx = _np.zeros((L, N)) - - #Overview of terms: - # objective: op_term = (S_inv * gate * S - target_op) - # jac: d(op_term) = (d (S_inv) * gate * S + S_inv * gate * dS ) - # d(op_term) = (-(S_inv * dS * S_inv) * gate * S + S_inv * gate * dS ) - - # objective: rho_term = (S_inv * rho - target_rho) - # jac: d(rho_term) = d (S_inv) * rho - # d(rho_term) = -(S_inv * dS * S_inv) * rho - - # objective: ET_term = (E.T * S - target_E.T) - # jac: d(ET_term) = E.T * dS - - #Overview of terms when frobenius_transform_target == True). Note that the objective - #expressions are identical to the above except for an additional overall minus sign and S <=> S_inv. - - # objective: op_term = (gate - S * target_op * S_inv) - # jac: d(op_term) = -(dS * target_op * S_inv + S * target_op * -(S_inv * dS * S_inv) ) - # d(op_term) = (-dS * target_op * S_inv + S * target_op * (S_inv * dS * S_inv) ) - - # objective: rho_term = (rho - S * target_rho) - # jac: d(rho_term) = - dS * target_rho - - # objective: ET_term = (E.T - target_E.T * S_inv) - # jac: d(ET_term) = - target_E.T * -(S_inv * dS * S_inv) - # d(ET_term) = target_E.T * (S_inv * dS * S_inv) - - #Distribute computation across processors - allDerivColSlice = slice(0, N) - derivSlices, myDerivColSlice, derivOwners, mySubComm = \ - _mpit.distribute_slice(allDerivColSlice, comm) - if mySubComm is not None: - _warnings.warn("Note: more CPUs(%d)" % comm.Get_size() - + " than gauge-opt derivative columns(%d)!" % N) # pragma: no cover - - n = _slct.length(myDerivColSlice) - wrtIndices = _slct.indices(myDerivColSlice) if (n < N) else None - my_jacMx = jacMx[:, myDerivColSlice] # just the columns I'm responsible for - - # S, and S_inv are shape (d,d) - #S = gauge_group_el.transform_matrix - S_inv = gauge_group_el.transform_matrix_inverse - dS = gauge_group_el.deriv_wrt_params(wrtIndices) # shape (d*d),n - dS.shape = (d, d, n) # call it (d1,d2,n) - dS = _np.rollaxis(dS, 2) # shape (n, d1, d2) - assert(dS.shape == (n, d, d)) - - # --- NOTE: ordering here, with running `start` index MUST - # correspond to those in Model.residuals, which in turn - # must correspond to those in ForwardSimulator.residuals - which - # currently orders as: gates, simplified_ops, preps, effects. - - # -- LinearOperator terms - # ------------------------- - for lbl, G in mdl_pre.operations.items(): - # d(op_term) = S_inv * (-dS * S_inv * G * S + G * dS) = S_inv * (-dS * G' + G * dS) - # Note: (S_inv * G * S) is G' (transformed G) - wt = item_weights.get(lbl, opWeight) - left = -1 * _np.dot(dS, mdl_post.operations[lbl].to_dense('minimal')) # shape (n,d1,d2) - right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # shape (d1,n,d2) -> (n,d1,d2) + if frobenius_transform_target: + gauge_group_el = gauge_group_el.inverse() + mdl_pre = full_target_model.copy() + mdl_post = mdl_pre.copy() + else: + mdl_pre = model.copy() + mdl_post = mdl_pre.copy() + mdl_post.transform_inplace(gauge_group_el) + + # Indices: Jacobian output matrix has shape (L, N) + start = 0 + d = mdl_pre.dim + N = gauge_group_el.num_params + L = mdl_pre.num_elements + + #Compute "extra" (i.e. beyond the model-element) rows of jacobian + if cptp_penalty_factor != 0: L += _cptp_penalty_size(mdl_pre) + if spam_penalty_factor != 0: L += _spam_penalty_size(mdl_pre) + + #Set basis for pentaly term calculation + if cptp_penalty_factor != 0 or spam_penalty_factor != 0: + mdl_pre.basis = mxBasis + mdl_post.basis = mxBasis + + jacMx = _np.zeros((L, N)) + + #Overview of terms: + # objective: op_term = (S_inv * gate * S - target_op) + # jac: d(op_term) = (d (S_inv) * gate * S + S_inv * gate * dS ) + # d(op_term) = (-(S_inv * dS * S_inv) * gate * S + S_inv * gate * dS ) + + # objective: rho_term = (S_inv * rho - target_rho) + # jac: d(rho_term) = d (S_inv) * rho + # d(rho_term) = -(S_inv * dS * S_inv) * rho + + # objective: ET_term = (E.T * S - target_E.T) + # jac: d(ET_term) = E.T * dS + + #Overview of terms when frobenius_transform_target == True). Note that the objective + #expressions are identical to the above except for an additional overall minus sign and S <=> S_inv. + + # objective: op_term = (gate - S * target_op * S_inv) + # jac: d(op_term) = -(dS * target_op * S_inv + S * target_op * -(S_inv * dS * S_inv) ) + # d(op_term) = (-dS * target_op * S_inv + S * target_op * (S_inv * dS * S_inv) ) + + # objective: rho_term = (rho - S * target_rho) + # jac: d(rho_term) = - dS * target_rho + + # objective: ET_term = (E.T - target_E.T * S_inv) + # jac: d(ET_term) = - target_E.T * -(S_inv * dS * S_inv) + # d(ET_term) = target_E.T * (S_inv * dS * S_inv) + + #Distribute computation across processors + allDerivColSlice = slice(0, N) + derivSlices, myDerivColSlice, derivOwners, mySubComm = \ + _mpit.distribute_slice(allDerivColSlice, comm) + if mySubComm is not None: + _warnings.warn("Note: more CPUs(%d)" % comm.Get_size() + + " than gauge-opt derivative columns(%d)!" % N) # pragma: no cover + + n = _slct.length(myDerivColSlice) + wrtIndices = _slct.indices(myDerivColSlice) if (n < N) else None + my_jacMx = jacMx[:, myDerivColSlice] # just the columns I'm responsible for + + # S, and S_inv are shape (d,d) + #S = gauge_group_el.transform_matrix + S_inv = gauge_group_el.transform_matrix_inverse + dS = gauge_group_el.deriv_wrt_params(wrtIndices) # shape (d*d),n + dS.shape = (d, d, n) # call it (d1,d2,n) + dS = _np.rollaxis(dS, 2) # shape (n, d1, d2) + assert(dS.shape == (n, d, d)) + + # --- NOTE: ordering here, with running `start` index MUST + # correspond to those in Model.residuals, which in turn + # must correspond to those in ForwardSimulator.residuals - which + # currently orders as: gates, simplified_ops, preps, effects. + + # -- LinearOperator terms + # ------------------------- + for lbl, G in mdl_pre.operations.items(): + # d(op_term) = S_inv * (-dS * S_inv * G * S + G * dS) = S_inv * (-dS * G' + G * dS) + # Note: (S_inv * G * S) is G' (transformed G) + wt = item_weights.get(lbl, opWeight) + left = -1 * _np.dot(dS, mdl_post.operations[lbl].to_dense('minimal')) # shape (n,d1,d2) + right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # shape (d1,n,d2) -> (n,d1,d2) + result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2) # shape (d1, d2, n) + result = result.reshape((d**2, n)) # must copy b/c non-contiguous + my_jacMx[start:start + d**2] = wt * result + start += d**2 + + # -- Instrument terms + # ------------------------- + for ilbl, Inst in mdl_pre.instruments.items(): + wt = item_weights.get(ilbl, opWeight) + for lbl, G in Inst.items(): + # same calculation as for operation terms + left = -1 * _np.dot(dS, mdl_post.instruments[ilbl][lbl].to_dense('minimal')) # (n,d1,d2) + right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # (d1,n,d2) -> (n,d1,d2) result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2) # shape (d1, d2, n) result = result.reshape((d**2, n)) # must copy b/c non-contiguous my_jacMx[start:start + d**2] = wt * result start += d**2 - # -- Instrument terms - # ------------------------- - for ilbl, Inst in mdl_pre.instruments.items(): - wt = item_weights.get(ilbl, opWeight) - for lbl, G in Inst.items(): - # same calculation as for operation terms - left = -1 * _np.dot(dS, mdl_post.instruments[ilbl][lbl].to_dense('minimal')) # (n,d1,d2) - right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # (d1,n,d2) -> (n,d1,d2) - result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2) # shape (d1, d2, n) - result = result.reshape((d**2, n)) # must copy b/c non-contiguous - my_jacMx[start:start + d**2] = wt * result - start += d**2 - - # -- prep terms - # ------------------------- - for lbl, rho in mdl_post.preps.items(): - # d(rho_term) = -(S_inv * dS * S_inv) * rho - # Note: (S_inv * rho) is transformed rho - wt = item_weights.get(lbl, spamWeight) - Sinv_dS = _np.dot(S_inv, dS) # shape (d1,n,d2) - result = -1 * _np.dot(Sinv_dS, rho.to_dense('minimal')) # shape (d,n) - my_jacMx[start:start + d] = wt * result + # -- prep terms + # ------------------------- + for lbl, rho in mdl_post.preps.items(): + # d(rho_term) = -(S_inv * dS * S_inv) * rho + # Note: (S_inv * rho) is transformed rho + wt = item_weights.get(lbl, spamWeight) + Sinv_dS = _np.dot(S_inv, dS) # shape (d1,n,d2) + result = -1 * _np.dot(Sinv_dS, rho.to_dense('minimal')) # shape (d,n) + my_jacMx[start:start + d] = wt * result + start += d + + # -- effect terms + # ------------------------- + for povmlbl, povm in mdl_pre.povms.items(): + for lbl, E in povm.items(): + # d(ET_term) = E.T * dS + wt = item_weights.get(povmlbl + "_" + lbl, spamWeight) + result = _np.dot(E.to_dense('minimal')[None, :], dS).T # shape (1,n,d2).T => (d2,n,1) + my_jacMx[start:start + d] = wt * result.squeeze(2) # (d2,n) start += d - # -- effect terms - # ------------------------- - for povmlbl, povm in mdl_pre.povms.items(): - for lbl, E in povm.items(): - # d(ET_term) = E.T * dS - wt = item_weights.get(povmlbl + "_" + lbl, spamWeight) - result = _np.dot(E.to_dense('minimal')[None, :], dS).T # shape (1,n,d2).T => (d2,n,1) - my_jacMx[start:start + d] = wt * result.squeeze(2) # (d2,n) - start += d - - # -- penalty terms -- Note: still use original gauge transform applied to `model` - # ------------------------- - if cptp_penalty_factor > 0 or spam_penalty_factor > 0: - if frobenius_transform_target: # reset back to non-target-tranform "mode" - gauge_group_el = original_gauge_group_el - mdl_pre = model.copy() - mdl_post = mdl_pre.copy() - mdl_post.transform_inplace(gauge_group_el) - - if cptp_penalty_factor > 0: - start += _cptp_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, - gauge_group_el, cptp_penalty_factor, - mdl_pre.basis, wrtIndices) - - if spam_penalty_factor > 0: - start += _spam_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, - gauge_group_el, spam_penalty_factor, - mdl_pre.basis, wrtIndices) - - #At this point, each proc has filled the portions (columns) of jacMx that - # it's responsible for, and so now we gather them together. - _mpit.gather_slices(derivSlices, derivOwners, jacMx, [], 1, comm) - #Note jacMx is completely filled (on all procs) - - if check_jac and (comm is None or comm.Get_rank() == 0): - def _mock_objective_fn(v): - return _objective_fn(gauge_group_el, False) - vec = gauge_group_el.to_vector() - _opt.check_jac(_mock_objective_fn, vec, jacMx, tol=1e-5, eps=1e-9, err_type='abs', - verbosity=1) - - return jacMx - - else: - # non-least-squares case where objective function returns a single float - # and (currently) there's no analytic jacobian - - assert gates_metric != "frobeniustt" - assert spam_metric != "frobeniustt" - # ^ PR #410 removed support for Frobenius transform-target metrics in this codepath. - - dim = int(_np.sqrt(mxBasis.dim)) - if n_leak > 0: - B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) - P = B @ B.T.conj() - if _np.linalg.norm(P.imag) > 1e-12: - msg = f"Attempting to run leakage-aware gauge optimization with basis {mxBasis}\n" - msg += "is resulting an orthogonal projector onto the computational subspace that\n" - msg += "is not real-valued. Try again with a different basis, like 'l2p1' or 'gm'." - raise ValueError(msg) - else: - P = P.real - transform_mx_arg = (P, _tools.matrixtools.IdentityOperator()) - # ^ The semantics of this tuple are defined by the frobeniusdist function - # in the ExplicitOpModelCalc class. - else: - transform_mx_arg = None - # ^ It would be equivalent to set this to a pair of IdentityOperator objects. - - def _objective_fn(gauge_group_el, oob_check): - mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) - ret = 0 + # -- penalty terms -- Note: still use original gauge transform applied to `model` + # ------------------------- + if cptp_penalty_factor > 0 or spam_penalty_factor > 0: + if frobenius_transform_target: # reset back to non-target-tranform "mode" + gauge_group_el = original_gauge_group_el + mdl_pre = model.copy() + mdl_post = mdl_pre.copy() + mdl_post.transform_inplace(gauge_group_el) if cptp_penalty_factor > 0: - mdl.basis = mxBasis # set basis for jamiolkowski iso - cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) - ret += _np.sum(cpPenaltyVec) + start += _cptp_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, + gauge_group_el, cptp_penalty_factor, + mdl_pre.basis, wrtIndices) if spam_penalty_factor > 0: - mdl.basis = mxBasis - spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) - ret += _np.sum(spamPenaltyVec) - - if target_model is None: - return ret - - if "frobenius" in gates_metric: - if spam_metric == gates_metric: - val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) - else: - wts = item_weights.copy() - wts['spam'] = 0.0 - for k in wts: - if k in mdl.preps or k in mdl.povms: - wts[k] = 0.0 - val = mdl.frobeniusdist(target_model, transform_mx_arg, wts, n_leak) - if "squared" in gates_metric: - val = val ** 2 - ret += val - - elif gates_metric == "fidelity": - # If n_leak==0, then subspace_entanglement_fidelity is just entanglement_fidelity - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - top = target_model.operations[opLbl].to_dense() - mop = mdl.operations[opLbl].to_dense() - ret += wt * (1.0 - _tools.subspace_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 - - elif gates_metric == "tracedist": - # If n_leak==0, then subspace_jtracedist is just jtracedist. - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - top = target_model.operations[opLbl].to_dense() - mop = mdl.operations[opLbl].to_dense() - ret += wt * _tools.subspace_jtracedist(top, mop, mxBasis, n_leak) + start += _spam_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, + gauge_group_el, spam_penalty_factor, + mdl_pre.basis, wrtIndices) - else: - raise ValueError("Invalid gates_metric: %s" % gates_metric) + #At this point, each proc has filled the portions (columns) of jacMx that + # it's responsible for, and so now we gather them together. + _mpit.gather_slices(derivSlices, derivOwners, jacMx, [], 1, comm) + #Note jacMx is completely filled (on all procs) - if "frobenius" in spam_metric and gates_metric == spam_metric: - # We already handled SPAM error in this case. Just return. - return ret + if check_jac and (comm is None or comm.Get_rank() == 0): + def _mock_objective_fn(v): + return _objective_fn(gauge_group_el, False) + vec = gauge_group_el.to_vector() + _opt.check_jac(_mock_objective_fn, vec, jacMx, tol=1e-5, eps=1e-9, err_type='abs', + verbosity=1) - if "frobenius" in spam_metric: - # SPAM and gates can have different choices for squared vs non-squared. - wts = item_weights.copy(); wts['gates'] = 0.0 - for k in wts: - if k in mdl.operations or k in mdl.instruments: - wts[k] = 0.0 - val = mdl.frobeniusdist(target_model, transform_mx_arg, wts) - if "squared" in spam_metric: - val = val ** 2 - ret += val - - elif spam_metric == "fidelity": - # Leakage-aware metrics NOT available - for preplabel, m_prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) - t_prep = target_model.preps[preplabel] - rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) - ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 - - for povmlabel in mdl.povms.keys(): - wt = item_weights.get(povmlabel, spamWeight) - fidelity = _tools.povm_fidelity(mdl, target_model, povmlabel) - ret += wt * (1.0 - fidelity)**2 - - elif spam_metric == "tracedist": - # Leakage-aware metrics NOT available. - for preplabel, m_prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) - t_prep = target_model.preps[preplabel] - rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) - ret += wt * _tools.tracedist(rhoMx1, rhoMx2) - - for povmlabel in mdl.povms.keys(): - wt = item_weights.get(povmlabel, spamWeight) - ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) - - else: - raise ValueError("Invalid spam_metric: %s" % spam_metric) + return jacMx - return ret + return _objective_fn, _jacobian_fn - _jacobian_fn = None - return _objective_fn, _jacobian_fn +def _create_objective_fn(model, target_model, item_weights: Optional[dict[str,float]]=None, + cptp_penalty_factor: float=0.0, spam_penalty_factor: float=0.0, + gates_metric="frobenius", spam_metric="frobenius", + method=None, comm=None, check_jac=False, n_leak: Optional[int]=None) -> tuple[GGElObjective, GGElJacobian]: + if item_weights is None: + item_weights = dict() + if method == 'ls': + return _legacy_create_least_squares_objective( + model, target_model, item_weights, cptp_penalty_factor, spam_penalty_factor, gates_metric, + spam_metric, comm, check_jac + ) + else: + return _legacy_create_scalar_objective( + model, target_model, item_weights, cptp_penalty_factor, spam_penalty_factor, gates_metric, + spam_metric, n_leak + ) def _cptp_penalty_size(mdl): diff --git a/pygsti/baseobjs/basis.py b/pygsti/baseobjs/basis.py index 9e5734a1e..8d13db56b 100644 --- a/pygsti/baseobjs/basis.py +++ b/pygsti/baseobjs/basis.py @@ -19,6 +19,8 @@ from typing import ( Union, Optional, + Any, + Sequence ) import numpy as _np @@ -263,6 +265,8 @@ def __init__(self, name: str, longname: str, real: bool, sparse: bool): self.longname = longname self.real = real # whether coefficients must be real (*not* whether elements are real - they're always complex) self.sparse = sparse # whether elements are stored as sparse vectors/matrices + self._is_hermitian = None + self._implies_leakage = None @property def dim(self): @@ -317,6 +321,31 @@ def elsize(self): if self.elshape is None: return 0 return int(_np.prod(self.elshape)) + # @property + # def ellookup(self) -> dict[str, Any]: + # raise NotImplementedError() + + # @property + # def elements(self) -> Sequence[Any]: + # raise NotImplementedError() + + # @property.setter + # def elements(self, val) -> None: + # self. + + # @property + # def labels(self) -> Sequence[Any]: + # raise ValueError() + + @property + def implies_leakage_modeling(self) -> bool: + if (not hasattr(self, '_implies_leakage')) or self._implies_leakage is None: + if 'I' not in self.labels: + return False + I = self.ellookup['I'] + self._implies_leakage = _np.linalg.matrix_rank(I) < I.shape[0] + return self._implies_leakage + @property def first_element_is_identity(self): """ @@ -356,6 +385,17 @@ def is_partial(self) -> bool: """ return not self.is_complete() + def is_hermitian(self) -> bool: + from pygsti.tools.matrixtools import is_hermitian as matrix_is_hermitian + if self.elndim != 2 or self.elshape[0] != self.elshape[1]: + return False + if not hasattr(self, '_is_hermitian') or self._is_hermitian is None: + tol = self.dim * _np.finfo(self.elements[0].dtype).eps + is_hermitian = all([matrix_is_hermitian(el, tol=tol) for el in self.elements]) + self._is_hermitian = is_hermitian + return self._is_hermitian + + @property def vector_elements(self): """ diff --git a/pygsti/baseobjs/basisconstructors.py b/pygsti/baseobjs/basisconstructors.py index 5ba4e40bd..f37473770 100644 --- a/pygsti/baseobjs/basisconstructors.py +++ b/pygsti/baseobjs/basisconstructors.py @@ -697,6 +697,7 @@ def gm_labels(matrix_dim): lblList.extend(["Z_{%d}" % (k) for k in range(1, d)]) return lblList + def lf_labels(matrix_dim: int) -> tuple[str,...]: if matrix_dim != 3: raise NotImplementedError() @@ -718,6 +719,7 @@ def lf_labels(matrix_dim: int) -> tuple[str,...]: ) return lbls + def lf_matrices(matrix_dim: int) -> list[_np.ndarray]: """ This basis is used to isolate the parts of Hilbert-Schmidt space that act on @@ -743,7 +745,6 @@ def lf_matrices(matrix_dim: int) -> list[_np.ndarray]: return leakage_basis_mxs - def qsim_matrices(matrix_dim): """ Get the elements of the QuantumSim basis with matrix dimension `matrix_dim`. diff --git a/pygsti/modelmembers/operations/linearop.py b/pygsti/modelmembers/operations/linearop.py index 7a323145c..a7097ba07 100644 --- a/pygsti/modelmembers/operations/linearop.py +++ b/pygsti/modelmembers/operations/linearop.py @@ -14,7 +14,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 optools as _ot, matrixtools as _mt from pygsti import SpaceT from typing import Any @@ -416,11 +416,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): diff --git a/pygsti/modelmembers/povms/effect.py b/pygsti/modelmembers/povms/effect.py index e1e0d46ee..b1aa31546 100644 --- a/pygsti/modelmembers/povms/effect.py +++ b/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): """ diff --git a/pygsti/models/explicitcalc.py b/pygsti/models/explicitcalc.py index 6717e1f66..2c1120a0f 100644 --- a/pygsti/models/explicitcalc.py +++ b/pygsti/models/explicitcalc.py @@ -172,6 +172,7 @@ def frobeniusdist(self, other_calc, transform_mx: Union[None, _np.ndarray, Trans if item_weights is None: item_weights = {} opWeight = item_weights.get('gates', 1.0) spamWeight = item_weights.get('spam', 1.0) + prepWeight = item_weights.get('prep', spamWeight) for opLabel, gate in self.operations.items(): wt = item_weights.get(opLabel, opWeight) @@ -179,7 +180,7 @@ def frobeniusdist(self, other_calc, transform_mx: Union[None, _np.ndarray, Trans nSummands += wt * (gate.dim)**2 for lbl, rhoV in self.preps.items(): - wt = item_weights.get(lbl, spamWeight) + wt = item_weights.get(lbl, prepWeight) d += wt * rhoV.frobeniusdist_squared(other_calc.preps[lbl], P, invP) nSummands += wt * rhoV.dim diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py index 5d276d9e8..7e3ae7c30 100644 --- a/pygsti/models/explicitmodel.py +++ b/pygsti/models/explicitmodel.py @@ -10,6 +10,8 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** +from __future__ import annotations + import collections as _collections import itertools as _itertools import uuid as _uuid @@ -1635,6 +1637,10 @@ def add_availability(opkey, op): availability, instrument_names=list(self.instruments.keys()), nonstd_instruments=self.instruments) + def copy(self) -> ExplicitOpModel: + c = super().copy() # <-- that is indeed an ExplicitOpModel + return c # type: ignore + def create_modelmember_graph(self): return _MMGraph({ 'preps': self.preps, diff --git a/pygsti/models/model.py b/pygsti/models/model.py index 5f6f8e940..27dba1742 100644 --- a/pygsti/models/model.py +++ b/pygsti/models/model.py @@ -10,6 +10,8 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** +from __future__ import annotations + import bisect as _bisect import copy as _copy import itertools as _itertools @@ -2309,7 +2311,7 @@ def _post_copy(self, copy_into, memo): copy_into._reinit_opcaches() super(OpModel, self)._post_copy(copy_into, memo) - def copy(self): + def copy(self) -> OpModel: """ Copy this model. @@ -2320,13 +2322,13 @@ def copy(self): """ self._clean_paramvec() # ensure _paramvec is rebuilt if needed if OpModel._pcheck: self._check_paramvec() - ret = Model.copy(self) + ret = Model.copy(self) # <-- don't be fooled! That's an OpModel. if self._param_bounds is not None and self.parameter_labels is not None: ret._clean_paramvec() # will *always* rebuild paramvec; do now so we can preserve param bounds assert _np.all(self.parameter_labels == ret.parameter_labels) # ensure ordering is the same ret._param_bounds = self._param_bounds.copy() if OpModel._pcheck: ret._check_paramvec() - return ret + return ret # type: ignore def create_modelmember_graph(self): """ diff --git a/pygsti/objectivefns/objectivefns.py b/pygsti/objectivefns/objectivefns.py index a88975f16..d69141b44 100644 --- a/pygsti/objectivefns/objectivefns.py +++ b/pygsti/objectivefns/objectivefns.py @@ -10,6 +10,8 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** +from typing import Optional + import itertools as _itertools import sys as _sys import time as _time @@ -29,6 +31,11 @@ from pygsti.models.model import OpModel as _OpModel from pygsti import SpaceT + +DEFAULT_MIN_PROB_CLIP = 1e-4 +DEFAULT_RADIUS = 1e-4 + + def _objfn(objfn_cls, model, dataset, circuits=None, regularization=None, penalties=None, op_label_aliases=None, comm=None, mem_limit=None, method_names=None, array_types=None, @@ -186,24 +193,25 @@ def create_from(cls, objective='logl', freq_weighted_chi2=False): """ if objective == "chi2": if freq_weighted_chi2: + mpc = RawFreqWeightedChi2Function.DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING builder = FreqWeightedChi2Function.builder( name='fwchi2', description="Freq-weighted sum of Chi^2", - regularization={'min_freq_clip_for_weighting': 1e-4}) + regularization={'min_freq_clip_for_weighting': mpc}) + else: + mpc = RawChi2Function.DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING builder = Chi2Function.builder( name='chi2', description="Sum of Chi^2", - regularization={'min_prob_clip_for_weighting': 1e-4}) + regularization={'min_prob_clip_for_weighting': mpc}) elif objective == "logl": builder = PoissonPicDeltaLogLFunction.builder( name='dlogl', description="2*Delta(log(L))", - regularization={'min_prob_clip': 1e-4, - 'radius': 1e-4}, - penalties={'cptp_penalty_factor': 0, - 'spam_penalty_factor': 0}) + regularization={'min_prob_clip': DEFAULT_MIN_PROB_CLIP, 'radius': DEFAULT_RADIUS}, + penalties={'cptp_penalty_factor': 0, 'spam_penalty_factor': 0}) elif objective == "tvd": builder = TVDFunction.builder( @@ -1660,6 +1668,9 @@ class RawChi2Function(RawObjectiveFunction): verbosity : int, optional Level of detail to print to stdout. """ + + DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING = DEFAULT_MIN_PROB_CLIP + def __init__(self, regularization=None, resource_alloc=None, name="chi2", description="Sum of Chi^2", verbosity=0): super().__init__(regularization, resource_alloc, name, description, verbosity) @@ -1678,7 +1689,7 @@ def chi2k_distributed_qty(self, objective_function_value): """ return objective_function_value - def set_regularization(self, min_prob_clip_for_weighting=1e-4): + def set_regularization(self, min_prob_clip_for_weighting: Optional[float]=None): """ Set regularization values. @@ -1692,7 +1703,11 @@ def set_regularization(self, min_prob_clip_for_weighting=1e-4): ------- None """ - self.min_prob_clip_for_weighting = min_prob_clip_for_weighting + if min_prob_clip_for_weighting is not None: + self.min_prob_clip_for_weighting = min_prob_clip_for_weighting + else: + self.min_prob_clip_for_weighting = RawChi2Function.DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING + return def lsvec(self, probs, counts, total_counts, freqs, intermediates=None): """ @@ -2286,7 +2301,6 @@ def _zero_freq_dterms_relaxed(self, total_counts, probs): # negative lsvec possible class RawFreqWeightedChi2Function(RawChi2Function): - """ The function `N(p-f)^2 / f` @@ -2307,6 +2321,9 @@ class RawFreqWeightedChi2Function(RawChi2Function): verbosity : int, optional Level of detail to print to stdout. """ + + DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING = RawChi2Function.DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING + def __init__(self, regularization=None, resource_alloc=None, name="fwchi2", description="Sum of freq-weighted Chi^2", verbosity=0): super().__init__(regularization, resource_alloc, name, description, verbosity) @@ -2326,7 +2343,7 @@ def chi2k_distributed_qty(self, objective_function_value): """ return objective_function_value # default is to assume the value *is* chi2_k distributed - def set_regularization(self, min_freq_clip_for_weighting=1e-4): + def set_regularization(self, min_freq_clip_for_weighting: Optional[float]=None): """ Set regularization values. @@ -2340,7 +2357,11 @@ def set_regularization(self, min_freq_clip_for_weighting=1e-4): ------- None """ - self.min_freq_clip_for_weighting = min_freq_clip_for_weighting + if min_freq_clip_for_weighting is not None: + self.min_freq_clip_for_weighting = min_freq_clip_for_weighting + else: + self.min_freq_clip_for_weighting = RawFreqWeightedChi2Function.DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING + return def _weights(self, p, f, total_counts): #Note: this could be computed once and cached? @@ -2744,8 +2765,10 @@ def chi2k_distributed_qty(self, objective_function_value): """ return 2 * objective_function_value # 2 * deltaLogL is what is chi2_k distributed - def set_regularization(self, min_prob_clip=1e-4, pfratio_stitchpt=None, pfratio_derivpt=None, - radius=1e-4, fmin=None): + def set_regularization(self, + min_prob_clip=DEFAULT_MIN_PROB_CLIP, pfratio_stitchpt=None, pfratio_derivpt=None, + radius=DEFAULT_RADIUS, fmin=None + ): """ Set regularization values. @@ -3144,7 +3167,7 @@ def chi2k_distributed_qty(self, objective_function_value): """ return 2 * objective_function_value # 2 * deltaLogL is what is chi2_k distributed - def set_regularization(self, min_prob_clip=1e-4, pfratio_stitchpt=None, pfratio_derivpt=None): + def set_regularization(self, min_prob_clip=DEFAULT_MIN_PROB_CLIP, pfratio_stitchpt=None, pfratio_derivpt=None): """ Set regularization values. @@ -5413,7 +5436,7 @@ def _array_types_for_method(cls, method_name, fsim): if method_name == 'dlsvec': return fsim._array_types_for_method('bulk_fill_timedep_dchi2') return super()._array_types_for_method(method_name, fsim) - def set_regularization(self, min_prob_clip_for_weighting=1e-4, radius=1e-4): + def set_regularization(self, min_prob_clip_for_weighting: Optional[float]=None, radius: float = DEFAULT_RADIUS): """ Set regularization values. @@ -5430,6 +5453,9 @@ def set_regularization(self, min_prob_clip_for_weighting=1e-4, radius=1e-4): ------- None """ + if min_prob_clip_for_weighting is None: + mpc = getattr(self.raw_objfn, 'DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING', DEFAULT_MIN_PROB_CLIP) + min_prob_clip_for_weighting = mpc self.min_prob_clip_for_weighting = min_prob_clip_for_weighting self.radius = radius # parameterizes "roundness" of f == 0 terms @@ -5572,7 +5598,7 @@ def _array_types_for_method(cls, method_name, fsim): if method_name == 'dlsvec': return fsim._array_types_for_method('bulk_fill_timedep_dloglpp') return super()._array_types_for_method(method_name, fsim) - def set_regularization(self, min_prob_clip=1e-4, radius=1e-4): + def set_regularization(self, min_prob_clip=DEFAULT_MIN_PROB_CLIP, radius=DEFAULT_RADIUS): """ Set regularization values. @@ -5828,7 +5854,7 @@ def _paramvec_norm_penalty(reg_factor, paramvec): return out -def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis, wrt_slice): +def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis, wrt_slice, error_tol=1e-4): """ Helper function - jacobian of CPTP penalty (sum of tracenorms of gates) Returns a (real) array of shape (len(mdl.operations), n_params). @@ -5842,7 +5868,7 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis #get sgn(chi-matrix) == d(|chi|_Tr)/dchi in std basis # so sgnchi == d(|chi_std|_Tr)/dchi_std chi = _tools.fast_jamiolkowski_iso_std(gate, op_basis) - assert(_np.linalg.norm(chi - chi.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(chi - chi.T.conjugate()) < error_tol), \ "chi should be Hermitian!" # Alt#1 way to compute sgnchi (evals) - works equally well to svd below @@ -5855,7 +5881,7 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis # _spl.sqrtm(_np.matrix(_np.dot(chi.T.conjugate(),chi))))) sgnchi = _tools.matrix_sign(chi) - assert(_np.linalg.norm(sgnchi - sgnchi.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(sgnchi - sgnchi.T.conjugate()) < error_tol), \ "sgnchi should be Hermitian!" # get d(gate)/dp in op_basis [shape == (nP,dim,dim)] @@ -5872,7 +5898,8 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis m_dgate_dp_std = _np.empty((nparams, mdl.dim, mdl.dim), 'complex') for p in range(nparams): # p indexes param m_dgate_dp_std[p] = _tools.fast_jamiolkowski_iso_std(dgate_dp[p], op_basis) # now "M(dGdp_std)" - assert(_np.linalg.norm(m_dgate_dp_std[p] - m_dgate_dp_std[p].T.conjugate()) < 1e-8) # check hermitian + assert(_np.linalg.norm(m_dgate_dp_std[p] - m_dgate_dp_std[p].T.conjugate()) < error_tol**2) # check hermitian + # Riley note: not clear to me why we're using the square of error_tol in the above. m_dgate_dp_std = _np.conjugate(m_dgate_dp_std) # so element-wise multiply # of complex number (einsum below) results in separately adding @@ -5883,7 +5910,7 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis #v = _np.einsum("ij,aij->a",sgnchi,MdGdp_std) v = _np.tensordot(sgnchi, m_dgate_dp_std, ((0, 1), (1, 2))) v *= prefactor * (0.5 / _np.sqrt(_tools.tracenorm(chi))) # add 0.5/|chi|_Tr factor - assert(_np.linalg.norm(v.imag) < 1e-4) + assert(_np.linalg.norm(v.imag) < error_tol) cp_penalty_vec_grad_to_fill[i, :] = 0.0 gate_gpinds_subset, within_wrtslice, within_gpinds = _slct.intersect_within(wrt_slice, gate.gpindices) @@ -5894,7 +5921,7 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis return len(mdl.operations) # the number of leading-dim indicies we filled in -def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_basis, wrt_slice): +def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_basis, wrt_slice, error_tol=1e-4): """ Helper function - jacobian of CPTP penalty (sum of tracenorms of gates) Returns a (real) array of shape ( _spam_penalty_size(mdl), n_params). @@ -5917,11 +5944,11 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_bas #get sgn(denMx) == d(|denMx|_Tr)/d(denMx) in std basis # dmDim = denMx.shape[0] denmx = _tools.vec_to_stdmx(prepvec, op_basis) - assert(_np.linalg.norm(denmx - denmx.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(denmx - denmx.T.conjugate()) < error_tol), \ "denMx should be Hermitian!" sgndm = _tools.matrix_sign(denmx) - assert(_np.linalg.norm(sgndm - sgndm.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(sgndm - sgndm.T.conjugate()) < error_tol), \ "sgndm should be Hermitian!" # get d(prepvec)/dp in op_basis [shape == (nP,dim)] @@ -5936,7 +5963,7 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_bas #v = _np.einsum("ij,aij,ab->b",sgndm,ddenMxdV,dVdp) v = _np.tensordot(_np.tensordot(sgndm, ddenmx_dv, ((0, 1), (1, 2))), dv_dp, (0, 0)) v *= prefactor * (0.5 / _np.sqrt(_tools.tracenorm(denmx))) # add 0.5/|denMx|_Tr factor - assert(_np.linalg.norm(v.imag) < 1e-4) + assert(_np.linalg.norm(v.imag) < error_tol) spam_penalty_vec_grad_to_fill[i, :] = 0.0 gpinds_subset, within_wrtslice, within_gpinds = _slct.intersect_within(wrt_slice, prepvec.gpindices) spam_penalty_vec_grad_to_fill[i, within_wrtslice] = v.real[within_gpinds] # slice or array index works! @@ -5953,11 +5980,11 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_bas #get sgn(EMx) == d(|EMx|_Tr)/d(EMx) in std basis emx = _tools.vec_to_stdmx(effectvec, op_basis) # dmDim = EMx.shape[0] - assert(_np.linalg.norm(emx - emx.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(emx - emx.T.conjugate()) < error_tol), \ "EMx should be Hermitian!" sgn_e = _tools.matrix_sign(emx) - assert(_np.linalg.norm(sgn_e - sgn_e.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(sgn_e - sgn_e.T.conjugate()) < error_tol), \ "sgnE should be Hermitian!" # get d(prepvec)/dp in op_basis [shape == (nP,dim)] @@ -5972,7 +5999,7 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_bas #v = _np.einsum("ij,aij,ab->b",sgnE,dEMxdV,dVdp) v = _np.tensordot(_np.tensordot(sgn_e, demx_dv, ((0, 1), (1, 2))), dv_dp, (0, 0)) v *= prefactor * (0.5 / _np.sqrt(_tools.tracenorm(emx))) # add 0.5/|EMx|_Tr factor - assert(_np.linalg.norm(v.imag) < 1e-4) + assert(_np.linalg.norm(v.imag) < error_tol) spam_penalty_vec_grad_to_fill[i, :] = 0.0 gpinds_subset, within_wrtslice, within_gpinds = _slct.intersect_within(wrt_slice, effectvec.gpindices) diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 83680912f..283ca42b9 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -1511,6 +1511,7 @@ def run(self, data, memlimit=None, comm=None, checkpoint=None, checkpoint_path=N return ret + class LinearGateSetTomography(_proto.Protocol): """ The linear gate set tomography protocol. @@ -1971,13 +1972,6 @@ def _load_pspec(processorspec_filename_or_obj): return processorspec_filename_or_obj -def _load_model(model_filename_or_obj): - if isinstance(model_filename_or_obj, str): - return _io.load_model(model_filename_or_obj) - else: - return model_filename_or_obj # assume a Model object - - def _load_pspec_or_model(processorspec_or_model_filename_or_obj): if isinstance(processorspec_or_model_filename_or_obj, str): # if a filename is given, just try to load a processor spec (can't load a model file yet) @@ -2010,25 +2004,6 @@ def _load_fiducials_and_germs(prep_fiducial_list_or_filename, return prep_fiducials, meas_fiducials, germs -def _load_dataset(data_filename_or_set, comm, verbosity): - """Loads a DataSet from the data_filename_or_set argument of functions in this module.""" - printer = _baseobjs.VerbosityPrinter.create_printer(verbosity, comm) - if isinstance(data_filename_or_set, str): - if comm is None or comm.Get_rank() == 0: - if _os.path.splitext(data_filename_or_set)[1] == ".pkl": - with open(data_filename_or_set, 'rb') as pklfile: - ds = _pickle.load(pklfile) - else: - ds = _io.read_dataset(data_filename_or_set, True, "aggregate", printer) - if comm is not None: comm.bcast(ds, root=0) - else: - ds = comm.bcast(None, root=0) - else: - ds = data_filename_or_set # assume a Dataset object - - return ds - - def _add_gaugeopt_and_badfit(results, estlbl, target_model, gaugeopt_suite, unreliable_ops, badfit_options, optimizer, resource_alloc, printer): tref = _time.time() @@ -2243,19 +2218,16 @@ def _add_badfit_estimates(results, base_estimate_label, badfit_options, #If this estimate is the target model then skip adding the diamond distance wildcard. if base_estimate_label != 'Target': try: - budget = _compute_wildcard_budget_1d_model(base_estimate, objfn_cache, mdc_objfn, parameters, - badfit_options, printer - 1, gaugeopt_suite) - + budget = _compute_wildcard_budget_1d_model(base_estimate, mdc_objfn, printer - 1, gaugeopt_suite) #if gaugeopt_labels is None then budget with be a PrimitiveOpsSingleScaleWildcardBudget #if it was non-empty though then we'll instead have budge as a dictionary with keys #corresponding to the elements of gaugeopt_labels. In this case let's make #base_estimate.extra_parameters['wildcard1d' + "_unmodeled_error"] a dictionary of #the serialized PrimitiveOpsSingleScaleWildcardBudget elements - if gaugeopt_suite is not None and gaugeopt_suite.gaugeopt_suite_names is not None: - gaugeopt_labels = gaugeopt_suite.gaugeopt_suite_names + if gaugeopt_suite is not None: + gaugeopt_labels = budget.keys() base_estimate.extra_parameters['wildcard1d' + "_unmodeled_error"] = {lbl: budget[lbl].to_nice_serialization() for lbl in gaugeopt_labels} - base_estimate.extra_parameters['wildcard1d' + "_unmodeled_active_constraints"] \ - = None + base_estimate.extra_parameters['wildcard1d' + "_unmodeled_active_constraints"] = None base_estimate.extra_parameters["unmodeled_error"] = {lbl: budget[lbl].to_nice_serialization() for lbl in gaugeopt_labels} base_estimate.extra_parameters["unmodeled_active_constraints"] = None @@ -2316,7 +2288,7 @@ def _add_badfit_estimates(results, base_estimate_label, badfit_options, gauge_opt_params.copy(), go_gs_final, gokey, comm, printer - 1) -def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, parameters, badfit_options, verbosity, gaugeopt_suite=None): +def _compute_wildcard_budget_1d_model(estimate, mdc_objfn, verbosity, gaugeopt_suite): """ Create a wildcard budget for a model estimate. This version of the function produces a wildcard estimate using the model introduced in https://doi.org/10.1038/s41534-023-00764-y. @@ -2326,20 +2298,12 @@ def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, paramete estimate : Estimate The estimate object containing the model and data to be used. - objfn_cache : ObjectiveFunctionCache - A cache for objective function evaluations. - mdc_objfn : ModelDatasetCircuitsStore An object that stores the model, dataset, and circuits to be used in the computation. parameters : dict Various parameters of the estimate at hand. - badfit_options : GSTBadFitOptions, optional - Options specifying what post-processing actions should be performed when - a fit is unsatisfactory. Contains detailed parameters for wildcard budget - creation. - verbosity : int, optional Level of detail printed to stdout. @@ -2351,12 +2315,10 @@ def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, paramete Returns ------- PrimitiveOpsWildcardBudget or dict - The computed wildcard budget. If gauge optimization is applied, a dictionary of - budgets keyed by gauge optimization labels is returned. + A dictionary of budgets keyed by gauge optimization labels. """ printer = _baseobjs.VerbosityPrinter.create_printer(verbosity, mdc_objfn.resource_alloc) - badfit_options = GSTBadFitOptions.cast(badfit_options) model = mdc_objfn.model ds = mdc_objfn.dataset global_circuits_to_use = mdc_objfn.global_circuits @@ -2378,46 +2340,54 @@ def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, paramete two_dlogl_threshold = _chi2.ppf(1 - percentile, max(ds_dof - nparams, 1)) redbox_threshold = _chi2.ppf(1 - percentile / nboxes, 1) - ref, reference_name = _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_suite) - - if gaugeopt_suite is None or gaugeopt_suite.gaugeopt_suite_names is None: - gaugeopt_labels = None - primitive_ops = list(ref.keys()) - if sum([v**2 for v in ref.values()]) < 1e-4: - _warnings.warn("Reference values for 1D wildcard budget are all near-zero!" - " This usually indicates an incorrect target model and will likely cause problems computing alpha.") - - else: - gaugeopt_labels = gaugeopt_suite.gaugeopt_suite_names - primitive_ops = list(ref[list(gaugeopt_labels)[0]].keys()) - if sum([v**2 for v in ref[list(gaugeopt_labels)[0]].values()]) < 1e-4: - _warnings.warn("Reference values for 1D wildcard budget are all near-zero!" - " This usually indicates an incorrect target model and will likely cause problems computing alpha.") - - if gaugeopt_labels is None: - wcm = _wild.PrimitiveOpsSingleScaleWildcardBudget(primitive_ops, [ref[k] for k in primitive_ops], - reference_name=reference_name) - _opt.wildcardopt.optimize_wildcard_bisect_alpha(wcm, mdc_objfn, two_dlogl_threshold, redbox_threshold, printer, + target = gaugeopt_suite.gaugeopt_target + if target is None: + target = estimate.models['target'] + + # Getting the dict of gauge-optimized models is very annoying. + if gaugeopt_suite.gaugeopt_argument_dicts is None: + gaugeopt_suite.gaugeopt_argument_dicts = dict() + # ^ We make that assignment because _compute_1d_reference_values_and_name + # shouldn't have to include reassign. + gop_dicts = gaugeopt_suite.gaugeopt_argument_dicts + gop_names = gaugeopt_suite.gaugeopt_suite_names + if gop_names is None: + gop_names = [gn for gn in gop_dicts if gn in estimate.models] + + goped_models = {lbl: estimate.models[lbl] for lbl in gop_names} + if len(goped_models) == 0: + final_iter_model = estimate.models['final iteration estimate'] + if isinstance(final_iter_model, _ExplicitOpModel): + goped_model = _alg.gaugeopt_to_target(final_iter_model, target) + goped_models['simple gaugeopt'] = goped_model + else: + goped_models['no gaugeopt'] = final_iter_model + + from pygsti.baseobjs.label import Label + ref : dict[str,dict[Label, float]] = _compute_1d_reference_values_and_name(target, goped_models, gaugeopt_suite) + reference_name = 'diamond distance' + + gaugeopt_labels = list(ref.keys()) + primitive_ops = list(ref[gaugeopt_labels[0]].keys()) + if sum([v**2 for v in ref[gaugeopt_labels[0]].values()]) < 1e-4: + _warnings.warn("Reference values for 1D wildcard budget are all near-zero!" + " This usually indicates an incorrect target model and will likely cause problems computing alpha.") + + wcm = {lbl:_wild.PrimitiveOpsSingleScaleWildcardBudget(primitive_ops, [ref[lbl][k] for k in primitive_ops], + reference_name=reference_name) for lbl in gaugeopt_labels} + for budget_to_optimize in wcm.values(): + _opt.wildcardopt.optimize_wildcard_bisect_alpha(budget_to_optimize, mdc_objfn, two_dlogl_threshold, redbox_threshold, printer, guess=0.1, tol=1e-3) # results in optimized wcm - else: - wcm = {lbl:_wild.PrimitiveOpsSingleScaleWildcardBudget(primitive_ops, [ref[lbl][k] for k in primitive_ops], - reference_name=reference_name) for lbl in gaugeopt_labels} - for budget_to_optimize in wcm.values(): - _opt.wildcardopt.optimize_wildcard_bisect_alpha(budget_to_optimize, mdc_objfn, two_dlogl_threshold, redbox_threshold, printer, - guess=0.1, tol=1e-3) # results in optimized wcm return wcm -def _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_suite=None): +def _compute_1d_reference_values_and_name(target_model, gopped_models, gaugeopt_suite): """ - Compute the reference values and name for the 1D wildcard budget. + Compute the reference values the 1D wildcard budget. Parameters ---------- - estimate : Estimate - The estimate object containing the models to be used. - badfit_options : GSTBadFitOptions Options specifying what post-processing actions should be performed when a fit is unsatisfactory. Contains detailed parameters for wildcard budget @@ -2429,101 +2399,77 @@ def _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_sui Returns ------- - dict or dict of dicts - The computed reference values for the 1D wildcard budget. If gauge optimization is applied, - a dictionary of reference values keyed by gauge optimization labels is returned. - - str - The name of the reference metric used. + dict of dicts + A dictionary of 1D wildcard budget reference values, keyed by gauge optimization labels. """ - if badfit_options.wildcard1d_reference == 'diamond distance': - if gaugeopt_suite is None or gaugeopt_suite.gaugeopt_suite_names is None: - final_model = estimate.models['final iteration estimate'] - target_model = estimate.models['target'] - - if isinstance(final_model, _ExplicitOpModel): - gaugeopt_model = _alg.gaugeopt_to_target(final_model, target_model) - operations_dict = gaugeopt_model.operations - targetops_dict = target_model.operations - preps_dict = gaugeopt_model.preps - targetpreps_dict = target_model.preps - povmops_dict = gaugeopt_model.povms - insts_dict = gaugeopt_model.instruments - targetinsts_dict = target_model.instruments - - else: - # Local/cloud noise models don't have default_gauge_group attribute and can't be gauge - # optimized - at least not easily. - gaugeopt_model = final_model - operations_dict = gaugeopt_model.operation_blks['gates'] - targetops_dict = target_model.operation_blks['gates'] - preps_dict = gaugeopt_model.prep_blks['layers'] - targetpreps_dict = target_model.prep_blks['layers'] - povmops_dict = {} # HACK - need to rewrite povm_diamonddist below to work - insts_dict = gaugeopt_model.instrument_blks['layers'] - targetinsts_dict = target_model.instrument_blks['layers'] - - dd = {} - for key, op in operations_dict.items(): - dd[key] = 0.5 * _tools.diamonddist(op.to_dense(), targetops_dict[key].to_dense()) - if dd[key] < 0: # indicates that diamonddist failed (cvxpy failure) - _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" - " Falling back to trace distance.") % str(key)) - dd[key] = _tools.jtracedist(op.to_dense(), targetops_dict[key].to_dense()) - - for key, op in insts_dict.items(): - inst_dd = 0.5* _tools.instrument_diamonddist(op, targetinsts_dict[key]) - if inst_dd < 0: # indicates that instrument_diamonddist failed - _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" - "No fallback presently available for instruments, so skipping.") % str(key)) - else: - dd[key] = inst_dd - - spamdd = {} - for key, op in preps_dict.items(): - spamdd[key] = _tools.tracedist(_tools.vec_to_stdmx(op.to_dense(), 'pp'), - _tools.vec_to_stdmx(targetpreps_dict[key].to_dense(), 'pp')) + assert isinstance(gopped_models, dict) + assert all([isinstance(k, str) for k in gopped_models]) + + def _memberdicts(_m) -> tuple[dict,dict,dict,dict]: + if isinstance(_m, _ExplicitOpModel): + ops = _m.operations + preps = _m.preps + povms = _m.povms + insts = _m.instruments + else: + ops = _m.operation_blks['gates'] + preps = _m.prep_blks['layers'] + povms = {} # disabled until povm_diamonddist can handle non-ExplicitOpModels. + insts = _m.instrument_blks['layers'] + return ops, preps, povms, insts + + target_ops, target_preps, target_povms, target_insts = _memberdicts(target_model) + + dd = {lbl: dict() for lbl in gopped_models} + for lbl, gaugeopt_model in gopped_models.items(): + + argdicts = gaugeopt_suite.gaugeopt_argument_dicts.get(lbl, dict()) + n_leak: Optional[int] = None + if isinstance(argdicts, list) and len(argdicts) > 0: + n_leak = argdicts[0].get('n_leak', n_leak) + + basis = gaugeopt_model.basis + udim = int(_np.round(_np.sqrt(basis.dim))) + I = _tools.matrixtools.IdentityOperator() + if n_leak is None and basis.first_element_is_identity: + P = I + elif n_leak is None: + P = _tools.superop_subspace_projector(basis) + elif n_leak > 0: + P = _tools.superop_subspace_projector(udim - n_leak, udim, basis) + else: + P = I - for key in povmops_dict.keys(): - spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) + ops, preps, _, insts = _memberdicts(gaugeopt_model) + mx_basis = gaugeopt_model.basis - dd['SPAM'] = sum(spamdd.values()) + for key, op in ops.items(): + dd[lbl][key] : float = 0.5 * _tools.diamonddist(op.to_dense() @ P, target_ops[key].to_dense() @ P, mx_basis=mx_basis) # type: ignore + if dd[lbl][key] < 0: # indicates that diamonddist failed (cvxpy failure) + _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" + " Falling back to trace distance.") % str(key)) + dd[lbl][key] = _tools.jtracedist(op.to_dense(), target_ops[key].to_dense()) + + for key, op in insts.items(): + inst_dd : float = 0.5* _tools.instrument_diamonddist(op, target_insts[key], mx_basis=mx_basis) # type: ignore + if inst_dd < 0: # indicates that instrument_diamonddist failed + _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" + "No fallback presently available for instruments, so skipping.") % str(key)) + else: + dd[lbl][key] = inst_dd - else: - #pull the target model from the gaugeopt suite. The GSTGaugeOptSuite doesn't currently - #support multiple targets, so this is safe for now. - target_model = gaugeopt_suite.gaugeopt_target - gaugeopt_models = [estimate.models[lbl] for lbl in gaugeopt_suite.gaugeopt_suite_names] - dd = {lbl: {} for lbl in gaugeopt_suite.gaugeopt_suite_names} - for gaugeopt_model, lbl in zip(gaugeopt_models, gaugeopt_suite.gaugeopt_suite_names): - for key, op in gaugeopt_model.operations.items(): - dd[lbl][key] = 0.5 * _tools.diamonddist(op.to_dense(), target_model.operations[key].to_dense()) - if dd[lbl][key] < 0: # indicates that diamonddist failed (cvxpy failure) - _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" - " Falling back to trace distance.") % str(key)) - dd[lbl][key] = _tools.jtracedist(op.to_dense(), target_model.operations[key].to_dense()) - - for key, op in gaugeopt_model.instruments.items(): - inst_dd = 0.5* _tools.instrument_diamonddist(op, target_model.instruments[key], gaugeopt_model.basis) - if inst_dd < 0: # indicates that instrument_diamonddist failed - _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" - "No fallback presently available for instruments, so skipping.") % str(key)) - else: - dd[lbl][key] = inst_dd - spamdd = {} - for key, op in gaugeopt_model.preps.items(): - spamdd[key] = _tools.tracedist(_tools.vec_to_stdmx(op.to_dense(), 'pp'), - _tools.vec_to_stdmx(target_model.preps[key].to_dense(), 'pp')) + spamdd = {} + for key, op in preps.items(): + rho1 = _tools.vec_to_stdmx(op.to_dense(), basis=mx_basis) + rho2 = _tools.vec_to_stdmx(target_preps[key].to_dense(), basis=mx_basis) + spamdd[key] = _tools.tracedist(rho1, rho2) - for key in gaugeopt_model.povms.keys(): - spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) + for key in target_povms: + spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) # type: ignore - dd[lbl]['SPAM'] = sum(spamdd.values()) + dd[lbl]['SPAM'] = sum(spamdd.values()) - return dd, 'diamond distance' - else: - raise ValueError("Invalid wildcard1d_reference value (%s) in bad-fit options!" - % str(badfit_options.wildcard1d_reference)) + return dd def _compute_robust_scaling(scale_typ, objfn_cache, mdc_objfn): @@ -3301,6 +3247,11 @@ def _add_param_preserving_gauge_opt(results: ModelEstimateResults, est_key: str, ggel = gop_dictorlist['_gaugeGroupEl'] model_implicit_gauge = transform_composed_model(est.models['final iteration estimate'], ggel) est.models[gop_name] = model_implicit_gauge + protocol = est.parameters['protocol'] + badfit_options = protocol.badfit_options + optimizer = protocol.optimizer + resource_alloc = _baseobjs.resourceallocation.ResourceAllocation() + _add_badfit_estimates(results, est_key, badfit_options, optimizer, resource_alloc, verbosity, gaugeopt_suite=gop_params) return diff --git a/pygsti/report/factory.py b/pygsti/report/factory.py index c8781f53f..aae5827bd 100644 --- a/pygsti/report/factory.py +++ b/pygsti/report/factory.py @@ -1230,7 +1230,7 @@ def construct_standard_report(results, title="auto", ws = ws or _ws.Workspace() advanced_options = advanced_options or {} - n_leak = advanced_options.get('n_leak', 0) + n_leak = advanced_options.get('n_leak', None) # ^ It would be preferable to store n_leak in a Basis object, or something similar. # We're using this for now since it's simple and gets the job done. linlogPercentile = advanced_options.get('linlog percentile', 5) diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index f92dc8d36..5111059ea 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -17,6 +17,7 @@ """ import importlib import warnings as _warnings +from typing import Union import numpy as _np import scipy.linalg as _spl @@ -40,6 +41,8 @@ FINITE_DIFF_EPS = 1e-7 +BasisLike = Union[str, _Basis] + def _null_fn(*arg): return None @@ -1016,14 +1019,23 @@ def maximum_trace_dist(gate, mx_basis): Maximum_trace_dist = _modf.opfn_factory(maximum_trace_dist) # init args == (model, op_label) + def leaky_maximum_trace_dist(gate, mx_basis): closestUOpMx = _alg.find_closest_unitary_opmx(gate) _tools.jamiolkowski_iso(closestUOpMx, mx_basis, mx_basis) - n_leak = 1 + # .... ^ What does that function call do ??? + if not isinstance(mx_basis, _Basis): + mx_basis = _Basis.cast(mx_basis, dim=gate.shape[0]) + if mx_basis.implies_leakage_modeling: + n_leak = None # we'll infer a detailed leakage model from mx_basis + else: + n_leak = 1 # assume one leakage level return _tools.subspace_jtracedist(gate, closestUOpMx, mx_basis, n_leak) + Leaky_maximum_trace_dist = _modf.opfn_factory(leaky_maximum_trace_dist) + def diamonddist_to_leakfree_cptp(op, ignore, mx_basis): import pygsti.tools.sdptools as _sdps prob, _, solvers = _sdps.diamond_distance_projection_model( @@ -1037,8 +1049,10 @@ def diamonddist_to_leakfree_cptp(op, ignore, mx_basis): continue return -1 + Diamonddist_to_leakfree_cptp = _modf.opsfn_factory(diamonddist_to_leakfree_cptp) + def subspace_diamonddist_to_leakfree_cptp(op, ignore, mx_basis): import pygsti.tools.sdptools as _sdps prob, _, solvers = _sdps.diamond_distance_projection_model( @@ -1052,49 +1066,41 @@ def subspace_diamonddist_to_leakfree_cptp(op, ignore, mx_basis): continue return -1 + SubspaceDiamonddist_to_leakfree_cptp = _modf.opsfn_factory(subspace_diamonddist_to_leakfree_cptp) + def subspace_diamonddist(op_a, op_b, basis): - dim_mixed = op_a.shape[0] - dim_pure = int(dim_mixed**0.5) - dim_pure_compsub = dim_pure - 1 - from pygsti.tools.leakage import leading_dxd_submatrix_basis_vectors - U = leading_dxd_submatrix_basis_vectors(dim_pure_compsub, dim_pure, basis) - P = U @ U.T.conj() - assert _np.linalg.norm(P - P.real) < 1e-10 - P = P.real - from pygsti.tools.optools import diamonddist - return diamonddist(op_a @ P, op_b @ P, basis) / 2 + from pygsti.tools.leakage import subspace_diamonddist as _subspace_diamonddist + return _subspace_diamonddist(op_a, op_b, basis) + SubspaceDiamonddist = _modf.opsfn_factory(subspace_diamonddist) -def pergate_leakrate_reduction(op, ignore, mx_basis, reduction): - assert op.shape == (9, 9) - lfb = _BuiltinBasis('l2p1', 9) - op_lfb = _tools.change_basis(op, mx_basis, lfb) - elinds = lfb.elindlookup - compinds = [elinds[sslbl] for sslbl in ['I','X','Y','Z'] ] - leakage_effect_superket = op_lfb[elinds['L'], compinds] - leakage_effect = _tools.vec_to_stdmx(leakage_effect_superket, 'pp') - leakage_rates = _np.linalg.eigvalsh(leakage_effect) - return reduction(leakage_rates) def pergate_leakrate_max(op, ignore, mx_basis): - return pergate_leakrate_reduction(op, ignore, mx_basis, max) + from pygsti.tools.leakage import gate_leakage_profile + evals, _ = gate_leakage_profile(op, mx_basis) + return evals[0] -def pergate_leakrate_min(op, ignore, mx_basis): - return pergate_leakrate_reduction(op, ignore, mx_basis, min) PerGateLeakRateMax = _modf.opsfn_factory(pergate_leakrate_max) + + +def pergate_leakrate_min(op, ignore, mx_basis): + from pygsti.tools.leakage import gate_leakage_profile + evals, _ = gate_leakage_profile(op, mx_basis) + return evals[-1] + + PerGateLeakRateMin = _modf.opsfn_factory(pergate_leakrate_min) + def pergate_seeprate(op, ignore, mx_basis): - assert op.shape == (9, 9) - lfb = _BuiltinBasis('l2p1', 9) - op_lfb = _tools.change_basis(op, mx_basis, lfb) - elinds = lfb.elindlookup - seeprate = op_lfb[elinds['I'], elinds['L']] - return seeprate + from pygsti.tools.leakage import gate_seepage_profile + seeprate, _ = gate_seepage_profile(op, mx_basis) + return seeprate[0] + PerGateSeepRate = _modf.opsfn_factory(pergate_seeprate) @@ -1169,10 +1175,17 @@ def entanglement_fidelity(a, b, mx_basis): Entanglement_fidelity = _modf.opsfn_factory(entanglement_fidelity) # init args == (model1, model2, op_label) + def subspace_entanglement_fidelity(a, b, mx_basis): - n_leak = 1 + if not isinstance(mx_basis, _Basis): + mx_basis = _Basis.cast(mx_basis, dim=a.shape[0]) + if mx_basis.implies_leakage_modeling: + n_leak = None # we'll infer a detailed leakage model from mx_basis + else: + n_leak = 1 # assume one leakage level return _tools.subspace_entanglement_fidelity(a, b, mx_basis, n_leak) + Subspace_entanglement_fidelity = _modf.opsfn_factory(subspace_entanglement_fidelity) @@ -1201,9 +1214,11 @@ def entanglement_infidelity(a, b, mx_basis): Entanglement_infidelity = _modf.opsfn_factory(entanglement_infidelity) # init args == (model1, model2, op_label) + def leaky_entanglement_infidelity(a, b, mx_basis): return 1 - subspace_entanglement_fidelity(a, b, mx_basis) + Leaky_entanglement_infidelity = _modf.opsfn_factory(leaky_entanglement_infidelity) @@ -1273,7 +1288,12 @@ def frobenius_diff(a, b, mx_basis): # assume vary model1, model2 fixed def leaky_gate_frob_dist(a, b, mx_basis): - n_leak = 1 + if not isinstance(mx_basis, _Basis): + mx_basis = _Basis.cast(mx_basis, dim=a.shape[0]) + if mx_basis.implies_leakage_modeling: + n_leak = None # we'll infer a detailed leakage model from mx_basis + else: + n_leak = 1 # assume one leakage level return _tools.subspace_superop_fro_dist(a, b, mx_basis, n_leak) @@ -1305,10 +1325,17 @@ def jtrace_diff(a, b, mx_basis): # assume vary model1, model2 fixed Jt_diff = _modf.opsfn_factory(jtrace_diff) # init args == (model1, model2, op_label) + def leaky_jtrace_diff(a, b, mx_basis): - n_leak = 1 + if not isinstance(mx_basis, _Basis): + mx_basis = _Basis.cast(mx_basis, dim=a.shape[0]) + if mx_basis.implies_leakage_modeling: + n_leak = None # we'll infer a detailed leakage model from mx_basis + else: + n_leak = 1 # assume one leakage level return _tools.subspace_jtracedist(a, b, mx_basis, n_leak) + Leaky_Jt_diff = _modf.opsfn_factory(leaky_jtrace_diff) @@ -1577,32 +1604,74 @@ def eigenvalue_nonunitary_avg_gate_infidelity(a, b, mx_basis): # init args == (model1, model2, op_label) -def eigenvalue_entanglement_infidelity(a, b, mx_basis): +def eigenvalue_entanglement_infidelity(a: _np.ndarray, b: _np.ndarray, mx_basis: BasisLike, is_tp=None, is_unitary=None, tol=1e-8) -> _np.floating: """ - Eigenvalue entanglement infidelity between a and b + Eigenvalue entanglement infidelity is the infidelity between certain diagonal matrices that + contain [see Note 1] the eigenvalues of J(a) and J(b), where J(.) is the Jamiolkowski + isomorphism map. This is equivalent to computing infidelity of (J(a), J(b)) while pretending + that they share an eigenbasis. Parameters ---------- - a : numpy.ndarray - The first process (transfer) matrix. + is_tp : bool, optional (default None) + Flag indicating that a and b are TP in their provided basis. If None (the default), + an explicit check is performed up to numerical tolerance `tol`. If True/False, then + the check is skipped and this is used as the result of the check as though it were + performed. - b : numpy.ndarray - The second process (transfer) matrix. + is_unitary : bool, optional (default None) + Flag indicating that b is unitary. If None (the default) an explicit check is performed + up to tolerance `tol`. If True/False, then the check is skipped and this is used as the + result of the check as though it were performed. - mx_basis : Basis or {'pp', 'gm', 'std'} - the basis that `a` and `b` are in. + Notes + ----- + [1] The eigenvalues are ordered using a min-weight matching algorithm. - Returns - ------- - float + [2] If a and b are trace-preserving (TP) and b is unitary, then we compute eigenvalue + entanglement fidelity by leveraging three facts: + + 1. If (x,y) share an eigenbasis and have (consistently ordered) eigenvalues (v(x),v(y)), + then Tr(x y^{\\dagger}) is the inner product of v(x) and v(y). + + 2. J(a) and J(b) share an eigenbasis of if `a` and `b` share an eigenbasis. + + 3. If `a` and `b` are TP and `b` is unitary, then their entanglement fidelity can be + expressed as |Tr( a.b^{\\dagger} )| / d^2. + + Then together, we see that if (a,b) are TP and `b` is unitary, then their eigenvalue + entanglement fidelity is expressible as the inner product of v(a) and v(b), where + v(a) and v(b) vectors holding are suitably-ordered eigenvalues of a and b. + """ d2 = a.shape[0] - evA = _np.linalg.eigvals(a) - evB = _np.linalg.eigvals(b) - _, pairs = _tools.minweight_match(evA, evB, lambda x, y: abs(x - y), - return_pairs=True) # just to get pairing - mlPl = abs(_np.sum([_np.conjugate(evB[j]) * evA[i] for i, j in pairs])) - return 1.0 - mlPl / float(d2) + + if is_unitary is None: + is_unitary = _np.allclose(_np.eye(d2), b @ b.T.conj(), atol=tol, rtol=tol) + if is_tp is None: + is_tp = _tools.is_trace_preserving(a, mx_basis, tol) and _tools.is_trace_preserving(b, mx_basis, tol) + + if is_unitary and is_tp: + evA = _np.linalg.eigvals(a) + evB = _np.linalg.eigvals(b) + _, pairs = _tools.minweight_match(evA, evB, lambda x, y: abs(x - y), + return_pairs=True) # just to get pairing + fid = abs(_np.sum([_np.conjugate(evB[j]) * evA[i] for i, j in pairs])) / d2 + else: + Ja = _tools.fast_jamiolkowski_iso_std(a, mx_basis) + Jb = _tools.fast_jamiolkowski_iso_std(b, mx_basis) + from pygsti.tools.optools import eigenvalue_fidelity + fid = eigenvalue_fidelity(Ja, Jb, gauge_invariant=True) + if fid < 0.9: + pass + # valsA, vecsA = _spl.eigh(Ja) + # valsB, vecsB = _spl.eigh(Jb) + # dissimilarity = lambda vec_x, vec_y : abs(1 - vec_x @ vec_y) + # _, pairs = _tools.minweight_match(vecsA.T.conj(), vecsB.T.conj(), dissimilarity, return_pairs=True) + # fid = abs(_np.sum([_np.conjugate(evB[j]) * evA[i] for i, j in pairs])) / d2 + # fid = _np.linalg.norm(_np.sqrt(evA) * _np.sqrt(evB), ord=1)**2 + + return 1.0 - fid Eigenvalue_entanglement_infidelity = _modf.opsfn_factory(eigenvalue_entanglement_infidelity) @@ -1643,7 +1712,8 @@ def eigenvalue_avg_gate_infidelity(a, b, mx_basis): def eigenvalue_diamondnorm(a, b, mx_basis): """ - Eigenvalue diamond distance between a and b + Eigenvalue diamond distance between `a` and `b`, assuming + that `a` and `b` are TP and that `b` is unitary. Parameters ---------- @@ -1663,8 +1733,8 @@ def eigenvalue_diamondnorm(a, b, mx_basis): d2 = a.shape[0] evA = _np.linalg.eigvals(a) evB = _np.linalg.eigvals(b) - return (d2 - 1.0) / d2 * _np.max(_tools.minweight_match(evA, evB, lambda x, y: abs(x - y), - return_pairs=False)) + temp = _tools.minweight_match(evA, evB, lambda x, y: abs(x - y), return_pairs=False) + return (d2 - 1.0) / d2 * _np.max(temp) Eigenvalue_diamondnorm = _modf.opsfn_factory(eigenvalue_diamondnorm) diff --git a/pygsti/report/section/gauge.py b/pygsti/report/section/gauge.py index 219bba41f..8a92ecf8f 100644 --- a/pygsti/report/section/gauge.py +++ b/pygsti/report/section/gauge.py @@ -167,7 +167,17 @@ def final_model_spam_vs_target_table(workspace, switchboard=None, confidence_lev @_Section.figure_factory(4) def final_gates_vs_target_table_gauge_var(workspace, switchboard=None, confidence_level=None, ci_brevity=1, **kwargs): - if kwargs.get('n_leak', 0) == 0: + + if switchboard is not None and 'mdl_target_grid' in switchboard: + basis_representative = switchboard['mdl_target_grid'][0].basis + else: + res_representative = kwargs['results'][list(kwargs['results'])[0]] + est_representative = res_representative.estimates[list(res_representative.estimates)[0]] + mdl_representative = est_representative.models[list(est_representative.models)[-1]] + basis_representative = mdl_representative.basis + n_leak_default = None if basis_representative.implies_leakage_modeling else 0 + + if kwargs.get('n_leak', n_leak_default) == 0: display = ('inf', 'agi', 'geni', 'trace', 'diamond', 'nuinf', 'nuagi') else: display = ('sub-inf', 'sub-trace', 'sub-diamond', 'plf-sub-diamond', 'leak-rate-max', 'leak-rate-min', 'seep-rate' ) diff --git a/pygsti/report/section/summary.py b/pygsti/report/section/summary.py index 4d6764763..b841db073 100644 --- a/pygsti/report/section/summary.py +++ b/pygsti/report/section/summary.py @@ -38,7 +38,17 @@ def final_model_fit_histogram(workspace, switchboard=None, linlog_percentile=5, @_Section.figure_factory() def final_gates_vs_target_table_insummary(workspace, switchboard=None, confidence_level=None, ci_brevity=1, show_unmodeled_error=False, **kwargs): - if kwargs.get('n_leak', 0) == 0: + + if switchboard is not None and 'mdl_target_grid' in switchboard: + basis_representative = switchboard['mdl_target_grid'][0].basis + else: + res_representative = kwargs['results'][list(kwargs['results'])[0]] + est_representative = res_representative.estimates[list(res_representative.estimates)[0]] + mdl_representative = est_representative.models[list(est_representative.models)[-1]] + basis_representative = mdl_representative.basis + n_leak_default = None if basis_representative.implies_leakage_modeling else 0 + + if kwargs.get('n_leak', n_leak_default) == 0: summary_display = ('inf', 'trace', 'diamond', 'evinf', 'evdiamond') else: summary_display = ('sub-inf', 'sub-trace', 'sub-diamond', 'plf-sub-diamond', 'leak-rate-max') diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py index f748f6a2b..5c4b68908 100644 --- a/pygsti/tools/leakage.py +++ b/pygsti/tools/leakage.py @@ -14,9 +14,12 @@ import copy from pygsti.tools import optools as pgot from pygsti.tools import basistools as pgbt +from pygsti.tools import matrixtools as pgmt from pygsti.tools.basistools import stdmx_to_vec +from pygsti.baseobjs import Label from pygsti.baseobjs.basis import TensorProdBasis, Basis, BuiltinBasis import numpy as np +import scipy.linalg as la import warnings from typing import Union, Dict, Optional, List, Any, TYPE_CHECKING @@ -26,20 +29,24 @@ from pygsti.processors import QubitProcessorSpec +BasisLike = Union[Basis, str] + -# Question: include the parenthetical in the heading below? NOTATION = \ """ -Default notation (deferential to text above) --------------------------------------------- +Default notation (possibly superceded by text above) +---------------------------------------------------- * H is a complex Hilbert space equipped with the standard basis. * C, the computational subspace, is the complex-linear span of the first dim(C) - standard basis vectors of H. + standard basis vectors of H. The orthogonal complement of C in H is H \\ C. * Given a complex Hilbert space, U, we write M[U] to denote the space of linear operators from U to U. Elements of M[U] have natural matrix representations. + We can interpret A ∈ M[C] as an element of M[H] by requiring that ker(A) and + ker(A^{†}) contain H \\ C. + * Given a space of linear operators, L, we write S[L] for the set of linear transformations ("superoperators") from L to L. @@ -74,15 +81,28 @@ def assign(fn): *A remark.* The Choi matrix for a superoperator G in S[H] is G(rho_t), where rho_t is the output of this function when n_leak=0. """ + NOTATION) -def tensorized_teststate_density(dim: int, n_leak: int) -> np.ndarray: - temp = np.eye(dim, dtype=np.complex128) - if n_leak > 0: - temp[-n_leak:,-n_leak:] = 0.0 - temp /= np.sqrt(dim - n_leak) - psi = pgbt.stdmx_to_stdvec(temp).ravel() - # +def tensorized_teststate_density(*args) -> np.ndarray: + if len(args) == 2: + dim : int = args[0] # type: ignore + n_leak : int = args[1] # type: ignore + args_ok = isinstance(dim, int) and isinstance(n_leak, int) and n_leak < dim + if not args_ok: + raise ValueError() + E = np.eye(dim, dtype=np.complex128) + if n_leak > 0: + E[-n_leak:,-n_leak:] = 0.0 + elif len(args) == 1: + basis : Basis = args[0] # type: ignore + args_ok = isinstance(basis, Basis) and basis.is_hermitian() and ('I' in basis.labels) + if not args_ok: + raise ValueError() + E = basis.ellookup['I'].copy() + else: + raise ValueError() + psi = pgbt.stdmx_to_stdvec(E).ravel() + psi /= la.norm(psi) + # In a standard leakage basis we have # |psi> = (|00> + |11> + ... + |dim - n_leak - 1>) / sqrt(dim - n_leak). - # rho_test = np.outer(psi, psi) return rho_test @@ -102,18 +122,20 @@ def tensorized_teststate_density(dim: int, n_leak: int) -> np.ndarray: *Warning!* At present, this function can only be used for gates over a single system (e.g., a single qubit), not for tensor products of such systems. """ + NOTATION) -def apply_tensorized_to_teststate(op_x: np.ndarray, op_y, op_basis: np.ndarray, n_leak: int=0) -> tuple[TensorProdBasis, np.ndarray, np.ndarray]: - dim = int(np.sqrt(op_x.shape[0])) - assert op_x.shape == (dim**2, dim**2) - assert op_y.shape == (dim**2, dim**2) +def apply_tensorized_to_teststate(op_x: np.ndarray, op_y: np.ndarray, op_basis: BasisLike, n_leak: Optional[int]=None) -> tuple[TensorProdBasis, np.ndarray, np.ndarray]: + udim = int(np.sqrt(op_x.shape[0])) + dim = udim**2 + assert op_x.shape == (dim, dim) + assert op_y.shape == (dim, dim) # op_x and op_y act on M[H]. # # We care about op_x and op_y only up to their action on the subspace - # U = {rho in M[H] : = 0 for all i >= dim - n_leak }. + # U = {rho in M[H] : = 0 for all i >= udim - n_leak }. # # It's easier to talk about this subspace (and related subspaces) if op_x and op_y are in # the standard basis. So the first thing we do is convert to that basis. - std_basis = BuiltinBasis('std', dim**2) + op_basis = Basis.cast(op_basis, dim=dim) + std_basis = BuiltinBasis('std', dim) op_x_std = pgbt.change_basis(op_x, op_basis, std_basis) op_y_std = pgbt.change_basis(op_y, op_basis, std_basis) @@ -125,12 +147,15 @@ def apply_tensorized_to_teststate(op_x: np.ndarray, op_y, op_basis: np.ndarray, # # for all sigma, rho in M[H]. The way we do this implicitly fixes a basis for S[H]⨂S[H] as # the tensor product basis. We'll make that explicit later on. - idle_gate = np.eye(dim**2, dtype=np.complex128) + idle_gate = np.eye(dim, dtype=np.complex128) lift_op_x_std = np.kron(op_x_std, idle_gate) lift_op_y_std = np.kron(op_y_std, idle_gate) # Now we'll compare these lifted operators by how they act on specific state in M[H]⨂M[H]. - rho_test = tensorized_teststate_density(dim, n_leak) + if n_leak is not None: + rho_test = tensorized_teststate_density(udim, n_leak) + else: + rho_test = tensorized_teststate_density(op_basis) # type: ignore # lift_op_x and lift_op_y only act on states in their superket representations, so we convert # rho_test to a superket representation in the induced tensor product basis for S[H]⨂S[H]. @@ -144,71 +169,25 @@ def apply_tensorized_to_teststate(op_x: np.ndarray, op_y, op_basis: np.ndarray, @lru_cache -def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis): - """ - Let "H" denote n^2 dimensional Hilbert-Schdmit space, and let "U" denote the d^2 - dimensional subspace of H spanned by vectors whose Hermitian matrix representations - are zero outside the leading d-by-d submatrix. - - This function returns a column-unitary matrix "B" where P = B B^{\\dagger} is the - orthogonal projector from H to U with respect to current_basis. We return B rather - than P only because it's simpler to get P from B than it is to get B from P. - - See below for this function's original use-case. - - Raison d'etre - ------------- - Suppose we canonically measure the distance between two process matrices (M1, M2) by - - D(M1, M2; H) = max || (M1 - M2) v || - v is in H, (Eq. 1) - tr(v) = 1, - v is positive - - for some norm || * ||. Suppose also that we want an analog of this distance when - (M1, M2) are restricted to the linear subspace U consisting of all vectors in H - whose matrix representations are zero outside of their leading d-by-d submatrix. - - One natural way to do this is via the function D(M1, M2; U) -- i.e., just replace - H in (Eq. 1) with the subspace U. Using P to denote the orthogonal projector onto U, - we claim that we can evaluate this function via the identity - - D(M1, M2; U) = D(M1 P, M2 P; H). (Eq. 2) - - To see why this is the case, consider a positive vector v and its projection u = P v. - Since a vector is called positive whenever its Hermitian matrix representation is positive - semidefinite (PSD), we need to show that u is positive. This can be seen by considering - block 2-by-2 partitions of the matrix representations of (u,v), where the leading block - is d-by-d: +@set_docstring( +""" +Here, H has dimension n and C ⊂ H has dimension d ≤ n. - mat(v) = [x11, x12] and mat(u) = [x11, 0] - [x21, x22] [ 0, 0]. - - In particular, u is positive if and only if x11 is PSD, and x11 must be PSD for v - to be positive. Furthermore, positivity of v requires that x22 is PSD, which implies +This function returns a column-unitary matrix B where P = B B^{\\dagger} is the +orthogonal projector from M[H] to M[C] with respect to current_basis. - 0 <= tr(u) = tr(x11) <= tr(v). - - Given this, it is easy to establish (Eq 2.) by considering how the following pair - of problems have the same optimal objective function value - - max || (M1 - M2) P v || and max || (M1 - M2) P v || - mat(v) = [x11, x12] mat(v) = [x11, x12] - [x21, x22] [x21, x22] - mat(v) is PSD x11 is PSD - tr(x11) + tr(x22) = 1 tr(x11) <= 1. - - In fact, this can be taken a little further! The whole argument goes through unchanged - if, instead of starting with the objective function || (M1 - M2) v ||, we started with - f((M1 - M2) v) and f satisfied the property that f(c v) >= f(v) whenever c is a scalar - greater than or equal to one. +If you only care about P, then you can call superop_subspace_projector instead. +""" + NOTATION) +def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis) -> np.ndarray: + """ + TODO: deprecate this function. """ assert d <= n + if d == n: + return np.eye(n**2) current_basis = Basis.cast(current_basis, dim=n**2) X = current_basis.create_transform_matrix('std') X = X.T.conj() - if d == n: - return X # we have to select a proper subset of columns in current_basis std_basis = BuiltinBasis(name='std', dim_or_statespace=n**2) label2ind = std_basis.elindlookup @@ -222,7 +201,90 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis): return submatrix_basis_vectors -CHOI_INDUCED_METRIC = \ +# TODO: document me +def projective_measurement_subspace_basis_vectors(E: np.ndarray, basis: Basis) -> np.ndarray: + """ + Our desired subspace is + { rho : ker(rho) contains ker(E) } . + Since E is (proportional to) a projector, this is equivalent to + { rho : E rho E = rho }. + """ + if not basis.is_hermitian(): + raise ValueError() + E = E.copy() + k = np.linalg.matrix_rank(E) + E *= (k/np.trace(E)) + if not pgmt.is_projector(E): + raise ValueError() + proj_elements = [ E @ B @ E for B in basis.elements ] + subspace_frame = np.column_stack([ pgbt.stdmx_to_vec(pB, basis) for pB in proj_elements ]) + # ^ a "frame" is an overcomplete basis. + subspace_frame = subspace_frame.real + # ^ Since basis.elements are Hermitian and E is Hermitian, we have that proj_elements are + # Hermitian and that their vectorizations are real. We cast to real here just to + # eliminate possible rounding errors. + U_full : np.ndarray = la.qr(subspace_frame, pivoting=True)[0] # type: ignore + U = U_full[:, :k**2] + return U + + +@lru_cache +@set_docstring( +NOTATION) +def superop_subspace_projector(*args) -> np.ndarray: + + num_positional_args = len(args) + if num_positional_args not in {3, 1}: + raise ValueError() + + if num_positional_args == 1: + basis = args[0] + if not isinstance(basis, Basis): + raise ValueError() + if not hasattr(basis, 'ellookup') or ('I' not in basis.ellookup): + raise ValueError() + dim = basis.dim + if basis.first_element_is_identity: + return np.eye(dim) + E = basis.ellookup['I'] + U = projective_measurement_subspace_basis_vectors(E, basis) + + else: + """ + This function returns the superoperator in S[H] that projects orthogonally from + M[H] to M[C], where H is n-dimensional and C ⊂ H is d-dimensional (d ≤ n). + + The action of this operator is easy to understand when M[H] and M[C] are viewed + as spaces of n-by-n matrices rather than spaces of length-n^2 vectors. + + For v in M[H] and u = P v, we have + + mat(v) = [x11, x12] and mat(u) = [x11, 0] + [x21, x22] [ 0, 0]. + + This characterization makes two facts about P apparent. First, P is positive + (i.e., it takes Hermitian psd operators to Hermitian psd operators). Second, + P is trace-non-increasing. + """ + d, n, basis = args + types_ok = isinstance(d, int) and isinstance(n, int) and isinstance(basis, Basis) + if not types_ok: + raise ValueError() + if d > n: + raise ValueError() + if n**2 != basis.elsize: + raise ValueError() + if d == n: + return np.eye(n**2) + E = np.eye(n) + E[d:, d:] = 0.0 + U = projective_measurement_subspace_basis_vectors(E, basis) + + P = U @ U.T + return P + + +CHOI_INDUCED_METRIC_TEMPLATE = \ """ The pair (op_x, op_y) represent some superoperators (X, Y) in S[H], using op_basis. @@ -236,49 +298,247 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis): """ + NOTATION -@set_docstring(CHOI_INDUCED_METRIC % 'entanglement fidelity') -def subspace_entanglement_fidelity(op_x, op_y, op_basis, n_leak=0): +@set_docstring(CHOI_INDUCED_METRIC_TEMPLATE % 'entanglement fidelity') +def subspace_entanglement_fidelity(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak: Optional[int] = None) -> float: ten_std_basis, temp1, temp2 = apply_tensorized_to_teststate(op_x, op_y, op_basis, n_leak) temp1_mx = pgbt.vec_to_stdmx(temp1, ten_std_basis, keep_complex=True) temp2_mx = pgbt.vec_to_stdmx(temp2, ten_std_basis, keep_complex=True) ent_fid = pgot.fidelity(temp1_mx, temp2_mx) - return ent_fid + return ent_fid # type: ignore -@set_docstring(CHOI_INDUCED_METRIC % 'jamiolkowski trace distance') -def subspace_jtracedist(op_x, op_y, op_basis, n_leak=0): +@set_docstring(CHOI_INDUCED_METRIC_TEMPLATE % 'jamiolkowski trace distance') +def subspace_jtracedist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak: Optional[int] = None) -> float: ten_std_basis, temp1, temp2 = apply_tensorized_to_teststate(op_x, op_y, op_basis, n_leak) temp1_mx = pgbt.vec_to_stdmx(temp1, ten_std_basis, keep_complex=True) temp2_mx = pgbt.vec_to_stdmx(temp2, ten_std_basis, keep_complex=True) j_dist = pgot.tracedist(temp1_mx, temp2_mx) - return j_dist + return j_dist # type: ignore -@set_docstring( + +PROJECTION_INDUCED_METRIC_TEMPLATE = \ """ The pair (op_x, op_y) represent some superoperators (X, Y) in S[H], using op_basis. -We return the Frobenius distance between op_x @ P and op_y @ P, where P is the -projector onto the computational subspace (i.e., C), of co-dimension n_leak. +We return the %s between op_x @ P and op_y @ P, where P is the +projector onto the computational subspace (i.e., C) of co-dimension n_leak. *Warning!* At present, this function can only be used for gates over a single system (e.g., a single qubit), not for tensor products of such systems. -""" + NOTATION) -def subspace_superop_fro_dist(op_x, op_y, op_basis, n_leak=0): +""" + NOTATION + + +@set_docstring(PROJECTION_INDUCED_METRIC_TEMPLATE % 'Frobenius distance') +def subspace_superop_fro_dist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak:Optional[int]=None) -> float: diff = op_x - op_y if n_leak == 0: - return np.linalg.norm(diff, 'fro') + return la.norm(diff, 'fro') # type: ignore + n = int(np.sqrt(op_x.shape[0])) + assert op_x.shape == op_y.shape == (n**2, n**2) + if n_leak is None: + P = superop_subspace_projector(op_basis) + else: + P = superop_subspace_projector(n - n_leak, n, op_basis) + return la.norm(diff @ P) # type: ignore + + +@set_docstring(PROJECTION_INDUCED_METRIC_TEMPLATE % 'diamond distance') +def subspace_diamonddist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak: Optional[int]=None) -> float: + """ + Here we give a brief motivating derivation for defining the subspace diamond norm in + the way that we have. This derivation won't convince a skeptic that our definition + is the best-possible. + + Suppose we canonically measure the distance between two superoperators (X, Y) by + + D(X, Y; H) = max || (X - Y) v || + v is in M[H], (Eq. 1) + tr(v) = 1, + v is positive + + for some norm || * ||. - d = int(np.sqrt(op_x.shape[0])) - assert op_x.shape == op_y.shape == (d**2, d**2) - B = leading_dxd_submatrix_basis_vectors(d-n_leak, d, op_basis) - P = B @ B.T.conj() - return np.linalg.norm(diff @ P) + We arrive at a natural analog of this metric when (X, Y) are restricted to M[C] + simply by replacing "H" in (Eq. 1) with "C". + + Using P to denote the orthogonal projector onto M[C], we claim that + + D(X, Y; C) = D(X P, Y P; H). (Eq. 2) + + Here's a proof of that claim: + + | It's easy to show that P is a positive trace-non-increasing map. In particular, + | if u = P v, then the matrix representations of u and v are + | + | mat(v) = [v11, v12] and mat(u) = [v11, 0] + | [v21, v22] [ 0, 0], + | + | where v11 and v22 are psd if v is positive. From here the claim follows once + | you've convinced yourself that the pair of problems below have the same optimal + | objective value + | + | max || (X - Y) P v || and max || (X - Y) P v || + | mat(v) = [v11, v12] mat(v) = [v11, v12] + | [v21, v22] [v21, v22] + | mat(v) is PSD v11 is PSD + | tr(v11) + tr(v22) = 1 tr(v11) <= 1. + + This can be taken a little further. The proof's argument goes through unchanged if, + instead of starting with the objective || (X - Y) v ||, we started with f((X - Y) v), + where f satisfies the property that f(c v) >= f(v) whenever c is a scalar >= 1. + """ + from pygsti.tools.optools import diamonddist + dim_mixed = op_x.shape[0] + dim_pure = int(dim_mixed**0.5) + if n_leak is None: + P = superop_subspace_projector(op_basis) + else: + P = superop_subspace_projector(dim_pure - n_leak, dim_pure, op_basis) + val : float = diamonddist(op_x @ P, op_y @ P, op_basis, return_x=False) / 2 # type: ignore + return val + + +TRANSPORT_PROFILES_TEMPLATE = \ +""" +This function returns a description of how a gate (in some process matrix representation) +transports population from SUB to the orthogonal complement of SUB in H. + +Underlying mathematics +---------------------- +Our subspace of interest is represented by the unique Hermitian operator E_SUB satisfying +⟨ E_SUB, ρ ⟩ = 1 and ρ = E_SUB ρ E_SUB for all densities ρ ∈ M[SUB]. %s + +The process matrix `op` is the mx_basis representation of a CPTP map G : M[H] ➞ M[H]. + +Consider an experimental protocol where a system is prepared in state ρ ∈ M[SUB], evolved +to G(ρ), and then measured with the 2-element POVM {E_SUB, 1 - E_SUB}. Runs of this protocol +that result in the `1 - E_SUB` measurement outcome are called transport events. + +The population transport profile of (G,SUB) is a full specification of the probabilties +of transport events, considering all possible choices of ρ ∈ M[SUB]. This profile is +encoded in the Hermitian operator + + E_transport := (E_SUB) G^{†}(1 - E_SUB) (E_SUB) ∈ M[SUB] ⊂ M[H]. + +This is a valid representation, since + + Pr{ G transports ρ | ρ ∈ M[SUB] } + = ⟨ 1 - E_SUB , G(ρ) ⟩ // definition of a transport event + = ⟨ G^{†}(1 - E_SUB) , ρ ⟩ // definition of the adjoint + = ⟨ E_transport , ρ ⟩. // using ρ = E_SUB ρ E_SUB + + +Function output and interpretation +---------------------------------- +This function returns a tuple with eigenvalues and eigenvectors of E_transport. + +Say this tuple is `(rates, states)`. If G is CPTP, then + + rates[-1] ≤ Pr{ G transports ρ | ρ ∈ M[SUB] } ≤ rates[0]. + +The upper bound is G's maximum transport of population (Max TOP) out of SUB, and +is acheived by applying G to the pure state `states[0]` ∈ SUB. + +The lower bound is G's minimum transport of population (Min TOP) out of SUB (assuming +population starts in SUB) and is acheived by applying G to `states[-1]` ∈ SUB. + +The length of `rates` and `states` is equal to the rank of E_SUB. + +""" + + +@set_docstring( +TRANSPORT_PROFILES_TEMPLATE.replace('SUB', 'sub') % '' + """ +Failure modes +------------- +We raise a ValueError if pygsti.tools.is_projector(E_sub, E_sub_tol) returns False. +""" + NOTATION ) +def pop_transport_profile(E_sub: np.ndarray, op: np.ndarray, mx_basis: Basis, E_sub_tol=1e-14) -> tuple[np.ndarray,list[np.ndarray]]: + n = int(np.sqrt(E_sub.size)) + if not pgmt.is_projector(E_sub, E_sub_tol): + msg = \ + f""" + The argument E_sub must be an orthogonal projector. The provided value was + + E_sub = {E_sub}; + + this failed the tests in pygsti.tools.is_projector at tolerance={E_sub_tol}. + """ + raise ValueError(msg) + + E_sub_perp_mat = np.eye(n) - E_sub + E_sub_perp_vec = pgbt.stdmx_to_vec(E_sub_perp_mat, mx_basis) + transport_E_vec = op.T @ E_sub_perp_vec + transport_E_mat = pgbt.vec_to_stdmx(transport_E_vec, mx_basis, keep_complex=True) + transport_E_mat = E_sub @ transport_E_mat @ E_sub + + rates, states = la.eigh(transport_E_mat) + dim_proj = int(np.round(np.trace(E_sub).real)) + ind = np.argsort(np.abs(rates))[::-1][:dim_proj] + rates = rates[ind] + states = [s for s in states.T[ind]] + + return rates, states + + +@set_docstring( +TRANSPORT_PROFILES_TEMPLATE.replace('SUB', 'C') % """This operator is +proportional to the element of mx_basis labeled 'I', since that basis element is assumed +be proportional to the orthogonal projector onto C. +""" + NOTATION ) +def gate_leakage_profile(op: np.ndarray, mx_basis: Basis) -> tuple[np.ndarray, list[np.ndarray]]: + mx_basis = Basis.cast(mx_basis, dim=int(np.sqrt(op.size))) + E_comp_mat = mx_basis.ellookup['I'] # type: ignore + n = int(np.sqrt(E_comp_mat.size)) + assert E_comp_mat.shape == (n, n) + dim_comp = np.linalg.matrix_rank(E_comp_mat, hermitian=True) + E_comp_mat = E_comp_mat * (dim_comp / np.trace(E_comp_mat)) + if dim_comp == n: + msg = \ + """ + The provided basis suggests that the computational subspace is equal to + the entire system Hilbert space. Returning with an empty leakage profile! + """ + warnings.warn(msg) + return np.empty((0,)), [] + + rates, states = pop_transport_profile(E_comp_mat, op, mx_basis) + return rates, states + + +@set_docstring( +TRANSPORT_PROFILES_TEMPLATE.replace('SUB', '[H \\ C]') % """This operator is +inferred from the element of mx_basis labeled 'I', since that basis element is assumed +be proportional to the orthogonal projector onto C. +""" + NOTATION ) +def gate_seepage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: + mx_basis = Basis.cast(mx_basis, dim=int(np.sqrt(op.size))) + E_comp_mat = mx_basis.ellookup['I'] # type: ignore + n = int(np.sqrt(E_comp_mat.size)) + assert E_comp_mat.shape == (n, n) + dim_comp = np.linalg.matrix_rank(E_comp_mat, hermitian=True) + E_comp_mat = E_comp_mat * (dim_comp / np.trace(E_comp_mat)) + E_leak_mat = np.eye(n) - E_comp_mat + if dim_comp == n: + msg = \ + """ + The provided basis suggests that the computational subspace is equal to + the entire system Hilbert space. Returning with an empty seepage profile! + """ + warnings.warn(msg) + return np.empty((0,)), [] + + rates, states = pop_transport_profile(E_leak_mat, op, mx_basis) + return rates, states # MARK: model construction -def leaky_qubit_model_from_pspec(ps_2level: QubitProcessorSpec, mx_basis: Union[str, Basis]='l2p1', levels_readout_zero=(0,)) -> ExplicitOpModel: +def leaky_qubit_model_from_pspec( + ps_2level: QubitProcessorSpec, mx_basis: Union[str, Basis]='l2p1', + levels_readout_zero=(0,), default_idle_gatename: Label = Label(()) + ) -> ExplicitOpModel: """ Return an ExplicitOpModel `m` whose (ideal) gates act on three-dimensional Hilbert space and whose members are represented in `mx_basis`, constructed as follows: @@ -322,10 +582,14 @@ def leaky_qubit_model_from_pspec(ps_2level: QubitProcessorSpec, mx_basis: Union[ from pygsti.modelmembers.povms import UnconstrainedPOVM from pygsti.modelmembers.states import FullState assert ps_2level.num_qubits == 1 + if ps_2level.idle_gate_names == ['{idle}']: + ps_2level.rename_gate_inplace('{idle}', default_idle_gatename) if isinstance(mx_basis, str): mx_basis = BuiltinBasis(mx_basis, 9) assert isinstance(mx_basis, Basis) + + ql = ps_2level.qubit_labels[0] Us = ps_2level.gate_unitaries rho0 = np.array([[1, 0, 0], [0, 0, 0], [0, 0, 0]], complex) @@ -333,7 +597,7 @@ def leaky_qubit_model_from_pspec(ps_2level: QubitProcessorSpec, mx_basis: Union[ E0[levels_readout_zero, levels_readout_zero] = 1 E1 = np.eye(3, dtype=complex) - E0 - ss = ExplicitStateSpace([0],[3]) + ss = ExplicitStateSpace([ql],[3]) tm_3level = ExplicitOpModel(ss, mx_basis) # type: ignore tm_3level.preps['rho0'] = FullState(stdmx_to_vec(rho0, mx_basis)) tm_3level.povms['Mdefault'] = UnconstrainedPOVM( @@ -348,9 +612,16 @@ def u2x2_to_9x9_superoperator(u2x2): return superop for gatename, unitary in Us.items(): - gatekey = (gatename, 0) if gatename != '{idle}' else ('Gi',0) + gatekey = gatename if isinstance(gatename, Label) else Label((gatename, ql)) tm_3level.operations[gatekey] = u2x2_to_9x9_superoperator(unitary) + from pygsti.models.gaugegroup import UnitaryGaugeGroup, DirectSumUnitaryGroup, TrivialGaugeGroup, _ExplicitStateSpace + ss_comp = _ExplicitStateSpace(ps_2level.qubit_labels, [2]) + ss_leak = _ExplicitStateSpace(['L'], [1]) + g_comp = UnitaryGaugeGroup(ss_comp, 'pp') + g_leak = TrivialGaugeGroup(ss_leak) + g_full = DirectSumUnitaryGroup((g_comp, g_leak), mx_basis) + tm_3level.default_gauge_group = g_full tm_3level.sim = 'map' # can use 'matrix', if that's preferred for whatever reason. return tm_3level @@ -359,14 +630,58 @@ def u2x2_to_9x9_superoperator(u2x2): def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: """ - Return a modified deep-copy of gopparams_dicts, where the modifications - strip out gauge optimization over the TPSpam gauge group and apply a - few other changes appropriate for leakage modeling. + goppparams_dicts is a list-of-dicts (LoDs) representation of a gauge optimization suite + suitable for models without leakage (e.g., a model of a 2-level system). + + This function returns a new gauge optimization suite (also in the LoDs representation) + by applying leakage-specific modifications to a deep-copy of gopparams_dicts. + + Example + ------- + Suppose we have a ModelEstimateResults object called `results` that includes a + CPTPLND estimate, and we want to update the models of that estimate to include + two types of leakage-aware gauge optimization. + + # + # Step 1: get the input to this function + # + estimate = results.estimates['CPTPLND'] + model = estimate.models['target'] + stdsuite = GSTGaugeOptSuite(gaugeopt_suite_names=('stdgaugeopt',)) + gopparams_dicts = stdsuite.to_dictionary(model)['stdgaugeopt'] + + # + # Step 2: use this function to build our GSTGaugeOptSuite. + # + s = lagoified_gopparams_dicts(gopparams_dicts) + t = lagoified_gopparams_dicts(gopparams_dicts) + t[-1]['gates_metric'] = 'fidelity' + t[-1]['spam_metric'] = 'fidelity' + specification = {'LAGO-std': s,'LAGO-custom': t} + gos = GSTGaugeOptSuite(gaugeopt_argument_dicts=specification) + + # + # Step 3: updating `estimate` requires that we modify `results`. + # + add_lago_models(results, 'CPTPLND', gos) + + After those lines execute, the `estimate.models` dict will have two new + key-value pairs, where the keys are 'LAGO-std' and 'LAGO-custom'. """ + from pygsti.models.gaugegroup import UnitaryGaugeGroup + tm = gopparams_dicts[0]['target_model'] + if tm.basis.implies_leakage_modeling: + # The basis carries enough information for us to infer how leakage should + # be modeled. + n_leak_arg = None + else: + # The basis doesn't say anything about how the leakage should be modeled. + # We'll assume a 3-level system with one leakage level. + n_leak_arg = 1 gopparams_dicts = [gp for gp in gopparams_dicts if 'TPSpam' not in str(type(gp['_gaugeGroupEl']))] gopparams_dicts = copy.deepcopy(gopparams_dicts) for inner_dict in gopparams_dicts: - inner_dict['n_leak'] = 1 + inner_dict['n_leak'] = n_leak_arg # ^ This function could accept n_leak as an argument instead. However, # downstream functions for gauge optimization only support n_leak=0 or 1. # @@ -374,30 +689,45 @@ def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: # about mismatches between an estimate and a target when restricted to the # computational subspace. We have code for evaluating the loss functions # themselves, but not their gradients. + inner_dict['gates_metric'] = 'fidelity' + inner_dict['spam_metric'] = 'fidelity' + inner_dict['item_weights'] = {'gates': 0.0, 'spam': 1.0} + gg = UnitaryGaugeGroup(tm.basis.state_space, tm.basis) + inner_dict['gauge_group'] = gg + inner_dict['_gaugeGroupEl'] = gg.compute_element(gg.initial_params) + # ^ We start with gauge optimization over the full unitary group, minimizing + # SPAM differences between the estimate and the target on the computational + # subspace. Our last step of gauge optimization (which is after this loop) + # includes gates. inner_dict['method'] = 'L-BFGS-B' # ^ We need this optimizer because it doesn't require a gradient oracle. - inner_dict['gates_metric'] = 'frobenius squared' - inner_dict['spam_metric'] = 'frobenius squared' - # ^ We have other metrics for leakage-aware gauge optimization, but they - # aren't really useful. inner_dict['convert_model_to']['to_type'] = 'full' # ^ The natural basis for Hilbert-Schmidt space in leakage modeling doesn't # have the identity matrix as its first element. This means we can't use # the full TP parameterization. There's no real harm in taking "full" as # our default because add_lago_models uses parameterization-preserving # gauge optimization. + inner_dict = inner_dict.copy() + gg = tm.default_gauge_group + # ^ The most likely scenario is that gg is a DirectSumGaugeGroup, consisting + # of a UnitaryGaugeGroup and a TrivialGaugeGroup. Rather than hard-code + # that choice, we go with the default gauge group of the target model. + if gg is not None: + inner_dict['gauge_group'] = gg + inner_dict['_gaugeGroupEl'] = gg.compute_element(gg.initial_params) + inner_dict['gates_metric'] = 'frobenius squared' + inner_dict['spam_metric'] = 'frobenius squared' + inner_dict['n_leak'] = n_leak_arg + inner_dict['item_weights'] = {'gates': 1.0, 'spam': 1.0} + gopparams_dicts.append(inner_dict) return gopparams_dicts -def std_lago_gopsuite(model: ExplicitOpModel) -> dict[str, list[dict]]: +def std_lago_gaugeopt_params(model: ExplicitOpModel) -> dict[str, list[dict]]: """ - Return a dictionary of the form {'LAGO': v}, where v is a - list-of-dicts representation of a gauge optimization suite. - - We construct v by getting the list-of-dicts representation of the - "stdgaugeopt" suite for `model`, and then changing some of its - options to be suitable for leakage-aware gauge optimization. These - changes are made in the `lagoified_gopparams_dicts` function. + Return a dictionary of the form {'LAGO': v}, where v is a list-of-dicts + representation of a gauge optimization suite suitable for leakage modeling, + obtained by modiftying the 'stdgaugeopt' suite induced by `model`. """ from pygsti.protocols.gst import GSTGaugeOptSuite std_gop_suite = GSTGaugeOptSuite(gaugeopt_suite_names=('stdgaugeopt',)) @@ -411,20 +741,23 @@ def add_lago_models(results: ModelEstimateResults, est_key: Optional[str] = None """ Update each estimate in results.estimates (or just results.estimates[est_key], if est_key is not None) with a model obtained by parameterization-preserving - gauge optimization. - - If no gauge optimization suite is provided, then we construct one by starting - with "stdgaugeopt" suite and then changing the settings to be suitable for leakage-aware gauge optimization. + + If no gauge optimization suite is provided, then we construct one by making + appropriate modifications to either the estimate's existing 'stdgaugeopt' suite + (if that exists) or to the 'stdgaugeopt' suite induced by the target model. """ from pygsti.protocols.gst import GSTGaugeOptSuite, _add_param_preserving_gauge_opt - if gos is None: - existing_est = results.estimates[est_key] - std_gos_lods = existing_est.goparameters['stdgaugeopt'] - lago_gos_lods = lagoified_gopparams_dicts(std_gos_lods) - gop_params = {'LAGO': lago_gos_lods} - gos = GSTGaugeOptSuite(gaugeopt_argument_dicts=gop_params) if isinstance(est_key, str): + if gos is None: + existing_est = results.estimates[est_key] + if 'stdgaugeopt' in existing_est.goparameters: + std_gos_lods = existing_est.goparameters['stdgaugeopt'] + lago_gos_lods = lagoified_gopparams_dicts(std_gos_lods) + gop_params = {'LAGO': lago_gos_lods} + else: + gop_params = std_lago_gaugeopt_params(results.estimates[est_key].models['target']) + gos = GSTGaugeOptSuite(gaugeopt_argument_dicts=gop_params) _add_param_preserving_gauge_opt(results, est_key, gos, verbosity) elif est_key is None: for est_key in results.estimates.keys(): @@ -473,10 +806,11 @@ def construct_leakage_report( for ek in est_key: assert isinstance(ek, str) add_lago_models(results, ek, verbosity=gaugeopt_verbosity) + est = results.estimates[ek] + n_leak_arg = None if est.models['LAGO'].basis.implies_leakage_modeling else 1 + from pygsti.report import construct_standard_report report = construct_standard_report( - results, advanced_options={'n_leak': 1}, **extra_report_kwargs + results, advanced_options={'n_leak': n_leak_arg}, **extra_report_kwargs ) return report, results - -print() diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py index 210ef23cd..04728cba2 100644 --- a/pygsti/tools/matrixtools.py +++ b/pygsti/tools/matrixtools.py @@ -65,7 +65,7 @@ def gram_matrix(m, adjoint=False): return out -def is_hermitian(mx, tol=1e-9): +def is_hermitian(mx: _np.ndarray, tol: float = 1e-9) -> bool: """ Test whether mx is a hermitian matrix. @@ -86,7 +86,17 @@ def is_hermitian(mx, tol=1e-9): if m != n: return False else: - return _np.all(_np.abs(mx - mx.T.conj()) <= tol) + return _np.all(_np.abs(mx - mx.T.conj()) <= tol) # type: ignore + + +def assert_hermitian(mat : _np.ndarray, tol: _np.floating) -> None: + hermiticity_error = _np.abs(mat - mat.T.conj()) + if _np.any(hermiticity_error > tol): + message = f""" + Input matrix 'mat' is not Hermitian, up to tolerance {tol}. + The absolute values of entries in (mat - mat^H) are \n{hermiticity_error}. + """ + raise ValueError(message) def is_pos_def(mx, tol=1e-9, attempt_cholesky=False): @@ -139,6 +149,18 @@ def is_valid_density_mx(mx, tol=1e-9): return abs(_np.trace(mx) - 1.0) < tol and is_pos_def(mx, tol) +def is_projector(mx, tol=1e-14): + # We use a stringent default tolerance since projectors tend to be + # computed to high accuracy. + if tol == _np.inf: + return True + if not is_hermitian(mx, tol): + return False + mx2 = mx @ mx + tol = _np.sqrt(mx.size) * tol + return _np.allclose(mx, mx2, atol=tol, rtol=tol) + + def nullspace(m, tol=1e-7): """ Compute the nullspace of a matrix. @@ -533,7 +555,6 @@ def matrix_sign(m): return sign - def print_mx(mx, width=9, prec=4, withbrackets=False): """ Print matrix in pretty format. diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 5dab09e8a..24c23ebac 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -86,18 +86,9 @@ def fidelity(a, b): # ^ use for checks that have no dimensional dependence; about 1e-12 for double precision. __VECTOR_TOL__ = (a.shape[0] ** 0.5) * __SCALAR_TOL__ # ^ use for checks that do have dimensional dependence (will naturally increase for larger matrices) - - def assert_hermitian(mat): - hermiticity_error = _np.abs(mat - mat.T.conj()) - if _np.any(hermiticity_error > __SCALAR_TOL__): - message = f""" - Input matrix 'mat' is not Hermitian, up to tolerance {__SCALAR_TOL__}. - The absolute values of entries in (mat - mat^H) are \n{hermiticity_error}. - """ - raise ValueError(message) - assert_hermitian(a) - assert_hermitian(b) + _mt.assert_hermitian(a, __SCALAR_TOL__) + _mt.assert_hermitian(b, __VECTOR_TOL__) def check_rank_one_density(mat): """ @@ -184,6 +175,32 @@ def psd_square_root(mat): return f +def eigenvalue_fidelity(a, b, gauge_invariant=True) -> _np.floating: + from pygsti.tools import minweight_match + tol = _np.finfo(a.dtype).eps ** 0.75 + _mt.assert_hermitian(a, tol) + _mt.assert_hermitian(b, tol) + if gauge_invariant: + valsA = _spl.eigvalsh(a) + valsB = _spl.eigvalsh(b) + dissimilarity = lambda x, y: abs(x - y) + _, pairs = minweight_match(valsA, valsB, dissimilarity, return_pairs=True) + else: + valsA, vecsA = _spl.eigh(a) + valsB, vecsB = _spl.eigh(b) + dissimilarity = lambda vec_x, vec_y : abs(1 - vec_x @ vec_y) + _, pairs = minweight_match(vecsA.T.conj(), vecsB.T.conj(), dissimilarity, return_pairs=True) + ind_a, ind_b = zip(*pairs) + arg_a = _np.maximum(valsA[list(ind_a)], 0) + arg_b = _np.maximum(valsB[list(ind_b)], 0) + f = _np.linalg.norm(arg_a**0.5 * arg_b**0.5, ord=1) + return f + + +def eigenvalue_infidelity(a, b, gauge_invariant=True) -> _np.floating: + return 1 - eigenvalue_fidelity(a, b, gauge_invariant) + + def frobeniusdist(a, b) -> _np.floating[Any]: """ Returns the frobenius distance between arrays: ||a - b||_Fro. @@ -377,22 +394,40 @@ def jtracedist(a, b, mx_basis='pp'): # Jamiolkowski trace distance: Tr(|J(a)-J return tracedist(JA, JB) +def is_trace_preserving(a: _np.ndarray, mx_basis: BasisLike='pp', tol=1e-14) -> bool: + dim = int(_np.round( a.size**0.5 )) + udim = int(_np.round( dim**0.5 )) + basis = _Basis.cast(mx_basis, dim=dim) + if basis.first_element_is_identity: + checkone = _np.isclose(a[0,0], 1.0, rtol=tol, atol=tol) + checktwo = _np.allclose(a[1:], 0.0, atol=tol) + return checkone and checktwo + # else, check that the adjoint of a is unital. + I_mat = _np.eye(udim) + I_vec = _bt.stdmx_to_vec(I_mat, basis) + if _np.isrealobj(a): + expect_I_vec = a.T @ I_vec + else: + expect_I_vec = a.T.conj() @ I_vec + atol = tol * udim + check = float(_np.linalg.norm(I_vec - expect_I_vec)) <= atol + return check + + def entanglement_fidelity(a, b, mx_basis: BasisLike='pp', is_tp=None, is_unitary=None): """ Returns the "entanglement" process fidelity between gate matrices. - This is given by: - - `F = Tr( sqrt{ sqrt(J(a)) * J(b) * sqrt(J(a)) } )^2` - - where J(.) is the Jamiolkowski isomorphism map that maps a operation matrix - to it's corresponding Choi Matrix. + The entanglement fidelity of (a, b) is defined as the fidelity of (J(a), J(b)), + where J(.) is the Jamiolkowski isomorphism map. - When the both of the input matrices a and b are TP, and - the target matrix b is unitary then we can use a more efficient - formula: + When the both of the input matrices a and b are TP, and the target matrix b + is unitary, then we can use the formula - `F= Tr(a @ b.conjugate().T)/d^2` + `F = Tr(a @ b.conjugate().T)/d^2` + + That formula avoids forming J(a) or J(b), and it avoids the eigenvalue computations + needed to evaluate fidelity in the general case. Parameters ---------- @@ -439,26 +474,14 @@ def entanglement_fidelity(a, b, mx_basis: BasisLike='pp', is_tp=None, is_unitary #if the tp flag isn't set we'll calculate whether it is true here if is_tp is None: - def is_tp_fn(x): - return _np.isclose(x[0, 0], 1.0) and _np.allclose(x[0, 1:d2], 0) - - is_tp= (is_tp_fn(a) and is_tp_fn(b)) + is_tp = is_trace_preserving(a, mx_basis) and is_trace_preserving(b, mx_basis) #if the unitary flag isn't set we'll calculate whether it is true here if is_unitary is None: - is_unitary= _np.allclose(_np.identity(d2, 'd'), _np.dot(b, b.conjugate().T)) + is_unitary = _np.allclose(_np.identity(d2, 'd'), b @ b.T.conj()) if is_tp and is_unitary: # then assume TP-like gates & use simpler formula - #old version, slower than einsum - #TrLambda = _np.trace(_np.dot(a, b.conjugate().T)) # same as using _np.linalg.inv(b) - - #Use einsum black magic to only calculate the diagonal elements - #if the basis is either pp or gm we know the elements are real-valued, so we - #don't need to take the conjugate - if mx_basis=='pp' or mx_basis=='gm': - TrLambda = _np.einsum('ij,ji->',a, b.T) - else: - TrLambda = _np.einsum('ij,ji->',a, b.conjugate().T) + TrLambda = _np.vdot(b, a) # == Trace(b.T.conj() @ a) == Trace(a @ b.T.conj()) return TrLambda / d2 JA = _jam.jamiolkowski_iso(a, mx_basis, mx_basis) @@ -637,6 +660,7 @@ def entanglement_infidelity(a, b, mx_basis: BasisLike = 'pp', is_tp=None, is_uni """ return 1 - entanglement_fidelity(a, b, mx_basis, is_tp, is_unitary) + def generator_infidelity(a, b, mx_basis = 'pp'): """ Returns the generator infidelity between a and b, where b is the "target" operation. @@ -687,6 +711,7 @@ def generator_infidelity(a, b, mx_basis = 'pp'): return _np.real_if_close(gen_infid) + def gateset_infidelity(model, target_model, itype='EI', weights=None, mx_basis=None, is_tp=None, is_unitary=None): """ diff --git a/pygsti/tools/sdptools.py b/pygsti/tools/sdptools.py index 071195a20..11da0dcc6 100644 --- a/pygsti/tools/sdptools.py +++ b/pygsti/tools/sdptools.py @@ -14,7 +14,7 @@ import numpy as np -from typing import Union, List, Tuple, TYPE_CHECKING +from typing import Union, List, Tuple, Optional, TYPE_CHECKING if TYPE_CHECKING: import cvxpy as cp ExpressionLike = Union[cp.Expression, np.ndarray] @@ -159,23 +159,30 @@ def cptp_superop_variable(purestate_dim: int, basis: BasisLike) -> Tuple[cp.Expr return X, constraints -def diamond_distance_projection_model(superop: np.ndarray, basis: Basis, leakfree: bool=False, seepfree: bool=False, n_leak: int=0, cptp: bool=True, subspace_diamond: bool=False): +def diamond_distance_projection_model(superop: np.ndarray, basis: BasisLike, leakfree: bool=False, seepfree: bool=False, n_leak: Optional[int]=ModuleNotFoundError, cptp: bool=True, subspace_diamond: bool=False): assert CVXPY_ENABLED dim_mixed = superop.shape[0] dim_pure = int(np.sqrt(dim_mixed)) assert dim_pure**2 == dim_mixed constraints = [] + from pygsti.tools.leakage import leading_dxd_submatrix_basis_vectors, projective_measurement_subspace_basis_vectors + + if not isinstance(basis, Basis): + basis = Basis.cast(basis, dim=dim_mixed) + if n_leak is None: + U = projective_measurement_subspace_basis_vectors(basis.ellookup['I'], basis) + else: + dim_pure_compsub = dim_pure - n_leak + U = leading_dxd_submatrix_basis_vectors(dim_pure_compsub, dim_pure, basis) # type: ignore + if cptp: proj_superop, cons = cptp_superop_variable(dim_pure, basis) constraints.extend(cons) else: proj_superop = cp.Variable((dim_mixed, dim_mixed)) + diamondnorm_arg = superop - proj_superop if (leakfree or seepfree or subspace_diamond): - assert n_leak == 1 - from pygsti.tools.leakage import leading_dxd_submatrix_basis_vectors - dim_pure_compsub = dim_pure - n_leak - U = leading_dxd_submatrix_basis_vectors(dim_pure_compsub, dim_pure, basis) P = U @ U.T.conj() I = np.eye(dim_mixed) if leakfree: diff --git a/test/unit/tools/test_matrixtools.py b/test/unit/tools/test_matrixtools.py index 3b2b20a75..0eb31c209 100644 --- a/test/unit/tools/test_matrixtools.py +++ b/test/unit/tools/test_matrixtools.py @@ -24,6 +24,25 @@ def test_is_pos_def(self): self.assertTrue(mt.is_pos_def(pos_mx)) self.assertFalse(mt.is_pos_def(non_pos_mx)) + def test_is_projector(self): + self.assertTrue(mt.is_projector(np.zeros((3, 3)))) + I2 = np.eye(2) + self.assertTrue(mt.is_projector(I2)) + I2[1,1] = 0.9999998 + self.assertFalse(mt.is_projector(I2)) + self.assertTrue( mt.is_projector(I2, tol=np.inf)) + A = np.array([ + [-0.653829-0.313923j, -0.129614+1.458021j, 0.783975+1.960258j], + [ 1.493431+1.801635j, -1.259066+1.315104j, 1.513924+0.35738j ], + [ 1.345875-1.208319j, 0.781311-0.004454j, 0.264456+0.656475j] + ]) + Q = np.linalg.qr(A)[0] + for i in range(1, 4): + U = Q[:,:i].reshape((3, i)) + P = U @ U.T.conj() + self.assertTrue(mt.is_projector(P)) + return + def test_mx_to_string(self): mx = np.array([[ 1, 1+2j], [1-2j, 3]], 'complex') @@ -192,3 +211,4 @@ def test_prime_factors(self): self.assertEqual(mt.prime_factors(7), [7]) self.assertEqual(mt.prime_factors(10), [2, 5]) self.assertEqual(mt.prime_factors(12), [2, 2, 3]) +