diff --git a/jupyter_notebooks/Tutorials/objects/advanced/Instruments.ipynb b/jupyter_notebooks/Tutorials/objects/advanced/Instruments.ipynb index 3e4093b39..989bf4acb 100644 --- a/jupyter_notebooks/Tutorials/objects/advanced/Instruments.ipynb +++ b/jupyter_notebooks/Tutorials/objects/advanced/Instruments.ipynb @@ -12,13 +12,15 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "import pygsti\n", "from pygsti.modelpacks import smq1Q_XYI as std\n", - "import numpy as np" + "import numpy as np\n", + "from pygsti.modelmembers.instruments import Instrument\n", + "from pprint import pprint" ] }, { @@ -31,58 +33,40 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "#Make a copy so we don't modify the original\n", - "target_model = std.target_model()\n", + "mdl_ideal = std.target_model()\n", "\n", - "#Create and add the ideal instrument\n", - "E0 = target_model.effects['0']\n", - "E1 = target_model.effects['1']\n", - " # Alternate indexing that uses POVM label explicitly\n", - " # E0 = target_model['Mdefault']['0'] # 'Mdefault' = POVM label, '0' = effect label\n", - " # E1 = target_model['Mdefault']['1']\n", - "Gmz_plus = np.dot(E0,E0.T) #note effect vectors are stored as column vectors\n", + "# Create and add the ideal instrument\n", + "E0 = mdl_ideal.effects['0']\n", + "E1 = mdl_ideal.effects['1']\n", + "Gmz_plus = np.dot(E0,E0.T) # note effect vectors are stored as column vectors\n", "Gmz_minus = np.dot(E1,E1.T)\n", - "target_model[('Iz',0)] = pygsti.modelmembers.instruments.Instrument({'p0': Gmz_plus, 'p1': Gmz_minus})\n", - "\n", - "#For later use, record the identity POVM vector\n", - "povm_ident = target_model.effects['0'] + target_model.effects['1'] " + "inst = Instrument({'p0': Gmz_plus, 'p1': Gmz_minus})\n", + "mdl_ideal[('Iz',0)] = inst" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In order to generate some simulated data later on, we'll now create a noisy version of `target_model` by depolarizing the state preparation, gates, and POVM, and also rotating the basis that is measured by the instrument and POVM." + "In order to generate some simulated data later on, we'll now create a noisy version of `mdl_ideal` by depolarizing the state preparation, gates, and POVM, and also rotating the basis that is measured by the instrument and POVM." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ - "mdl_noisy = target_model.depolarize(op_noise=0.01, spam_noise=0.01)\n", - "mdl_noisy.effects.depolarize(0.01) #because above call only depolarizes the state prep, not the POVMs\n", - "\n", - "# add a rotation error to the POVM\n", - "Uerr = pygsti.rotation_gate_mx([0,0.02,0])\n", - "E0 = np.dot(mdl_noisy['Mdefault']['0'].T,Uerr).T\n", - "E1 = povm_ident - E0\n", - "mdl_noisy.povms['Mdefault'] = pygsti.modelmembers.povms.UnconstrainedPOVM({'0': E0, '1': E1},\n", - " evotype='default')\n", - "\n", - "# Use the same rotated effect vectors to \"rotate\" the instrument Iz too\n", - "E0 = mdl_noisy.effects['0']\n", - "E1 = mdl_noisy.effects['1']\n", - "Gmz_plus = np.dot(E0,E0.T)\n", - "Gmz_minus = np.dot(E1,E1.T)\n", - "mdl_noisy[('Iz',0)] = pygsti.modelmembers.instruments.Instrument({'p0': Gmz_plus, 'p1': Gmz_minus})\n", - "\n", - "#print(mdl_noisy) #print the model" + "mdl_noisy = mdl_ideal.depolarize(op_noise=0.005, spam_noise=0.01)\n", + "mdl_noisy = mdl_noisy.rotate(max_rotate=0.025, seed=2048)\n", + "mdl_noisy.effects.depolarize(0.01)\n", + "mdl_ideal.convert_members_inplace('CPTPLND')\n", + "mdl_noisy.convert_members_inplace('CPTPLND')" ] }, { @@ -95,20 +79,66 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{('p0', '0'): 0.5000000000000003,\n", + " ('p0', '1'): 1.4597476004180728e-17,\n", + " ('p1', '0'): 5.073268178784763e-18,\n", + " ('p1', '1'): 0.5}\n", + "\n", + "{('p0', '0'): 0.481367923790963,\n", + " ('p0', '1'): 0.0024189342904068554,\n", + " ('p1', '0'): 0.0025810657095932145,\n", + " ('p1', '1'): 0.5136320762090365}\n", + "\n" + ] + } + ], "source": [ - "dict(target_model.probabilities( pygsti.circuits.Circuit(( ('Gxpi2',0) , ('Iz',0) )) ))" + "for mdl in [mdl_ideal, mdl_noisy]:\n", + " pprint(dict(mdl.probabilities( pygsti.circuits.Circuit(( ('Gxpi2',0) , ('Iz',0) )) )))\n", + " print()" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 16, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{('p0', 'p0', '0'): 0.5000000000000002,\n", + " ('p0', 'p0', '1'): -1.6613637599827203e-18,\n", + " ('p0', 'p1', '0'): -1.1185571585378685e-17,\n", + " ('p0', 'p1', '1'): 0.4999999999999999,\n", + " ('p1', 'p0', '0'): 1.625883976416346e-17,\n", + " ('p1', 'p0', '1'): 1.8939874217871704e-34,\n", + " ('p1', 'p1', '0'): 4.5374861265545935e-34,\n", + " ('p1', 'p1', '1'): 1.625883976416347e-17}\n", + "\n", + "{('p0', 'p0', '0'): 0.4775917799207901,\n", + " ('p0', 'p0', '1'): 0.002399958693069308,\n", + " ('p0', 'p1', '0'): 0.0025624981205110728,\n", + " ('p0', 'p1', '1'): 0.5099371259816915,\n", + " ('p1', 'p0', '0'): 0.0038579004920827257,\n", + " ('p1', 'p0', '1'): 1.938643463358172e-05,\n", + " ('p1', 'p1', '0'): 1.8156751786106926e-05,\n", + " ('p1', 'p1', '1'): 0.003613193605435256}\n", + "\n" + ] + } + ], "source": [ - "dict(target_model.probabilities( pygsti.circuits.Circuit(( ('Iz',0), ('Gxpi2',0) , ('Iz',0) )) ))" + "for mdl in [mdl_ideal, mdl_noisy]:\n", + " pprint(dict(mdl.probabilities( pygsti.circuits.Circuit(( ('Iz',0), ('Gxpi2',0) , ('Iz',0) )) )))\n", + " print()" ] }, { @@ -120,11 +150,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 17, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "probs = {('0',): 0.5000000000000004, ('1',): 0.5}\n", + "probs['0'] = 0.5000000000000004\n", + "probs[('0',)] = 0.5000000000000004\n" + ] + } + ], "source": [ - "probs = target_model.probabilities( pygsti.circuits.Circuit([('Gxpi2',0)]) )\n", + "probs = mdl_ideal.probabilities( pygsti.circuits.Circuit([('Gxpi2',0)]) )\n", "print(\"probs = \",dict(probs))\n", "print(\"probs['0'] = \", probs['0']) #This works...\n", "print(\"probs[('0',)] = \", probs[('0',)]) # and so does this." @@ -135,158 +175,175 @@ "metadata": {}, "source": [ "## Performing tomography\n", - "\n", - "### Simulated data generation\n", - "Now let's perform tomography on a model that includes instruments. First, we'll generate some data using `mdl_noisy` in exactly the same way as we would for any other model:" + "Now let's perform tomography on a model that includes instruments. \n", + "First, we'll build an experiment design. This notebook's experiment design makes a minor modification to the default design for our working modelpack (smq1Q_XYI). " ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "germs = std.germs()\n", - "germs += [pygsti.circuits.Circuit([('Iz', 0)])] # add the instrument as a germ.\n", - "\n", - "prep_fiducials = std.prep_fiducials()\n", - "meas_fiducials = std.meas_fiducials()\n", - "max_lengths = [1] # keep it simple & fast\n", - "\n", - "lsgst_list = pygsti.circuits.create_lsgst_circuits(\n", - " mdl_noisy,prep_fiducials,meas_fiducials,germs,max_lengths)\n", - "\n", - "#print(\"LinearOperator sequences:\")\n", - "#print(lsgst_list) #note that this contains LGST strings with \"Iz\"\n", - "\n", - "#Create the DataSet\n", - "ds = pygsti.data.simulate_data(mdl_noisy,lsgst_list,1000,'multinomial',seed=2018)\n", - "\n", - "#Write it to a text file to demonstrate the format:\n", - "pygsti.io.write_dataset(\"../../tutorial_files/intermediate_meas_dataset.txt\",ds)" + "germs += [pygsti.circuits.Circuit([('Iz', 0)])]\n", + "ed = std.create_gst_experiment_design(max_max_length=8, germs=germs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Notice the format of [intermediate_meas_dataset.txt](../../tutorial_files/intermediate_meas_dataset.txt), which includes a column for each distinct outcome tuple. Since not all experiments contain data for all outcome tuples, the `\"--\"` is used as a placeholder. Now that the data is generated, we run LGST or LSGST just like we would for any other model:\n", - "\n", - "### LGST" + "### Simulated data generation\n", + "Next, we generate data using `mdl_noisy` in exactly the same way as we would for any other model. We write the data to disk so you can see how datasets look when they include measurement data." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 19, "metadata": {}, "outputs": [], "source": [ - "#Run LGST\n", - "mdl_lgst = pygsti.run_lgst(ds, prep_fiducials, meas_fiducials, target_model)\n", - "#print(mdl_lgst)\n", - "\n", - "#Gauge optimize the result to the true data-generating model (mdl_noisy),\n", - "# and compare. Mismatch is due to finite sample noise.\n", - "mdl_lgst_opt = pygsti.gaugeopt_to_target(mdl_lgst,mdl_noisy)\n", - "print(mdl_noisy.strdiff(mdl_lgst_opt))\n", - "print(\"Frobdiff after GOpt = \",mdl_noisy.frobeniusdist(mdl_lgst_opt))" + "ds = pygsti.data.simulate_data(mdl_noisy, ed.all_circuits_needing_data, 10_000, 'multinomial', seed=2018)\n", + "pygsti.io.write_dataset(\"../../tutorial_files/intermediate_meas_dataset.txt\", ds)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Long-sequence GST\n", - "Instruments just add parameters to a `Model` like gates, state preparations, and POVMs do. The total number of parameters in our model is \n", + "Notice the format of [intermediate_meas_dataset.txt](../../tutorial_files/intermediate_meas_dataset.txt), which includes a column for each distinct outcome tuple. Since not all experiments contain data for all outcome tuples, the `\"--\"` is used as a placeholder. Now that the data is generated, we run LGST or LSGST just like we would for any other model:\n", "\n", - "$4$ (prep) + $2\\times 4$ (2 effects) + $5\\times 16$ (3 gates and 2 instrument members) $ = 92$." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "target_model.num_params" + "### Running the GST fit and generating a report" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 20, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-- Std Practice: Iter 1 of 2 (full TP) --: \n", + " --- Iterative GST: [##################################################] 100.0% 592 circuits ---\n", + "-- Std Practice: Iter 2 of 2 (CPTPLND) --: \n", + " --- Iterative GST: [##################################################] 100.0% 592 circuits ---\n" + ] + } + ], "source": [ - "#Run long sequence GST\n", - "results = pygsti.run_long_sequence_gst(ds,target_model,prep_fiducials,meas_fiducials,germs,max_lengths)" + "from pygsti.protocols import StandardGST, ProtocolData\n", + "gst = StandardGST(\n", + " modes=('full TP', 'CPTPLND'), target_model=mdl_ideal, verbosity=2\n", + ")\n", + "pd = ProtocolData(ed, ds)\n", + "res = gst.run(pd)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 21, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running idle tomography\n", + "Computing switchable properties\n", + "Found standard clifford compilation from smq1Q_XYI\n", + "Found standard clifford compilation from smq1Q_XYI\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/models/model.py:145: UserWarning:\n", + "\n", + "Model.num_modeltest_params could not obtain number of *non-gauge* parameters - using total instead\n", + "\n", + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/tools/optools.py:152: UserWarning:\n", + "\n", + "The input matrix a is not trace-1 up to tolerance 1.4901161193847656e-08. Beware result!\n", + "\n", + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/tools/optools.py:154: UserWarning:\n", + "\n", + "The input matrix b is not trace-1 up to tolerance 1.4901161193847656e-08. Beware result!\n", + "\n", + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/tools/optools.py:179: UserWarning:\n", + "\n", + "\n", + " Input matrix is not PSD up to tolerance 1.4901161193847656e-08.\n", + " We'll project out the bad eigenspaces to only work with the PSD part.\n", + " \n", + "\n", + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/report/workspacetables.py:2741: RuntimeWarning:\n", + "\n", + "divide by zero encountered in log\n", + "\n", + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/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" + ] + } + ], "source": [ - "#Compare estimated model (after gauge opt) to data-generating one\n", - "mdl_est = results.estimates['GateSetTomography'].models['go0']\n", - "mdl_est_opt = pygsti.gaugeopt_to_target(mdl_est,mdl_noisy)\n", - "print(\"Frobdiff after GOpt = \", mdl_noisy.frobeniusdist(mdl_est_opt))" + "from pygsti.report import construct_standard_report\n", + "report_dir = '../../tutorial_files/cptp-instrument-report'\n", + "report_object = construct_standard_report(res, title='CPTP intrument GST')\n", + "report_object.write_html(report_dir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The same analysis can be done for a trace-preserving model, whose instruments are constrained to *add* to a perfectly trace-preserving map. The number of parameters in the model are now: \n", + "**Thats it!** You've done tomography on a model with intermediate measurments (instruments).\n", "\n", - "$3$ (prep) + $1\\times 4$ (effect and complement) + $3\\times 12$ (3 gates) + $(2\\times 16 - 3)$ (TP instrument) $ = 71$" + "As a bonus, the code below checks for violation of complete positivity of the instrument operations from the two fits. We see that the CPTP fit has no violation, while the full TP fit has small violations." ] }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mdl_targetTP = target_model.copy()\n", - "mdl_targetTP.set_all_parameterizations(\"full TP\")\n", - "print(\"POVM type = \",type(mdl_targetTP[\"Mdefault\"]),\" Np=\",mdl_targetTP[\"Mdefault\"].num_params)\n", - "print(\"Instrument type = \",type(mdl_targetTP[(\"Iz\",0)]),\" Np=\",mdl_targetTP[(\"Iz\",0)].num_params)\n", - "print(\"Number of model parameters = \", mdl_targetTP.num_params)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "resultsTP = pygsti.run_long_sequence_gst(ds,mdl_targetTP,prep_fiducials,meas_fiducials,germs,max_lengths)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#Again compare estimated model (after gauge opt) to data-generating one\n", - "mdl_est = resultsTP.estimates['GateSetTomography'].models['go0']\n", - "mdl_est_opt = pygsti.gaugeopt_to_target(mdl_est,mdl_noisy)\n", - "print(\"Frobdiff after GOpt = \", mdl_noisy.frobeniusdist(mdl_est_opt))" - ] - }, - { - "cell_type": "markdown", + "execution_count": 22, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Iz:0\n", + "p0: 0.0009697360237423556\n", + "p1: 2.6230206864367674e-05\n", + "\n", + "Iz:0\n", + "p0: 0\n", + "p1: 0\n", + "\n" + ] + } + ], "source": [ - "**Thats it!** You've done tomography on a model with intermediate measurments (instruments)." + "for mdl in [res.estimates['full TP'].models['stdgaugeopt'], res.estimates['CPTPLND'].models['stdgaugeopt']]:\n", + " for instlbl, inst in mdl.instruments.items():\n", + " print(instlbl)\n", + " for ioplbl, iop in inst.items():\n", + " violation = max(0, pygsti.tools.sum_of_negative_choi_eigenvalues_gate(iop.to_dense(), 'pp'))\n", + " print(ioplbl + ': ' + str(violation) )\n", + " print()" ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "pgdev311", "language": "python", "name": "python3" }, @@ -300,7 +357,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.13" + "version": "3.11.13" } }, "nbformat": 4, diff --git a/pygsti/baseobjs/__init__.py b/pygsti/baseobjs/__init__.py index 44eedfded..4e78d505f 100644 --- a/pygsti/baseobjs/__init__.py +++ b/pygsti/baseobjs/__init__.py @@ -13,7 +13,7 @@ from .smartcache import SmartCache from .verbosityprinter import VerbosityPrinter from .profiler import Profiler -from .basis import Basis, BuiltinBasis, ExplicitBasis, TensorProdBasis, DirectSumBasis +from .basis import Basis, BuiltinBasis, ExplicitBasis, TensorProdBasis, DirectSumBasis, BasisLike from .label import Label, CircuitLabel from .nicelyserializable import NicelySerializable from .mongoserializable import MongoSerializable diff --git a/pygsti/modelmembers/instruments/__init__.py b/pygsti/modelmembers/instruments/__init__.py index 8ef2c0ed2..e3331013b 100644 --- a/pygsti/modelmembers/instruments/__init__.py +++ b/pygsti/modelmembers/instruments/__init__.py @@ -13,8 +13,18 @@ from .instrument import Instrument from .tpinstrument import TPInstrument from .tpinstrumentop import TPInstrumentOp +from ..operations.cptrop import RootConjOperator, SummedOperator -from pygsti.tools import optools as _ot +import warnings as _warnings +import scipy.linalg as _la +import numpy as _np +from pygsti.tools import optools as _ot, basistools as _bt +from types import NoneType +from pygsti.baseobjs.label import Label +from pygsti.baseobjs.basis import Basis +from pygsti.modelmembers import operations as _op +from pygsti.modelmembers import povms as _pv +from pygsti.modelmembers.povms.basepovm import _BasePOVM # Avoid circular import import pygsti.modelmembers as _mm @@ -37,7 +47,7 @@ def instrument_type_from_op_type(op_type): # Limited set (only matching what is in convert) instr_conversion = { - 'auto': 'full', + 'auto': 'full TP', 'static unitary': 'static unitary', 'static clifford': 'static clifford', 'static': 'static', @@ -45,19 +55,29 @@ def instrument_type_from_op_type(op_type): 'full TP': 'full TP', 'full CPTP': 'full CPTP', 'full unitary': 'full unitary', + 'GLND': 'full TP', + # ^ It's pretty harmless to associate GLND operations with "full TP" + # instruments. In both cases we're relaxing positivity constraints. + 'CPTPLND': 'CPTPLND' } instr_type_preferences = [] for typ in op_type_preferences: - instr_type = None - if _ot.is_valid_lindblad_paramtype(typ): - # Lindblad types are passed through as TP only (matching current convert logic) - instr_type = "full TP" - else: - instr_type = instr_conversion.get(typ, None) + instr_type = instr_conversion.get(typ, None) - if instr_type is None: - continue + if instr_type is None and _ot.is_valid_lindblad_paramtype(typ): + # NOTE: need to update the message below if more lindblad + # types are added as keys to the instr_conversion dict. + msg = \ + f""" + Operation type {typ} is a Lindblad parameterization, but + is neither 'GLND' or 'CPTPLND'. That means you might be + asking for a reduced-order model that's a subset of all + CPTP models. We don't support that parameterization at this + time. We're falling back to a 'full TP' parameterization! + """ + _warnings.warn(msg) + instr_type = 'full TP' # non-CPTPLND falls back to full TP. if instr_type not in instr_type_preferences: instr_type_preferences.append(instr_type) @@ -110,26 +130,153 @@ def convert(instrument, to_type, basis, ideal_instrument=None, flatten_structure The converted instrument, usually a distinct object from the object passed as input. """ - to_types = to_type if isinstance(to_type, (tuple, list)) else (to_type,) # HACK to support multiple to_type values + if not isinstance(to_type, str): + if len(to_type) > 1: + raise ValueError(f"Expected to_type to be a string, but got {to_type}") + to_type = to_type[0] + assert isinstance(to_type, str) + destination_types = {'full TP': TPInstrument} - NoneType = type(None) - - for to_type in to_types: - try: - if isinstance(instrument, destination_types.get(to_type, NoneType)): - return instrument - - if to_type == "full TP": - return TPInstrument(list(instrument.items()), instrument.evotype, instrument.state_space) - elif to_type in ("full", "static", "static unitary"): - from ..operations import convert as _op_convert - ideal_items = dict(ideal_instrument.items()) if (ideal_instrument is not None) else {} - members = [(k, _op_convert(g, to_type, basis, ideal_items.get(k, None), flatten_structure)) - for k, g in instrument.items()] - return Instrument(members, instrument.evotype, instrument.state_space) + + if isinstance(instrument, destination_types.get( to_type, NoneType ) ): + return instrument + + if to_type == "full TP": + inst_arrays = dict() + for k, v in instrument.items(): + if hasattr(v, 'to_dense'): + inst_arrays[k] = v.to_dense('HilbertSchmidt') else: - raise ValueError("Cannot convert an instrument to type %s" % to_type) - except: - pass # try next to_type + inst_arrays[k] = v + return TPInstrument(list(inst_arrays.items()), instrument.evotype, instrument.state_space) + + if to_type in ("full", "static", "static unitary"): + from ..operations import convert as _op_convert + ideal_items = dict(ideal_instrument.items()) if (ideal_instrument is not None) else {} + members = [(k, _op_convert(g, to_type, basis, ideal_items.get(k, None), flatten_structure)) + for k, g in instrument.items()] + return Instrument(members, instrument.evotype, instrument.state_space) + + if to_type == 'CPTPLND': + op_arrays = {k: v.to_dense('HilbertSchmidt') for (k,v) in instrument.items()} + inst = cptp_instrument(op_arrays, basis) + return inst + + raise ValueError("Cannot convert an instrument to type %s" % to_type) + + +def cptp_instrument(op_arrays: dict[str, _np.ndarray], basis: Basis, + error_tol: float = 1e-6, trunc_tol: float = 1e-7, + povm_errormap=None) -> Instrument: + """ + Construct a CPTPLND-parameterized instrument from a dictionary of dense superoperator arrays. + + Each operator in `op_arrays` must be a completely positive trace-reducing (CPTR) map. + The function Kraus-decomposes each operator, polar-decomposes each Kraus operator into + a unitary and a positive semidefinite part, then assembles the result using + CPTPLND-parameterized unitaries and a shared `ComposedPOVM` that enforces the + completeness (trace-preservation) constraint across all instrument outcomes. + + The construction follows these steps: + + 1. For each outcome label, compute a minimal Kraus decomposition and polar-decompose + each Kraus operator ``K = u p`` (``p`` PSD, ``u`` unitary). + 2. Represent each Kraus term as the composition of a ``RootConjOperator`` (encoding + ``ρ ↦ p² ρ p²``, parameterized by the POVM effect for ``p²``) and a + CPTPLND-parameterized unitary channel. + 3. Collect all POVM effects into a shared ``ComposedPOVM`` whose error map is + constrained to be CPTPLND. This enforces that the effects sum to ≤ I, which + ensures the full instrument is CPTR. + 4. If an outcome requires more than one Kraus term, wrap the summands in a + ``SummedOperator``; otherwise use the single term directly. + + Parameters + ---------- + op_arrays : dict of str → numpy.ndarray + Mapping from outcome label to dense superoperator matrix (in the given basis). + Each matrix must represent a CPTR map; the function raises or warns if the + Kraus decomposition is inconsistent with this. + basis : Basis or str + The operator basis in which the superoperators are expressed (e.g. ``'pp'``). + error_tol : float, optional + Tolerance for eigenvalue errors in the minimal Kraus decomposition. Eigenvalues + below ``-error_tol`` cause an error to be raised. Default ``1e-6``. + trunc_tol : float, optional + Truncation tolerance for the minimal Kraus decomposition. Eigenvalues between + ``-error_tol`` and ``trunc_tol`` are silently set to zero. Default ``1e-7``. + povm_errormap : LinearOperator, optional + A CPTPLND-parameterized error map to use as the shared POVM error map. + If ``None`` (default), the identity channel is used as the starting point + and converted to CPTPLND parameterization. + + Returns + ------- + Instrument + An ``Instrument`` whose member operations are CPTR maps sharing a common set + of CPTPLND parameters. The instrument is not itself a ``TPInstrument`` because + the completeness constraint is enforced implicitly through the shared POVM. + """ + unitaries = dict() + effects = dict() + + # Step 1. Build CPTPLND-parameterized unitaries and static POVM effects + # from Kraus decompositions of the provided operators. + for oplbl, op in op_arrays.items(): + krausops = _ot.minimal_kraus_decomposition(op, basis, error_tol, trunc_tol) + cur_effects = [] + cur_unitaries = [] + for K in krausops: + # Define F(rho) := K rho K^\\dagger. If we compute a polar decomposition + # K = u p, then we can express F(rho) = u (p rho p) u^\\dagger, which is + # a composition of a unitary channel and the "RootConjOperator" of p@p. + u, p = _la.polar(K, side='right') + u_linop = _op.StaticUnitaryOp(u, basis) # type: ignore + u_cptp = _op.convert(u_linop, 'CPTPLND', basis) + cur_unitaries.append(u_cptp) + E_superket = _bt.stdmx_to_vec(p @ p, basis) + E = _pv.StaticPOVMEffect(E_superket) + cur_effects.append(E) + effects[oplbl] = cur_effects + unitaries[oplbl] = cur_unitaries + + # Step 2. Build a CPTPLND-parameterized POVM from the static POVM effects. + # + # We can't use povms.convert(...) because we might have more + # effects than the Hilbert space dimension. + # + base_povm_effects = {} + for oplbl, cur_effects in effects.items(): + for i, E in enumerate(cur_effects): + elbl = Label((oplbl, i)) + base_povm_effects[elbl] = E + base_povm = _BasePOVM(base_povm_effects, preserve_sum=False) + udim = int(_np.round(basis.dim ** 0.5)) + if povm_errormap is None: + I_ideal = _op.StaticUnitaryOp(_np.eye(udim), basis) # type: ignore + I_cptplnd = _op.convert(I_ideal, 'CPTPLND', basis) + povm_errormap = I_cptplnd.factorops[1] + povm_cptp = _pv.ComposedPOVM(povm_errormap, base_povm) + + # Step 3. Assemble the CPTPLND-parameterized unitaries and POVM + # effects to represent each operator as a completely-positive + # trace-reducing map. + # + inst_ops = dict() + for lbl in op_arrays: + cur_effects = effects[lbl] + cur_unitaries = unitaries[lbl] + + op_summands = [] + for i, U_cptplnd in enumerate(cur_unitaries): + E_cptp = povm_cptp[Label((lbl, i))] + rcop = RootConjOperator(E_cptp, basis) + cptr = _op.ComposedOp([rcop, U_cptplnd]) + op_summands.append(cptr) + if len(op_summands) == 1: + op = op_summands[0] + else: + op = SummedOperator(op_summands, basis) - raise ValueError("Could not convert instrument to to type(s): %s" % str(to_types)) + inst_ops[lbl] = op + inst = Instrument(inst_ops) + return inst diff --git a/pygsti/modelmembers/instruments/instrument.py b/pygsti/modelmembers/instruments/instrument.py index 2e0a5a3a0..38e01ba5c 100644 --- a/pygsti/modelmembers/instruments/instrument.py +++ b/pygsti/modelmembers/instruments/instrument.py @@ -19,6 +19,9 @@ from pygsti.baseobjs import statespace as _statespace from pygsti.tools import matrixtools as _mt from pygsti.tools import slicetools as _slct +from pygsti.tools import basistools as _bt +from pygsti.tools import optools as _ot +from pygsti.baseobjs.basis import Basis as _Basis from pygsti.baseobjs.label import Label as _Label from pygsti.baseobjs.statespace import StateSpace as _StateSpace @@ -378,13 +381,14 @@ def acton(self, state): staterep = state._rep outcome_probs_and_states = _collections.OrderedDict() + for lbl, element in self.items(): output_rep = element._rep.acton(staterep) - output_unnormalized_state = output_rep.to_dense() - prob = output_unnormalized_state[0] * state.dim**0.25 - output_normalized_state = output_unnormalized_state / prob # so [0]th == 1/state_dim**0.25 - outcome_probs_and_states[lbl] = (prob, _state.StaticState(output_normalized_state, self.evotype, - self.state_space)) + unnormalized_state_array = output_rep.to_dense() + prob = _ot.superket_trace(unnormalized_state_array, element.basis) + output_state_array = unnormalized_state_array / prob + output_state = _state.StaticState(output_state_array, self.evotype, self.state_space) + outcome_probs_and_states[lbl] = (prob, output_state) return outcome_probs_and_states diff --git a/pygsti/modelmembers/operations/__init__.py b/pygsti/modelmembers/operations/__init__.py index d0f340c4f..e585e36ea 100644 --- a/pygsti/modelmembers/operations/__init__.py +++ b/pygsti/modelmembers/operations/__init__.py @@ -39,11 +39,13 @@ from .stochasticop import StochasticNoiseOp from .lindbladcoefficients import LindbladCoefficientBlock as _LindbladCoefficientBlock from .affineshiftop import AffineShiftOp +from .cptrop import RootConjOperator, SummedOperator from pygsti.baseobjs import statespace as _statespace from pygsti.tools import basistools as _bt from pygsti.tools import optools as _ot from pygsti import SpaceT + def create_from_unitary_mx(unitary_mx, op_type, basis='pp', stdname=None, evotype='default', state_space=None): """ TODO: docstring - note that op_type can be a list/tuple of types in order of precedence """ op_type_preferences = verbose_type_from_op_type(op_type) diff --git a/pygsti/modelmembers/operations/cptrop.py b/pygsti/modelmembers/operations/cptrop.py new file mode 100644 index 000000000..886ad30e6 --- /dev/null +++ b/pygsti/modelmembers/operations/cptrop.py @@ -0,0 +1,202 @@ +""" +Operators used in CPTR (completely positive trace-reducing) parameterizations. + +A completely positive trace-reducing (CPTR) map is a completely positive map +whose Kraus operators satisfy sum_i K_i^† K_i ≤ I (i.e., trace is non-increasing). +Every CPTR map can be decomposed as a CPTP map followed by a "root-conjugation" +by a positive semidefinite matrix E with 0 ≤ E ≤ I. + +The two classes here, `RootConjOperator` and `SummedOperator`, are `LinearOperator` +subclasses used as building blocks when constructing CPTPLND-parameterized instruments +via :func:`pygsti.modelmembers.instruments.cptp_instrument`. +""" +#*************************************************************************************************** +# Copyright 2015, 2019, 2025 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights +# in this software. +# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except +# in compliance with the License. You may obtain a copy of the License at +# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. +#*************************************************************************************************** + +import numpy as _np +from pygsti.pgtypes import SpaceT +from pygsti.baseobjs.basis import Basis +from pygsti.modelmembers.povms.effect import POVMEffect +from pygsti.modelmembers.operations.linearop import LinearOperator +from pygsti.tools import optools as _ot + +from typing import Union +BasisLike = Union[Basis, str] + + +class RootConjOperator(LinearOperator): + """ + A linear operator parameterized by a PSD matrix E with 0 ≤ E ≤ I, acting as ρ ↦ E^½ ρ E^½. + + Every CPTR (completely positive trace-reducing) map can be expressed as the + composition of a CPTP map with a ``RootConjOperator``. This class implements + the "root-conjugation" factor: given a POVM effect E (representing the PSD matrix + in superket form), the superoperator for ρ ↦ E^½ ρ E^½ is constructed and kept + up to date as the effect's parameters change. + + The effect E encodes the constraint 0 ≤ E ≤ I; its parameters are inherited + directly as the parameters of this operator. + + Parameters + ---------- + effect : POVMEffect + A POVM effect whose superket encodes the PSD matrix E. The effect's + parameter vector is shared (via gpindices) with this operator. + basis : Basis or str + The operator basis in which the superoperator is represented. + + Attributes + ---------- + EIGTOL_WARNING : float + Eigenvalue tolerance below which a warning is issued during construction. + EIGTOL_ERROR : float + Eigenvalue tolerance below which an error is raised during construction. + """ + + EIGTOL_WARNING = 1e-10 + EIGTOL_ERROR = 1e-8 + + def __init__(self, effect: POVMEffect, basis: BasisLike): + self._basis = Basis.cast(basis, effect.dim) + self._effect = effect + self._state_space = effect.state_space + self._evotype = effect.evotype + + dim = self._state_space.dim + self._rep = self._evotype.create_dense_superop_rep( + _np.zeros((dim, dim)), self._basis, self._state_space + ) + self._update_rep_base() + LinearOperator.__init__(self, self._rep, self._evotype) + self.init_gpindices() + + def submembers(self): + return [self._effect] + + def _update_rep_base(self): + # Analogous to TPInstrumentOp._construct_matrix. + self._rep.base.flags.writeable = True + assert(self._rep.base.shape == (self.dim, self.dim)) + effect_superket = self._effect.to_dense() + mx = _ot.rootconj_superop(effect_superket, self._basis) + self._rep.base[:] = mx + self._rep.base.flags.writeable = False + self._rep.base_has_changed() + + def deriv_wrt_params(self, wrt_filter=None): + return LinearOperator.deriv_wrt_params(self, wrt_filter) + + def has_nonzero_hessian(self): + # Not affine in its parameters. + return True + + def from_vector(self, v, close=False, dirty_value=True): + for sm, local_inds in zip(self.submembers(), self._submember_rpindices): + sm.from_vector(v[local_inds], close, dirty_value) + self._update_rep_base() + + @property + def num_params(self): + return len(self.gpindices_as_array()) + + def to_vector(self): + v = _np.empty(self.num_params, 'd') + for param_op, local_inds in zip(self.submembers(), self._submember_rpindices): + v[local_inds] = param_op.to_vector() + return v + + def to_dense(self, on_space: SpaceT = 'HilbertSchmidt') -> _np.ndarray: + assert on_space in ('HilbertSchmidt', 'minimal') + out = self._rep.base.copy() + out.flags.writeable = True + return out + + +class SummedOperator(LinearOperator): + """ + A linear operator whose superoperator is the sum of a list of constituent operators. + + ``SummedOperator`` wraps a sequence of ``LinearOperator`` objects and exposes + their pointwise sum as a single ``LinearOperator``. It is used inside + :func:`~pygsti.modelmembers.instruments.cptp_instrument` when a CPTR map requires + more than one Kraus term: each term contributes a separate summand, and the full + map is their sum. + + The parameter vector is the concatenation (with shared-index deduplication via + gpindices) of the submembers' parameter vectors. The constituent operators' ``_rep`` + objects are linked directly, so calling ``from_vector`` on a submember automatically + updates this operator's representation. + + Parameters + ---------- + operators : list of LinearOperator + The operators to sum. All must share the same evotype, state space, and + superoperator dimension. + basis : Basis or str + The operator basis (used for bookkeeping; the operators themselves must + already be expressed in this basis). + + Notes + ----- + ``deriv_wrt_params`` is not implemented for this class. + """ + + def __init__(self, operators, basis: BasisLike): + op = operators[0] + self._basis = Basis.cast(basis, op.dim) + self._operators = operators + self._state_space = op.state_space + self._evotype = op.evotype + self._subreps = [op._rep for op in self._operators] + self._rep = self._evotype.create_sum_rep(self._subreps, self._state_space) + LinearOperator.__init__(self, self._rep, self._evotype) + self.init_gpindices() + # NOTE: No _update_rep_base analogue is needed here. Each constituent + # operator's from_vector(...) updates its own attached OpRep, and the + # sum rep reflects those changes automatically. + + def submembers(self): + out = [] + hit = set() + for op in self._operators: + temp = op.submembers() + for sm in temp: + if id(temp) not in hit: + hit.add(id(temp)) + out.append(sm) + return out + + def deriv_wrt_params(self, wrt_filter=None): + raise NotImplementedError() + + def has_nonzero_hessian(self): + return any(op.has_nonzero_hessian() for op in self._operators) + + def from_vector(self, v, close=False, dirty_value=True): + for sm, local_inds in zip(self.submembers(), self._submember_rpindices): + sm.from_vector(v[local_inds], close, dirty_value) + + @property + def num_params(self): + return len(self.gpindices_as_array()) + + def to_vector(self): + v = _np.empty(self.num_params, 'd') + for param_op, local_inds in zip(self.submembers(), self._submember_rpindices): + v[local_inds] = param_op.to_vector() + return v + + def to_dense(self, on_space: SpaceT = 'HilbertSchmidt') -> _np.ndarray: + assert on_space in ('HilbertSchmidt', 'minimal') + out = self._operators[0].to_dense('HilbertSchmidt') + if not out.flags.writeable: + out = out.copy() + for op in self._operators[1:]: + out += op.to_dense('HilbertSchmidt') + return out diff --git a/pygsti/modelmembers/operations/denseop.py b/pygsti/modelmembers/operations/denseop.py index e5f600182..cce6e0088 100644 --- a/pygsti/modelmembers/operations/denseop.py +++ b/pygsti/modelmembers/operations/denseop.py @@ -255,7 +255,6 @@ def __complex__(self): return complex(self._ptr) class DenseOperator(DenseOperatorInterface, _KrausOperatorInterface, _LinearOperator): """ - TODO: update docstring An operator that behaves like a dense super-operator matrix. This class is the common base class for more specific dense operators. @@ -278,10 +277,29 @@ class DenseOperator(DenseOperatorInterface, _KrausOperatorInterface, _LinearOper The state space for this operation. If `None` a default state space with the appropriate number of qubits is used. - Attributes + Properties ---------- - base : numpy.ndarray - Direct access to the underlying process matrix data. + DenseOperatorInterface is deprecated and this class' dependence on it + will be removed in the immediate future. For the time being, the + implementation of __getattr__ in DenseOperatorInterface means that + DenseOperator has every property of its attached OpRepDenseSuperop, + which is stored in DenseOperator._rep. + + Private attributes + ------------------ + _evotype (required by ModelMember; positional arg in DenseOperator.__init__) + _rep (required by LinearOperator; either Cython or Python "OpRepDenseSuperop") + _basis (indirect requirement of _rep) + _ptr (required by DenseOperatorInterface, implemented as self._rep.base) + _state_space (required by ModelMember; inferrable from None as __init__ kwarg) + + Other private attributes + ------------------------ + ModelMember._gpindices + ModelMember._paramlbls + ModelMember._param_bounds + ModelMember._dirty + ModelMember._submember_rpindices """ @classmethod @@ -393,44 +411,9 @@ def _is_similar(self, other, rtol, atol): @property def kraus_operators(self): """A list of this operation's Kraus operators as numpy arrays.""" - # Let I index a basis element, rho be a d x d matrix, and (I // d, I mod d) := (i,ii) be the "2D" - # index corresponding to I. - # op(rho) = sum_IJ choi_IJ BI rho BJ_dag = sum_IJ (evecs_IK D_KK evecs_inv_KJ) BI rho BJ_dag - # Note: evecs can be and are assumed chosen to be orthonormal, so evecs_inv = evecs^dag - # if {Bi} is the set of matrix units then ... - # = sum_IJK(i',j') (evecs_IK D_KK evecs_inv_KJ) unitI_ii' rho_i'j' unitJ_j'j - # = sum_IJK(i',j') (evecs_(Ia,Ib)K D_KK evecs_inv_K(Ja,Jb)) unitI_ii' rho_i'j' unitJ_j'j - # using fact that sum(i') unitI_ii' ==> i'=Ib and delta_i,Ia factor - # = sum_IJK (evecs_(Ia,Ib)K D_KK evecs_inv_K(Ja,Jb)) rho_IbJa * delta_i,Ia, delta_j,Jb - # = sum_K D_KK [ sum_(Ib,Ja) evector[K]_(i,Ib) rho_IbJa dual_evector[K]_(Ja,j) ] - # -> let reshaped K-th evector be called O_K and dual O_K^dag - # = sum_K D_KK O_K rho O_K^dag assert(self._basis is not None), "Kraus operator functionality requires specifying a superoperator basis" - superop_mx = self.to_dense("HilbertSchmidt"); d = int(_np.round(_np.sqrt(superop_mx.shape[0]))) - std_basis = _Basis.cast('std', superop_mx.shape[0]) - choi_mx = _jt.jamiolkowski_iso(superop_mx, self._basis, std_basis) * d # see NOTE below - # NOTE: multiply by `d` (density mx dimension) to un-normalize choi_mx as given by - # jamiolkowski_iso. Here we *want* the trace of choi_mx to be `d`, not 1, so that - # op(rho) = sum_IJ choi_IJ BI rho BJ_dag is true. - - #CHECK 1 (to unit test?) REMOVE - #tmp_std = _bt.change_basis(superop_mx, self._basis, 'std') - #B = _bt.basis_matrices('std', superop_mx.shape[0]) - #check_superop = sum([ choi_mx[i,j] * _np.kron(B[i], B[j].conjugate()) for i in range(d*d) for j in range(d*d)]) - #assert(_np.allclose(check_superop, tmp_std)) - - evals, evecs = _np.linalg.eigh(choi_mx) - assert(_np.allclose(evecs @ _np.diag(evals) @ (evecs.conjugate().T), choi_mx)) - TOL = 1e-7 # consider lowering this tolerance as it leads to errors of this order in the Kraus decomp - if any([ev <= -TOL for ev in evals]): - raise ValueError("Cannot compute Kraus decomposition of non-positive-definite superoperator!") - kraus_ops = [evecs[:, i].reshape(d, d) * _np.sqrt(ev) for i, ev in enumerate(evals) if abs(ev) > TOL] - - #CHECK 2 (to unit test?) REMOVE - #std_superop = sum([_ot.unitary_to_std_process_mx(kop) for kop in kraus_ops]) - #assert(_np.allclose(std_superop, tmp_std)) - - return kraus_ops + kops = self._kraus_operators() + return kops def set_kraus_operators(self, kraus_operators): """ diff --git a/pygsti/modelmembers/operations/linearop.py b/pygsti/modelmembers/operations/linearop.py index 6edd84459..a1bc6695b 100644 --- a/pygsti/modelmembers/operations/linearop.py +++ b/pygsti/modelmembers/operations/linearop.py @@ -14,6 +14,7 @@ import numpy as _np +from pygsti.baseobjs.basis import Basis as _Basis 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 @@ -103,6 +104,16 @@ def shape(self) -> tuple[int, int]: # are given broad freedom to define semantics of self._rep. return (self.dim, self.dim) + def _kraus_operators(self): + """A list of this operation's Kraus operators as numpy arrays.""" + basis = getattr(self, '_basis', None) + if not isinstance(basis, _Basis): + msg = "Kraus operator functionality requires specifying a superoperator basis" + raise NotImplementedError(msg) + mx = self.to_dense('HilbertSchmidt') + kops = _ot.minimal_kraus_decomposition(mx, basis, 1e-7, 1e-7) + return kops + def set_dense(self, m): """ Set the dense-matrix value of this operation. diff --git a/pygsti/modelmembers/povms/__init__.py b/pygsti/modelmembers/povms/__init__.py index 8dab11ce4..a183cebf7 100644 --- a/pygsti/modelmembers/povms/__init__.py +++ b/pygsti/modelmembers/povms/__init__.py @@ -618,7 +618,7 @@ def _objfn(v): tol=1e-13) if not soln.success and soln.fun > 1e-6: # not "or" because success is often not set correctly raise ValueError("Failed to find an errorgen such that lab', conj_basis_mxs, Lbar, F1) # only a >= b nonzero (F1) - dVdp += _np.einsum('mal,mb,ab->lab', conj_basis_mxs, L, F1) # ditto - dVdp += _np.einsum('bml,ma,ab->lab', conj_basis_mxs, Lbar, F2) # only b > a nonzero (F2) - dVdp += _np.einsum('mbl,ma,ab->lab', conj_basis_mxs, L, F2.conjugate()) # ditto + dVdp = _np.einsum('aml,mb,ab->lab', conj_basis_mxs, Lbar, F1) # only a >= b nonzero (F1) + dVdp += _np.einsum('mal,mb,ab->lab', conj_basis_mxs, L, F1) # ditto + dVdp += _np.einsum('bml,ma,ab->lab', conj_basis_mxs, Lbar, F2) # only b > a nonzero (F2) + dVdp += _np.einsum('mbl,ma,ab->lab', conj_basis_mxs, L, F2.conj()) # ditto dVdp.shape = [dVdp.shape[0], nP] # jacobian with respect to "p" params, # which don't include normalization for TP-constraint @@ -363,28 +363,3 @@ def has_nonzero_hessian(self): bool """ return True - - def hessian_wrt_params(self, wrt_filter1=None, wrt_filter2=None): - """ - Construct the Hessian of this state vector with respect to its parameters. - - This function returns a tensor whose first axis corresponds to the - flattened operation matrix and whose 2nd and 3rd axes correspond to the - parameters that are differentiated with respect to. - - Parameters - ---------- - wrt_filter1 : list or numpy.ndarray - List of parameter indices to take 1st derivatives with respect to. - (None means to use all the this operation's parameters.) - - wrt_filter2 : list or numpy.ndarray - List of parameter indices to take 2nd derivatives with respect to. - (None means to use all the this operation's parameters.) - - Returns - ------- - numpy array - Hessian with shape (dimension, num_params1, num_params2) - """ - raise NotImplementedError("TODO: add hessian computation for CPTPState") diff --git a/pygsti/modelpacks/_modelpack.py b/pygsti/modelpacks/_modelpack.py index 796f4a5df..5c8fc29b0 100644 --- a/pygsti/modelpacks/_modelpack.py +++ b/pygsti/modelpacks/_modelpack.py @@ -360,7 +360,8 @@ def create_gst_experiment_design(self, max_max_length, qubit_labels=None, fpr=Fa """ for k in kwargs.keys(): if k not in ('germ_length_limits', 'keep_fraction', 'keep_seed', 'include_lgst', 'nest', 'circuit_rules', - 'op_label_aliases', 'dscheck', 'action_if_missing', 'verbosity', 'add_default_protocol'): + 'op_label_aliases', 'dscheck', 'action_if_missing', 'verbosity', 'add_default_protocol', + 'germs'): raise ValueError("Invalid argument '%s' to StandardGSTDesign constructor" % k) if qubit_labels is None: qubit_labels = self._sslbls @@ -380,11 +381,13 @@ def create_gst_experiment_design(self, max_max_length, qubit_labels=None, fpr=Fa else: max_lengths_list = list(_gen_max_length(max_max_length)) + germs = kwargs.get('germs', self.germs(qubit_labels, lite)) + return _gst.StandardGSTDesign( self.processor_spec(qubit_labels), self.prep_fiducials(qubit_labels), self.meas_fiducials(qubit_labels), - self.germs(qubit_labels, lite), + germs, max_lengths_list, kwargs.get('germ_length_limits', None), fidpairs, diff --git a/pygsti/models/model.py b/pygsti/models/model.py index 2b4be7278..7fd3f474c 100644 --- a/pygsti/models/model.py +++ b/pygsti/models/model.py @@ -1099,7 +1099,6 @@ def _get_shift(j): return _bisect.bisect_left(indices_to_remove, j) #rebuild the model index to model member map if needed. self._build_index_mm_map() - def _init_virtual_obj(self, obj): """ Initializes a "virtual object" - an object (e.g. LinearOperator) that *could* be a @@ -1224,8 +1223,6 @@ def set_parameter_value(self, index, val, close=False): """ self.set_parameter_values([index], [val], close) - - def set_parameter_values(self, indices, values, close=False): """ diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py index d43ac70f8..85a20e347 100644 --- a/pygsti/tools/matrixtools.py +++ b/pygsti/tools/matrixtools.py @@ -14,7 +14,7 @@ import itertools as _itertools import warnings as _warnings -from typing import Protocol, Any, runtime_checkable, TypeVar, Optional, Union +from typing import Protocol, Any, runtime_checkable, TypeVar, Optional, Union, Literal import numpy as _np import scipy.linalg as _spl @@ -150,6 +150,54 @@ def is_valid_density_mx(mx, tol=1e-9): return abs(_np.trace(mx) - 1.0) < tol and is_pos_def(mx, tol) +def project_onto_simplex(vec: _np.ndarray, ord:Literal[1,2,'inf']) -> _np.ndarray: + """ + Return the solution, `proj`, to the optimization problem + + min{ ||vec - proj||_ord : 0 <= proj, sum(proj) == 1 }, + + where ||*||_ord refers to the standard vector `ord`-norm. + """ + vec = vec.ravel() + if ord == 2: + # Sort in descending order + u = _np.sort(vec)[::-1] + css = _np.cumsum(u) + # Find the largest index rho such that + # u[rho] > (css[rho] - 1) / (rho + 1) + inds = _np.arange(vec.size) + 1 + cond = u - (css - 1) / inds > 0 + if not _np.any(cond): + # fallback: if no positive entries, project uniformly on first coordinate + theta = (css[0] - 1) / 1 + else: + rho = _np.where(cond)[0][-1] + theta = (css[rho] - 1) / (rho + 1) + # The projection is max(flat - theta, 0) + proj = _np.maximum(vec - theta, 0) + else: + raise NotImplementedError() + return proj + + +def project_onto_densities(herm_mx: _np.ndarray, ord:Literal['fro','nuc']='fro') -> _np.ndarray: + if ord != 'fro': + raise NotImplementedError() + assert_hermitian(herm_mx, 1e-8) + u, v, udag = eigendecomposition(herm_mx, assume_hermitian=True) + proj_v = project_onto_simplex(v, 2) + density_mx = u @ (proj_v[:, None] * udag) + return density_mx + + +def clip_eigenvalues(herm_mx: _np.ndarray, lower: float, upper:float) -> _np.ndarray: + assert_hermitian(herm_mx, 1e-8) + u, v, udag = eigendecomposition(herm_mx, assume_hermitian=True) + proj_v = _np.clip(v, lower, upper) + clipped_mx = u @ (proj_v[:, None] * udag) + return clipped_mx + + def nullspace(m, tol=1e-7): """ Compute the nullspace of a matrix. diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 4551375cf..738c592ab 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -387,7 +387,6 @@ def diamonddist(a, b, mx_basis='pp', return_x=False): `dm = trace( |(J(a)-J(b)).T * W| )`. """ assert _sdps.CVXPY_ENABLED - import cvxpy as _cp assert a.shape == b.shape mx_basis = _bt.create_basis_for_matrix(a, mx_basis) @@ -396,33 +395,22 @@ def diamonddist(a, b, mx_basis='pp', return_x=False): # It will convert c a "single-block" basis representation # when mx_basis has multiple blocks. J = _jam.fast_jamiolkowski_iso_std(c, mx_basis, normalized=False) - prob, vars = _sdps.diamond_norm_model_jamiolkowski(J) + problem, vars = _sdps.diamond_norm_model_jamiolkowski(J) + objval, varvals = _sdps.solve_sdp(problem) - objective_val = -2 - varvals = [_np.zeros_like(J), None, None] - sdp_solvers = ['MOSEK', 'CLARABEL', 'CVXOPT'] - for i, solver in enumerate(sdp_solvers): - try: - with _warnings.catch_warnings(): - _warnings.filterwarnings('ignore','.*Solution may be inaccurate.*', UserWarning) - prob.solve(solver=solver) - objective_val = prob.value - varvals = [v.value for v in vars] - break - except (AssertionError, _cp.SolverError) as e: - if solver != 'MOSEK': - msg = f"Received error {e} when trying to use solver={solver}." - if i + 1 == len(sdp_solvers): - failure_msg = "Out of solvers. Returning -2 for diamonddist." - else: - failure_msg = f"Trying {sdp_solvers[i+1]} next." - msg += f'\n{failure_msg}' - _warnings.warn(msg, _CVXPYFailure) + if _np.isnan(objval): + objval = -2 + varvals = [_np.zeros_like(J), None, None] + else: + # the returned `varvals` is a dict keyed by + # variable name. We care about the variables + # as ordered by `vars`. + varvals = [v.value for v in vars] if return_x: - return objective_val, varvals + return objval, varvals else: - return objective_val + return objval def jtracedist(a, b, mx_basis='pp'): # Jamiolkowski trace distance: Tr(|J(a)-J(b)|) @@ -478,6 +466,15 @@ def is_trace_preserving(a: _np.ndarray, mx_basis: BasisLike='pp', tol=1e-8) -> b return check +def superket_trace(superket: _np.ndarray, basis: _Basis): + if basis.first_element_is_identity: + t = superket.ravel()[0] + else: + mx = _bt.vec_to_stdmx(superket, basis, keep_complex=True) + t = _np.real(_np.trace(mx)) + return t + + def entanglement_fidelity(a, b, mx_basis: BasisLike='pp', is_tp=None, is_unitary=None): """ Returns the "entanglement" process fidelity between gate matrices. @@ -591,6 +588,117 @@ def tensorized_with_eye(op: _np.ndarray, op_basis: _Basis, ten_basis: Optional[_ return ten_op, ten_basis +def rootconj_superop(effect_superket: _np.ndarray, basis: _Basis, abstol_warn: float=1e-10, abstol_error: float=1e-8) -> _np.ndarray: + """ + Let E denote the Hermitian matrix representation effect_superket, where 0 ≤ E ≤ 1. + + This function returns the array representation (in `basis`) of the map that takes + a Hermitian matrix ρ to the Hermitian matrix E½ ρ E½. + """ + effect_mat = _bt.vec_to_stdmx(effect_superket, basis, keep_complex=True) + vecs, vals, inv_vecs = _mt.eigendecomposition(effect_mat, assume_hermitian=True) + + msg = f'Eigenvalues {vals} fall outside [0.0, 1.0], up to tolerance %s.' + if _np.any(vals < 0.0 - abstol_error) or _np.any(vals > 1.0 + abstol_error): + raise ValueError( msg % abstol_error ) + if _np.any(vals < 0.0 - abstol_warn) or _np.any(vals > 1.0 + abstol_warn ): + _warnings.warn( msg % abstol_warn ) + vals = _np.clip(vals, a_min=0.0, a_max=1.0) + + rooteffect_mat = (vecs * _np.sqrt(vals)[_np.newaxis, :]) @ inv_vecs + mx_std = _np.kron(rooteffect_mat, rooteffect_mat.T) + mx : _np.ndarray = _bt.change_basis(mx_std, 'std', basis, expect_real=True) # type: ignore + return mx + + +def partial_trace(mx: _np.ndarray, dims: tuple[int,...], axis: int) -> _np.ndarray: + """ + + Notes + ----- + This implementation is stolen from CVXPY, which was stolen from some Julia library. + """ + if mx.ndim < 2 or mx.shape[0] != mx.shape[1]: + raise ValueError("partial_trace only supports 2-d square arrays.") + if axis < 0 or axis >= len(dims): + msg = f"Invalid axis argument, should be between 0 and {len(dims)}, got {axis}." + raise ValueError(msg) + if mx.shape[0] != _np.prod(dims): + raise ValueError("Dimension of system doesn't correspond to dimension of subsystems.") + + def _term(j: int) -> _np.ndarray: + a = _sps.coo_matrix(([1.0], ([0], [0]))) + b = _sps.coo_matrix(([1.0], ([0], [0]))) + for (i_axis, dim) in enumerate(dims): + if i_axis == axis: + v = _sps.coo_matrix(([1], ([j], [0])), shape=(dim, 1)) + a = _sps.kron(a, v.T) + b = _sps.kron(b, v) + else: + eye_mat = _sps.eye_array(dim) + a = _sps.kron(a, eye_mat) + b = _sps.kron(b, eye_mat) + return a @ mx @ b + + return _np.sum([_term(j) for j in range(dims[axis])]) + + +def trace_effect(op: _np.ndarray, op_basis: BasisLike, on_space: SpaceT = 'HilbertSchmidt'): + """ + Let `op` be the array representation of a superoperator G in `op_basis`, + where G maps from and to the space of order-d Hermitian operators. + + The trace effect of G is the Heritian operator E that satifies + + trace(G(ρ)) = trace(E ρ) for all order-d Hermitian matrices ρ. + + If on_space='HilbertSchmidt', then this function returns a superket representation + of E in `op_basis`. If on_space='Hilbert', then we return E itself. + """ + d = int(_np.round(op.size ** 0.25)) + assert op.shape == (d**2, d**2) + basis = op_basis if isinstance(op_basis, _Basis) else _Basis.cast(op_basis, dim=d**2) + vecI = _bt.stdmx_to_vec(_np.eye(d), basis) + vecE = op.T.conj() @ vecI + if on_space == 'HilbertSchmidt': + return vecE + else: + E = _bt.vec_to_stdmx(vecE, op_basis) + return E + + +def minimal_kraus_decomposition(op_x: _np.ndarray, op_basis: _Basis, error_tol:float=1e-6, trunc_tol:float=1e-7) -> list[_np.ndarray]: + """ + The array `op_x` represents a completely positive superoperator X on + Hilbert-Schmidt space, using `op_basis` as the basis for that space. + + A Kraus decomposition of X is a set of square matrices, {K_i}_i, that satisfy + + X(ρ) = \\sum_i K_i ρ K_i^\\dagger. + + The matrices appearing in any such set are called Kraus operators of X. + + This function returns a minimal-length list of Kraus operators of X. + """ + d = int(_np.round(op_x.size ** 0.25)) + d2 = d**2 + assert op_x.shape == (d2, d2) + from pygsti.tools.jamiolkowski import jamiolkowski_iso + choi_mx : _np.ndarray = jamiolkowski_iso(op_x, op_basis, 'std', normalized=True) * d # type: ignore + + evecs, evals, _ = _mt.eigendecomposition(choi_mx, assume_hermitian=True) + if any([ev < -error_tol for ev in evals]): + raise ValueError("Cannot compute Kraus decomposition of non-positive-definite superoperator!") + keep = evals >= trunc_tol + evals = evals[keep] + evecs = evecs[:, keep] + out = [] + for i, ev in enumerate(evals): + temp = _np.sqrt(ev) * evecs[:, i].reshape((d, d), order='F') + out.append(temp) + return out + + def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): """ Computes the average gate fidelity (AGF) between two gates. @@ -1111,7 +1219,7 @@ def povm_diamonddist(model, target_model, povmlbl, _premultiplier=None): return diamonddist(povm_mx, target_povm_mx, target_model.basis) except AssertionError as e: assert '`dim` must be a perfect square' in str(e) - return _np.NaN + return _np.nan def instrument_infidelity(a, b, mx_basis): diff --git a/pygsti/tools/sdptools.py b/pygsti/tools/sdptools.py index af9338364..dd3e05435 100644 --- a/pygsti/tools/sdptools.py +++ b/pygsti/tools/sdptools.py @@ -13,13 +13,15 @@ from __future__ import annotations import numpy as np +import warnings -from typing import Union, List, Tuple, TYPE_CHECKING +from typing import Any, Union, List, Tuple, Sequence, TYPE_CHECKING if TYPE_CHECKING: import cvxpy as cp ExpressionLike = Union[cp.Expression, np.ndarray] -from pygsti.baseobjs.basis import Basis, BasisLike +from pygsti.baseobjs import Basis, BasisLike +from pygsti.tools.matrixtools import assert_hermitian from pygsti.tools.basistools import stdmx_to_vec from pygsti.tools.jamiolkowski import jamiolkowski_iso @@ -34,6 +36,70 @@ CVXPY_ENABLED = False +SDP_SOLVER_PRIORITY = ['MOSEK', 'CLARABEL', 'CVXOPT'] + + +def solve_sdp(prob: cp.Problem, **kwargs) -> tuple[np.floating, dict[str, np.ndarray]]: + + objective_val : np.floating = np.array(np.nan).item() + varvals : dict[str, np.ndarray] = dict() + for i, solver in enumerate(SDP_SOLVER_PRIORITY): + try: + with warnings.catch_warnings(): + warnings.filterwarnings('ignore','.*Solution may be inaccurate.*', UserWarning) + prob.solve(solver=solver, **kwargs) + objective_val = prob.value # type: ignore + varvals.update({k: v.value for (k, v) in prob.var_dict.items()}) # type: ignore + break + except (AssertionError, cp.SolverError) as e: + if solver != 'MOSEK': + msg = f"Received error {e} when trying to use solver={solver}." + if i + 1 == len(SDP_SOLVER_PRIORITY): + failure_msg = "Out of solvers. Returning NaN." + else: + failure_msg = f"Trying {SDP_SOLVER_PRIORITY[i+1]} next." + msg += f'\n{failure_msg}' + warnings.warn(msg) + + return objective_val, varvals + + +def povm_projection_model(effects: Sequence[np.ndarray], independent_complement: bool=False) -> cp.Problem: + N = round(effects[0].size ** 0.5) + cumsum_E = np.zeros((N, N), dtype=complex) + for E in effects: + assert_hermitian(E, tol=1e-8) + cumsum_E += E + if np.real(np.trace(cumsum_E)) < (N - 0.01): + # we're going to assume the complement POVM effect was omitted. + complement_E = np.eye(N) - np.sum(effects, axis=0) + effects = [E for E in effects] # shallow-copy + effects.append( complement_E ) + + proj_effects, constraints = povm_effect_variables(N, len(effects), independent_complement) + objective_expr = cp.Constant(0.0) + for E_proj, E in zip( proj_effects, effects ): + objective_expr = objective_expr + cp.sum_squares(E_proj - E) + objective = cp.Minimize(objective_expr) + + prob = cp.Problem(objective, constraints) + return prob + + +def povm_effect_variables(N: int, num_effects: int, independent_complement: bool) -> tuple[list[cp.Expression], list[cp.Constraint]]: + num_vars = num_effects if independent_complement else num_effects - 1 + vars : list[cp.Expression] = [cp.Variable(shape=(N, N), hermitian=True) for _ in range(num_vars)] # type: ignore + cons : list[cp.Constraint] = [v >> 0 for v in vars] + if num_vars == num_effects: + expr : cp.Expression = cp.sum(vars) # type: ignore + cons.append( expr == np.eye(N) ) + else: + expr : cp.Expression = np.eye(N) - cp.sum(vars) # type: ignore + cons.append( expr >> 0 ) + vars.append( expr ) + return vars, cons + + def diamond_norm_model_jamiolkowski(J: ExpressionLike) -> tuple[cp.Problem, List[cp.Variable]]: # return a model for computing the diamond norm. # diff --git a/test/unit/modelmembers/test_instrument.py b/test/unit/modelmembers/test_instrument.py new file mode 100644 index 000000000..2ae9f4329 --- /dev/null +++ b/test/unit/modelmembers/test_instrument.py @@ -0,0 +1,222 @@ +import numpy as np + +import pygsti.modelmembers.operations as op +from pygsti.modelmembers import instruments as inst +from pygsti.modelmembers.instruments import cptp_instrument +from pygsti.modelmembers.instruments import Instrument, TPInstrument +from pygsti.modelmembers.operations import RootConjOperator, SummedOperator +from pygsti.modelmembers import povms as pv +from pygsti.modelpacks.legacy import std1Q_XYI as std +from pygsti.models.gaugegroup import FullGaugeGroupElement +from ..util import BaseCase +from .test_operation import ImmutableDenseOpBase + + +class InstrumentTestBase(BaseCase): + """Shared setUp and test methods for Instrument and TPInstrument. + + Subclasses must define a class attribute ``constructor`` pointing to the + instrument class under test. + """ + __test__ = False + constructor: type + + def setUp(self): + self.n_elements = 32 + self.model = std.target_model() + E = self.model.povms['Mdefault']['0'] + Erem = self.model.povms['Mdefault']['1'] + self.Gmz_plus = np.dot(E, E.T) + self.Gmz_minus = np.dot(Erem, Erem.T) + self.instrument: Instrument = self.constructor({'plus': self.Gmz_plus, 'minus': self.Gmz_minus}) + self.model.instruments['Iz'] = self.instrument + + def test_num_elements(self): + self.assertEqual(self.instrument.num_elements, self.n_elements) + + def test_copy(self): + inst_copy = self.instrument.copy() + self.assertIsInstance(inst_copy, type(self.instrument)) + self.assertEqual(list(inst_copy.keys()), list(self.instrument.keys())) + for key in self.instrument.keys(): + actual = inst_copy[key].to_dense() + expected = self.instrument[key].to_dense() + self.assertArraysEqual(actual, expected) + + def test_to_string(self): + inst_str = str(self.instrument) + self.assertIsInstance(inst_str, str) + self.assertIn("Instrument with elements:", inst_str) + for key in self.instrument.keys(): + self.assertIn(str(key), inst_str) + + def test_transform(self): + T = FullGaugeGroupElement( + np.array([[1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 0, 1], + [0, 0, 1, 0]], 'd')) + T_mx = T.transform_matrix + T_inv = T.transform_matrix_inverse + originals = {k: v.to_dense().copy() for k, v in self.instrument.items()} + self.instrument.transform_inplace(T) + for key, E_orig in originals.items(): + expected = T_mx @ E_orig @ T_inv + self.assertArraysAlmostEqual(self.instrument[key].to_dense(), expected) + + def test_simplify_operations(self): + gates = self.instrument.simplify_operations(prefix="ABC") + self.assertEqual(len(gates), len(self.instrument)) + expected_keys = ["ABC_" + k for k in self.instrument.keys()] + self.assertEqual(list(gates.keys()), expected_keys) + for gate, orig in zip(gates.values(), self.instrument.values()): + self.assertIs(gate, orig) + + def test_convert_raises_on_unknown_parameterization(self): + with self.assertRaises(ValueError): + inst.convert(self.instrument, "H+S", self.model.basis) + + +class InstrumentInstanceTester(InstrumentTestBase): + __test__ = True + constructor = inst.Instrument + + +class TPInstrumentInstanceTester(InstrumentTestBase): + __test__ = True + constructor = inst.TPInstrument + + def test_raise_on_modify(self): + with self.assertRaises(ValueError): + self.instrument['plus'] = None + + +class CPTPInstrumentTester(BaseCase): + """Tests for cptp_instrument and the CPTR operation primitives.""" + + def setUp(self): + model = std.target_model() + E = model.povms['Mdefault']['0'] + Erem = model.povms['Mdefault']['1'] + self.Gmz_plus = np.dot(E, E.T) + self.Gmz_minus = np.dot(Erem, Erem.T) + self.basis = model.basis + self.op_arrays = {'plus': self.Gmz_plus, 'minus': self.Gmz_minus} + self.instrument = cptp_instrument(self.op_arrays, self.basis) + + # --- cptp_instrument factory tests --- + + def test_construction_returns_instrument(self): + self.assertIsInstance(self.instrument, Instrument) + + def test_construction_keys(self): + self.assertSetEqual(set(self.instrument.keys()), {'plus', 'minus'}) + + def test_to_from_vector_roundtrip(self): + v0 = self.instrument.to_vector() + self.assertEqual(len(v0), self.instrument.num_params) + # Perturb and restore; dense output should change then return. + v_perturbed = v0 + 1e-3 * np.random.default_rng(0).standard_normal(len(v0)) + self.instrument.from_vector(v_perturbed) + v1 = self.instrument.to_vector() + self.assertFalse(np.allclose(v0, v1)) + self.instrument.from_vector(v0) + v2 = self.instrument.to_vector() + np.testing.assert_allclose(v0, v2, atol=1e-12) + + def test_convert_to_cptplnd(self): + plain = inst.Instrument({'plus': self.Gmz_plus, 'minus': self.Gmz_minus}) + converted = inst.convert(plain, 'CPTPLND', self.basis) + self.assertIsInstance(converted, Instrument) + self.assertSetEqual(set(converted.keys()), {'plus', 'minus'}) + + def test_member_dense_shapes(self): + dim = self.basis.dim + for op in self.instrument.values(): + mx = op.to_dense('HilbertSchmidt') + self.assertEqual(mx.shape, (dim, dim)) + + def test_sum_of_members_is_cptp_like(self): + # The sum of all instrument members' superoperators should equal the + # superoperator of the total channel (which is CPTP), i.e. the first + # row of the sum must be [1, 0, 0, 0] in the pp basis. + total = sum(op.to_dense('HilbertSchmidt') for op in self.instrument.values()) + # First row encodes trace preservation: tr(rho) is preserved. + np.testing.assert_allclose(total[0, 0], 1.0, atol=1e-6) + np.testing.assert_allclose(total[0, 1:], 0.0, atol=1e-6) + + # --- RootConjOperator unit tests --- + + def test_root_conj_operator_to_dense_shape(self): + # Build a real StaticPOVMEffect from the plus POVM effect + effect = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + rcop = RootConjOperator(effect, self.basis) + mx = rcop.to_dense('HilbertSchmidt') + dim = self.basis.dim + self.assertEqual(mx.shape, (dim, dim)) + + def test_root_conj_operator_num_params(self): + effect = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + rcop = RootConjOperator(effect, self.basis) + # StaticPOVMEffect has zero parameters; RootConjOperator inherits them. + self.assertEqual(rcop.num_params, 0) + + def test_root_conj_operator_has_nonzero_hessian(self): + effect = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + rcop = RootConjOperator(effect, self.basis) + self.assertTrue(rcop.has_nonzero_hessian()) + + # --- SummedOperator unit tests --- + + def test_summed_operator_to_dense_equals_sum(self): + e1 = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + e2 = pv.StaticPOVMEffect(self.Gmz_minus[:, 0]) + r1 = RootConjOperator(e1, self.basis) + r2 = RootConjOperator(e2, self.basis) + sop = SummedOperator([r1, r2], self.basis) + expected = r1.to_dense() + r2.to_dense() + np.testing.assert_allclose(sop.to_dense(), expected, atol=1e-12) + + def test_summed_operator_num_params(self): + e1 = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + e2 = pv.StaticPOVMEffect(self.Gmz_minus[:, 0]) + r1 = RootConjOperator(e1, self.basis) + r2 = RootConjOperator(e2, self.basis) + sop = SummedOperator([r1, r2], self.basis) + self.assertEqual(sop.num_params, 0) + + def test_summed_operator_deriv_raises(self): + e1 = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + r1 = RootConjOperator(e1, self.basis) + sop = SummedOperator([r1], self.basis) + with self.assertRaises(NotImplementedError): + sop.deriv_wrt_params() + + +class TPInstrumentOpTester(ImmutableDenseOpBase, BaseCase): + n_params = 28 + + @staticmethod + def build_gate(): + Gmz_plus = np.array([[0.5, 0, 0, 0.5], + [0, 0, 0, 0], + [0, 0, 0, 0], + [0.5, 0, 0, 0.5]]) + Gmz_minus = np.array([[0.5, 0, 0, -0.5], + [0, 0, 0, 0], + [0, 0, 0, 0], + [-0.5, 0, 0, 0.5]]) + evotype = 'default' + instrument = TPInstrument({ + 'plus': op.FullArbitraryOp(Gmz_plus, 'pp', evotype), + 'minus': op.FullArbitraryOp(Gmz_minus, 'pp', evotype) + }) + return instrument['plus'] + + def test_vector_conversion(self): + self.gate.to_vector() # now to_vector is allowed + + def test_deriv_wrt_params_shape(self): + super(TPInstrumentOpTester, self).test_deriv_wrt_params() + deriv = self.gate.deriv_wrt_params([0]) + self.assertEqual(deriv.shape[1], 1) diff --git a/test/unit/modelmembers/test_operation.py b/test/unit/modelmembers/test_operation.py index a9861d951..a0990eedb 100644 --- a/test/unit/modelmembers/test_operation.py +++ b/test/unit/modelmembers/test_operation.py @@ -10,7 +10,6 @@ import pygsti.tools.optools as gt from pygsti.models.modelconstruction import create_spam_vector, create_operation from pygsti.evotypes import Evotype -from pygsti.modelmembers.instruments import TPInstrument from pygsti.modelmembers.states import FullState from pygsti.models import ExplicitOpModel from pygsti.baseobjs import statespace, basisconstructors as bc @@ -219,6 +218,7 @@ def test_rotate(self): self.gate.rotate([0.01, 0.02, 0.03], 'gm') # TODO assert correctness + class ImmutableDenseOpBase(DenseOpBase): def test_raises_on_set_value(self): M = np.asarray(self.gate) # gate as a matrix @@ -430,6 +430,7 @@ def test_first_row_read_only(self): with self.assertRaises(ValueError): self.gate[0][1:2] = [0] + class AffineShiftOpTester(DenseOpBase, BaseCase): n_params = 3 @@ -773,36 +774,6 @@ def test_constructor_raises_on_state_space_label_mismatch(self): op.EmbeddedOp(state_space, ['Q0', 'Q1'], op.FullArbitraryOp(mx, evotype=evotype, state_space=None)) -class TPInstrumentOpTester(ImmutableDenseOpBase, BaseCase): - n_params = 28 - - @staticmethod - def build_gate(): - # XXX can this be constructed directly? EGN: what do you mean? - Gmz_plus = np.array([[0.5, 0, 0, 0.5], - [0, 0, 0, 0], - [0, 0, 0, 0], - [0.5, 0, 0, 0.5]]) - Gmz_minus = np.array([[0.5, 0, 0, -0.5], - [0, 0, 0, 0], - [0, 0, 0, 0], - [-0.5, 0, 0, 0.5]]) - evotype = 'default' - inst = TPInstrument({'plus': op.FullArbitraryOp(Gmz_plus, evotype=evotype), 'minus': op.FullArbitraryOp( - Gmz_minus, evotype=evotype)}) - return inst['plus'] - - def test_vector_conversion(self): - self.gate.to_vector() # now to_vector is allowed - - def test_deriv_wrt_params(self): - super(TPInstrumentOpTester, self).test_deriv_wrt_params() - - # XXX does this check anything meaningful? EGN: yes, this checks that when I give deriv_wrt_params a length-1 list it's return value has the right shape. - deriv = self.gate.deriv_wrt_params([0]) - self.assertEqual(deriv.shape[1], 1) - - class StochasticNoiseOpTester(BaseCase): def test_instance(self): state_space = statespace.default_space_for_dim(4) diff --git a/test/unit/objects/test_instrument.py b/test/unit/objects/test_instrument.py deleted file mode 100644 index b11850a5c..000000000 --- a/test/unit/objects/test_instrument.py +++ /dev/null @@ -1,74 +0,0 @@ -import numpy as np - -from pygsti.modelmembers import instruments as inst -from pygsti.modelpacks.legacy import std1Q_XYI as std -from pygsti.models.gaugegroup import FullGaugeGroupElement -from ..util import BaseCase - - -class InstrumentMethodBase(object): - def test_num_elements(self): - self.assertEqual(self.instrument.num_elements, self.n_elements) - - def test_copy(self): - inst_copy = self.instrument.copy() - # TODO assert correctness - - def test_to_string(self): - inst_str = str(self.instrument) - # TODO assert correctness - - def test_transform(self): - T = FullGaugeGroupElement( - np.array([[1, 0, 0, 0], - [0, 1, 0, 0], - [0, 0, 0, 1], - [0, 0, 1, 0]], 'd')) - self.instrument.transform_inplace(T) - # TODO assert correctness - - def test_simplify_operations(self): - gates = self.instrument.simplify_operations(prefix="ABC") - # TODO assert correctness - - def test_constructor_raises_on_non_none_param_conflict(self): - with self.assertRaises(AssertionError): - self.constructor(["Non-none-matrices"], 'default', None, False, ["Non-none-items"]) # can't both be non-None - - def test_constructor_raises_on_bad_op_matrices_type(self): - with self.assertRaises(ValueError): - self.constructor("foobar") # op_matrices must be a list or dict - - def test_convert_raises_on_unknown_basis(self): - with self.assertRaises(ValueError): - inst.convert(self.instrument, "foobar", self.model.basis) - - -class InstrumentInstanceBase(object): - def setUp(self): - # Initialize standard target model for instruments - # XXX can instruments be tested independently of a model? EGN: yes, I was just lazy; but they should also be tested within a model. - self.n_elements = 32 - - self.model = std.target_model() - E = self.model.povms['Mdefault']['0'] - Erem = self.model.povms['Mdefault']['1'] - self.Gmz_plus = np.dot(E, E.T) - self.Gmz_minus = np.dot(Erem, Erem.T) - # XXX is this used? - self.povm_ident = self.model.povms['Mdefault']['0'] + self.model.povms['Mdefault']['1'] - self.instrument = self.constructor({'plus': self.Gmz_plus, 'minus': self.Gmz_minus}) - self.model.instruments['Iz'] = self.instrument - super(InstrumentInstanceBase, self).setUp() - - -class InstrumentInstanceTester(InstrumentMethodBase, InstrumentInstanceBase, BaseCase): - constructor = inst.Instrument - - -class TPInstrumentInstanceTester(InstrumentMethodBase, InstrumentInstanceBase, BaseCase): - constructor = inst.TPInstrument - - def test_raise_on_modify(self): - with self.assertRaises(ValueError): - self.instrument['plus'] = None # can't set value of a TP Instrument element