Ensure Consistency of Global Parameter Indices.

pygsti/models/model.py

@@ -665,21 +665,25 @@ def _iter_parameterized_objs(self):
     def _check_paramvec(self, debug=False):
         if debug: print("---- Model._check_paramvec ----")
 
+        ops_paramvec = self._model_paramvec_to_ops_paramvec(self._paramvec)
+        if debug: 
+            print(f'{ops_paramvec=}')
+            print(f'{self._paramvec=}')
         TOL = 1e-8
         for lbl, obj in self._iter_parameterized_objs():
             if debug: print(lbl, ":", obj.num_params, obj.gpindices)
             w = obj.to_vector()
             msg = "None" if (obj.parent is None) else id(obj.parent)
             assert(obj.parent is self), "%s's parent is not set correctly (%s)!" % (lbl, msg)
             if obj.gpindices is not None and len(w) > 0:
-                if _np.linalg.norm(self._paramvec[obj.gpindices] - w) > TOL:
+                if _np.linalg.norm(ops_paramvec[obj.gpindices] - w) > TOL:
                     if debug: print(lbl, ".to_vector() = ", w, " but Model's paramvec = ",
-                                    self._paramvec[obj.gpindices])
-                    raise ValueError("%s is out of sync with paramvec!!!" % lbl)
+                                    ops_paramvec[obj.gpindices])
+                    raise ValueError(f"{str(lbl)} is out of sync with paramvec!!!")
             if not self.dirty and obj.dirty:
-                raise ValueError("%s is dirty but Model.dirty=False!!" % lbl)
+                raise ValueError(f"{str(lbl)} is dirty but Model.dirty=False!!")
 
-    def _clean_paramvec(self):
+    def _clean_paramvec(self, debug=False):
         """ Updates _paramvec corresponding to any "dirty" elements, which may
             have been modified without out knowing, leaving _paramvec out of
             sync with the element's internal data.  It *may* be necessary
@@ -708,7 +712,6 @@ def _clean_paramvec(self):
         if self.dirty:  # if any member object is dirty (ModelMember.dirty setter should set this value)
             TOL = 1e-8
             ops_paramvec = self._model_paramvec_to_ops_paramvec(self._paramvec)
-
             #Note: lbl args used *just* for potential debugging - could strip out once
             # we're confident this code always works.
             def clean_single_obj(obj, lbl):  # sync an object's to_vector result w/_paramvec
@@ -723,7 +726,8 @@ def clean_single_obj(obj, lbl):  # sync an object's to_vector result w/_paramvec
                         else:
                             raise e  # we don't know what went wrong.
                     chk_norm = _np.linalg.norm(ops_paramvec[obj.gpindices] - w)
-                    #print(lbl, " is dirty! vec = ", w, "  chk_norm = ",chk_norm)
+                    if debug:
+                        print(f"{lbl} is dirty! vec = {w}, chk_norm = {chk_norm} gpindices = {obj.gpindices}")
                     if (not _np.isfinite(chk_norm)) or chk_norm > TOL:
                         ops_paramvec[obj.gpindices] = w
                     obj.dirty = False
@@ -984,9 +988,10 @@ def _rebuild_paramvec(self):
                 max_index_processed_so_far = max(max_index_processed_so_far, max_existing_index)
                 insertion_point = max_index_processed_so_far + 1
                 if num_new_params > 0:
+                    memo = set()
                     # If so, before allocating anything, make the necessary space in the parameter arrays:
                     for _, o in self._iter_parameterized_objs():
-                        o.shift_gpindices(insertion_point, num_new_params, self)
+                        o.shift_gpindices(insertion_point, num_new_params, self, memo)
                     w = _np.insert(w, insertion_point, _np.empty(num_new_params, 'd'))
                     wl = _np.insert(wl, insertion_point, _np.empty(num_new_params, dtype=object))
                     wb = _np.insert(wb, insertion_point, _default_param_bounds(num_new_params), axis=0)
@@ -1327,7 +1332,13 @@ def param_interposer(self, interposer):
             self._paramvec = self._model_paramvec_to_ops_paramvec(self._paramvec)
         self._param_interposer = interposer
         if interposer is not None:  # add new interposer
+            #the clean_paramvec call below happens at a point between when we've set the interposer 
+            #attribute, but before we've used it to set the model's parameter vector, which will results
+            #in an edge case in the parameter vector integrity logic, so disable temporarily.
+            pcheck = self._pcheck
+            OpModel._pcheck = False
             self._clean_paramvec()
+            OpModel._pcheck = pcheck
             self._paramvec = self._ops_paramvec_to_model_paramvec(self._paramvec)
 
     def _model_paramvec_to_ops_paramvec(self, v):

test/unit/extras/interpygate/test_construction.py

@@ -166,7 +166,6 @@ def setUpClass(cls):
         cls.model['Gxpi2',0] = x_gate
         cls.model['Gypi2',0] = y_gate
 
-
     def test_gpindices(self):
         model = self.model.copy()
         model['rho0'].set_gpindices(slice(0,4),model)
@@ -178,15 +177,16 @@ def test_gpindices(self):
 
     def test_circuit_probabilities(self):
         datagen_model = self.model.copy()
-        datagen_params = datagen_model.to_vector()
+        datagen_params = datagen_model.to_vector().copy()
         datagen_params[-2:] = [1.1,1.1]
         datagen_model.from_vector(datagen_params)
         probs = datagen_model.probabilities( (('Gxpi2',0),))
         self.assertAlmostEqual(probs['0'],0.8247240241650917)
 
     def test_germ_selection(self):
         datagen_model = self.model.copy()
-        datagen_params = datagen_model.to_vector()
+
+        datagen_params = datagen_model.to_vector().copy()
         datagen_params[-2:] = [1.1,1.1]
         datagen_model.from_vector(datagen_params)
 
@@ -198,7 +198,15 @@ def test_germ_selection(self):
         self.assertEqual(final_germs, [pygsti.circuits.circuit.Circuit('Gxpi2:0')])
 
 
+    def test_modifying_view_of_model_params_desynchronizes_the_global_parameter_vector(self):
+        datagen_model = self.model.copy()
+
+        datagen_params = datagen_model.to_vector()
+        datagen_params[-2:] = [1.1,1.1]
 
+        with self.assertRaises(ValueError):
+            datagen_model._check_paramvec()
+        # Modify a copy of the parameter array instead!
 
 
 

test/unit/objects/test_composed_spam.py

@@ -47,12 +47,12 @@ def test_forward_simulation(self):
         indep_mdl['rho0'] = pure_vec
         indep_mdl['G0'] = noise_op
         indep_mdl['Mdefault'] = self.base_povm
-        indep_mdl.num_params  # triggers paramvec rebuild
-        
+        indep_mdl._rebuild_paramvec()  # trigger paramvec rebuild directly.
+
         composed_mdl = ExplicitOpModel(['Q0'], evotype='default')
         composed_mdl['rho0'] = self.vec
         composed_mdl['Mdefault'] = self.base_povm
-        composed_mdl.num_params  # triggers paramvec rebuild
+        composed_mdl._rebuild_paramvec()  # trigger paramvec rebuild directly.
 
         # Sanity check
         indep_circ = Circuit(['rho0', 'G0', 'Mdefault'])
@@ -167,7 +167,7 @@ def test_forward_simulation(self):
         indep_mdl['rho0'] = self.base_prep
         indep_mdl['G0'] = noise_op.copy()
         indep_mdl['Mdefault'] = pure_povm
-        indep_mdl._clean_paramvec()
+        indep_mdl._rebuild_paramvec()
 
         composed_mdl = ExplicitOpModel(['Q0'], evotype='default')
         composed_mdl['rho0'] = self.base_prep

test/unit/objects/test_model_updates_after_creation.py

@@ -0,0 +1,181 @@
+
+from typing import Literal
+import numpy as np
+from pygsti.baseobjs import qubitgraph as _qgraph
+from pygsti.baseobjs import QubitSpace
+from pygsti.models import modelconstruction as pgmc
+from pygsti.processors import QubitProcessorSpec
+from pygsti.modelmembers.states import ComposedState, ComputationalBasisState
+from pygsti.modelmembers.povms import ComposedPOVM
+from pygsti.modelmembers.operations import LindbladErrorgen, ExpErrorgenOp
+from pygsti.baseobjs.errorgenbasis import CompleteElementaryErrorgenBasis
+from pygsti.tools import slicetools as _slct
+from pygsti.modelmembers.operations import ComposedOp, EmbeddedOp
+from pygsti.algorithms import BuiltinBasis
+from pygsti.modelmembers.operations import create_from_unitary_mx
+
+from ..util import BaseCase
+
+#region Create Model
+
+
+def make_spam(num_qubits):
+    state_space = QubitSpace(num_qubits)
+    max_weights = {'H': 1, 'S': 1, 'C': 1, 'A': 1}
+    egbn_hamiltonian_only = CompleteElementaryErrorgenBasis(BuiltinBasis("PP", 4),
+                                                            state_space,
+                                                            ('H', ),
+                                                            max_weights)
+
+    rho_errgen_rates = {ell: 0.0 for ell in egbn_hamiltonian_only.labels}
+    rho_lindblad = LindbladErrorgen.from_elementary_errorgens(rho_errgen_rates,
+                                                              parameterization='H',
+                                                              state_space=state_space,
+                                                              evotype='densitymx')
+    rho_errorgen = ExpErrorgenOp(rho_lindblad)
+    rho_ideal = ComputationalBasisState([0] * num_qubits)
+    rho = ComposedState(rho_ideal, rho_errorgen)
+
+    povm_errgen_rates = {ell: 0.0 for ell in egbn_hamiltonian_only.labels}
+    povm_linblad = LindbladErrorgen.from_elementary_errorgens(povm_errgen_rates,
+                                                              parameterization='H',
+                                                              state_space=state_space,
+                                                              evotype='densitymx')
+
+    measure = ComposedPOVM(ExpErrorgenOp(povm_linblad))
+
+    return rho, measure
+
+
+def make_target_model(num_qubits, independent_gates: bool = True):
+    ps_geometry = _qgraph.QubitGraph.common_graph(
+        num_qubits, geometry='line',
+        directed=True, all_directions=True,
+        qubit_labels=tuple(range(num_qubits))
+    )
+    u_ecr = 1 / np.sqrt(2) * np.array([[0, 0, 1, 1j],
+                                       [0, 0, 1j, 1],
+                                       [1, -1j, 0, 0],
+                                       [-1j, 1, 0, 0]])
+
+    gatenames = ["Gxpi2", "Gi", "Gecr"]
+    ps = QubitProcessorSpec(
+        num_qubits=num_qubits,
+        gate_names=gatenames,
+        nonstd_gate_unitaries={'Gecr': u_ecr},
+        geometry=ps_geometry
+    )
+    gateerrs = {}
+    basis = BuiltinBasis("PP", QubitSpace(1))
+    egb1 = CompleteElementaryErrorgenBasis(basis, QubitSpace(1), ('H', 'S'))
+    for gn in gatenames[:-1]:
+        gateerrs[gn] = {ell: 0 for ell in egb1.labels}
+    egb2 = CompleteElementaryErrorgenBasis(basis, QubitSpace(2), ('H', 'S'))
+    gateerrs['Gecr'] = {ell: 0 for ell in egb2.labels}
+
+    tmn = pgmc.create_crosstalk_free_model(ps, lindblad_error_coeffs=gateerrs, independent_gates=independent_gates)
+
+    return tmn
+#endregion Create Model
+
+
+class TestRebuildParamVec(BaseCase):
+
+    def __init__(self, methodName="runTest"):
+        super().__init__(methodName)
+        self.num_qubits = 2
+
+    def setup_model_for_test(self, where: Literal["op", "sp", "povm"]):
+
+        model = make_target_model(self.num_qubits)
+
+        paramlbls = model._paramlbls.copy()
+        obj_of_interest = None
+        if where in ["povm", "sp"]:
+            rho, povm = make_spam(self.num_qubits)
+            if where == "sp":
+                model.prep_blks['layers']['rho0'] = rho
+                obj_of_interest = rho
+            else:
+                model.povm_blks['layers']['Mdefault'] = povm
+                obj_of_interest = povm
+
+        elif where == "op":
+            gate_name = "Gypi2"
+            egb1 = CompleteElementaryErrorgenBasis(BuiltinBasis("PP", QubitSpace(1)),
+                                                   QubitSpace(1), ('H', 'S'))
+
+            new_gate_errgen_rates = {ell: 0.0 for ell in egb1.labels}
+            my_lindbladian = LindbladErrorgen.from_elementary_errorgens(new_gate_errgen_rates,
+                                                                        parameterization='auto',
+                                                                        state_space=QubitSpace(1),
+                                                                        evotype='densitymx')
+
+            comp = ComposedOp([create_from_unitary_mx(np.eye(2), "static standard",
+                                                      stdname=gate_name),
+                               ExpErrorgenOp(my_lindbladian)])
+            comp1 = comp.copy()
+
+            model.operation_blks["gates"][(gate_name, 0)] = comp
+            model.operation_blks["gates"][(gate_name, 1)] = comp1
+
+            # Add the embedded op in a layer as well.
+
+            embedded = EmbeddedOp(QubitSpace(self.num_qubits), [0], comp)
+            embedded1 = EmbeddedOp(QubitSpace(self.num_qubits), [1], comp1)
+            model.operation_blks["layers"][(gate_name, 0)] = embedded
+            model.operation_blks["layers"][(gate_name, 1)] = embedded1
+
+            obj_of_interest = (comp, comp1)
+        else:
+            raise ValueError("Unexpected value for the type of new layer to add to the model.")
+
+        model._rebuild_paramvec()
+
+        return paramlbls, model, model._paramlbls.copy(), obj_of_interest
+
+    def test_add_state_prep_after_creation_of_implicit_noisy_model(self):
+
+        original_lbls, model, modified_lbls, rho = self.setup_model_for_test("sp")
+
+        inds = _slct.indices_as_array(rho._gpindices)
+        avail_inds = np.arange(model.num_params)
+        cross_check = np.where(avail_inds[:, None] != inds[None, :], 1, 0)
+
+        totals = np.sum(cross_check, axis=1)
+        self.assertEqual(len(totals), model.num_params)
+
+        used_inds = np.where(totals == len(inds))
+        self.assertArraysEqual(original_lbls, modified_lbls[used_inds])
+
+    def test_add_povm_after_creation_of_implicit_noisy_model(self):
+
+        original_lbls, model, modified_lbls, povm = self.setup_model_for_test("povm")
+
+        inds = _slct.indices_as_array(povm._gpindices)
+        avail_inds = np.arange(model.num_params)
+        cross_check = np.where(avail_inds[:, None] != inds[None, :], 1, 0)
+
+        totals = np.sum(cross_check, axis=1)
+        self.assertEqual(len(totals), model.num_params)
+
+        used_inds = np.where(totals == len(inds))
+        self.assertArraysEqual(original_lbls, modified_lbls[used_inds])
+
+    def test_add_gate_operation_after_creation_of_implicit_noisy_model(self):
+
+        original_lbls, model, modified_lbls, (comp, comp1) = self.setup_model_for_test("op")
+
+        excluded_inds = np.array(list(_slct.to_array(comp.gpindices)) + list(_slct.to_array(comp1.gpindices)))
+        avail_inds = np.arange(model.num_params)
+        cross_check = np.where(avail_inds[:, None] != excluded_inds[None, :], 1, 0)
+        totals = np.sum(cross_check, axis=1)
+        self.assertEqual(len(totals), model.num_params)
+
+        norm = np.linalg.norm(comp.gpindices_as_array() - comp1.gpindices_as_array(), 2)
+        tol = 2**-53  # episilon machine for double precision.
+        msg = "Composed Op copy is treated as the same object as the original in creation of gpindices."
+        self.assertGreater(norm, tol, msg)
+
+        used_inds = np.where(totals == len(excluded_inds))
+        self.assertArraysEqual(original_lbls, modified_lbls[used_inds])