Expire Dense Interfaces

pygsti/algorithms/contract.py

@@ -19,6 +19,7 @@
 from pygsti import tools as _tools
 from pygsti.modelmembers import operations as _op
 from pygsti.modelmembers import povms as _povm
+from pygsti.modelmembers import states as _state
 
 
 def contract(model, to_what, dataset=None, maxiter=1000000, tol=0.01, use_direct_cp=True, method="Nelder-Mead",
@@ -192,7 +193,7 @@ def _contract_to_cp_direct(model, verbosity, tp_also=False, maxiter=100000, tol=
     printer.log(("--- Contract to %s (direct) ---" % ("CPTP" if tp_also else "CP")), 1)
 
     for (opLabel, gate) in model.operations.items():
-        new_op = gate.copy()
+        new_op = gate.to_dense()
         if(tp_also):
             for k in range(new_op.shape[1]): new_op[0, k] = 1.0 if k == 0 else 0.0
 
@@ -294,16 +295,18 @@ def _contract_to_cp_direct(model, verbosity, tp_also=False, maxiter=100000, tol=
         #else: print "contract_to_cp_direct success in %d iterations" % it  #DEBUG
 
         printer.log("Direct CP contraction of %s gate gives frobenius diff of %g" %
-                    (opLabel, _tools.frobeniusdist(mdl.operations[opLabel], gate)), 2)
+                    (opLabel, _tools.frobeniusdist(mdl.operations[opLabel].to_dense(), gate.to_dense())), 2)
 
     #mdl.log("Choi-Truncate to %s" % ("CPTP" if tp_also else "CP"), { 'maxiter': maxiter } )
     distance = mdl.frobeniusdist(model)
     printer.log(('The closest legal point found was distance: %s' % str(distance)), 1)
 
     if tp_also:  # TP also constrains prep vectors
         op_dim = mdl.dim
-        for rhoVec in list(mdl.preps.values()):
-            rhoVec[0, 0] = 1.0 / op_dim**0.25
+        for lbl, rhoVec in mdl.preps.items():
+            new_rho_vec = rhoVec.to_dense()
+            new_rho_vec[0] = 1.0 / op_dim**0.25
+            mdl.preps[lbl] = new_rho_vec
 
     mdl._need_to_rebuild = True
     return distance, mdl
@@ -315,13 +318,18 @@ def _contract_to_tp(model, verbosity):
     #printer.log('', 2)
     printer.log("--- Contract to TP ---", 1)
     mdl = model.copy()
-    for gate in list(mdl.operations.values()):
-        gate[0, 0] = 1.0
-        for k in range(1, gate.shape[1]): gate[0, k] = 0.0
+    for lbl, gate in mdl.operations.items():
+        new_gate_mx = gate.to_dense()
+        new_gate_mx[0, 0] = 1.0
+        for k in range(1, new_gate_mx.shape[1]): 
+            new_gate_mx[0, k] = 0.0
+        mdl.operations[lbl] = new_gate_mx
 
     op_dim = mdl.dim
-    for rhoVec in list(mdl.preps.values()):
-        rhoVec[0, 0] = 1.0 / op_dim**0.25
+    for lbl, rhoVec in mdl.preps.items():
+        new_rho_vec = rhoVec.to_dense()
+        new_rho_vec[0] = 1.0 / op_dim**0.25
+        mdl.preps[lbl] = new_rho_vec
 
     mdl._need_to_rebuild = True
     distance = mdl.frobeniusdist(model)
@@ -363,20 +371,20 @@ def _contract_to_valid_spam(model, verbosity=0):
 
     # rhoVec must be positive semidefinite and trace = 1
     for prepLabel, rhoVec in mdl.preps.items():
-        vec = rhoVec.copy()
+        vec = rhoVec.to_dense()
 
         #Ensure trace == 1.0 (maybe later have this be optional)
         firstElTarget = 1.0 / firstElTrace  # TODO: make this function more robust
         # to multiple rhovecs -- can only use on ratio,
         # so maybe take average of ideal ratios for each rhoVec
         # and apply that?  The function works fine now for just one rhovec.
-        if abs(firstElTarget - vec[0, 0]) > TOL:
-            r = firstElTarget / vec[0, 0]
+        if abs(firstElTarget - vec[0]) > TOL:
+            r = firstElTarget / vec[0]
             vec *= r  # multiply rhovec by factor
             for povmLbl in list(mdl.povms.keys()):
                 scaled_effects = []
                 for ELabel, EVec in mdl.povms[povmLbl].items():
-                    scaled_effects.append((ELabel, EVec / r))
+                    scaled_effects.append((ELabel, EVec.to_dense() / r))
                 # Note: always creates an unconstrained POVM
                 mdl.povms[povmLbl] = _povm.UnconstrainedPOVM(scaled_effects, mdl.evotype, mdl.state_space)
 
@@ -386,11 +394,11 @@ def _contract_to_valid_spam(model, verbosity=0):
         lowEval = min([ev.real for ev in _np.linalg.eigvals(mx)])
         while(lowEval < -TOL):
             # only element with trace (even for multiple qubits) -- keep this constant and decrease others
-            idEl = vec[0, 0]
-            vec /= 1.00001; vec[0, 0] = idEl
+            idEl = vec[0]
+            vec /= 1.00001; vec[0] = idEl
             lowEval = min([ev.real for ev in _np.linalg.eigvals(_tools.ppvec_to_stdmx(vec))])
 
-        diff += _np.linalg.norm(model.preps[prepLabel] - vec)
+        diff += _np.linalg.norm(model.preps[prepLabel].to_dense() - vec)
         mdl.preps[prepLabel] = vec
 
     # EVec must have eigenvals between 0 and 1 <==> positive semidefinite and trace <= 1
@@ -399,6 +407,7 @@ def _contract_to_valid_spam(model, verbosity=0):
         for ELabel, EVec in mdl.povms[povmLbl].items():
             #if isinstance(EVec, _povm.ComplementPOVMEffect):
             #    continue #don't contract complement vectors
+            EVec = EVec.to_dense()
             evals, evecs = _np.linalg.eig(_tools.ppvec_to_stdmx(EVec))
             if(min(evals) < 0.0 or max(evals) > 1.0):
                 if all([ev > 1.0 for ev in evals]):
@@ -410,7 +419,7 @@ def _contract_to_valid_spam(model, verbosity=0):
                     if ev > 1.0: evals[k] = 1.0
                 mx = _np.dot(evecs, _np.dot(_np.diag(evals), _np.linalg.inv(evecs)))
                 vec = _tools.stdmx_to_ppvec(mx)
-                diff += _np.linalg.norm(model.povms[povmLbl][ELabel] - vec)
+                diff += _np.linalg.norm(model.povms[povmLbl][ELabel].to_dense() - vec)
                 scaled_effects.append((ELabel, vec))
             else:
                 scaled_effects.append((ELabel, EVec))  # no scaling
@@ -423,12 +432,12 @@ def _contract_to_valid_spam(model, verbosity=0):
     printer.log("--- Contract to valid SPAM ---", 1)
     printer.log(("Sum of norm(deltaE) and norm(deltaRho) = %g" % diff), 1)
     for (prepLabel, rhoVec) in model.preps.items():
-        printer.log("  %s: %s ==> %s " % (prepLabel, str(_np.transpose(rhoVec)),
-                                          str(_np.transpose(mdl.preps[prepLabel]))), 2)
+        printer.log("  %s: %s ==> %s " % (prepLabel, str(_np.transpose(rhoVec.to_dense())),
+                                          str(_np.transpose(mdl.preps[prepLabel].to_dense()))), 2)
     for povmLbl, povm in model.povms.items():
         printer.log(("  %s (POVM)" % povmLbl), 2)
         for ELabel, EVec in povm.items():
-            printer.log("  %s: %s ==> %s " % (ELabel, str(_np.transpose(EVec)),
-                                              str(_np.transpose(mdl.povms[povmLbl][ELabel]))), 2)
+            printer.log("  %s: %s ==> %s " % (ELabel, str(_np.transpose(EVec.to_dense())),
+                                              str(_np.transpose(mdl.povms[povmLbl][ELabel].to_dense()))), 2)
 
     return mdl  # return contracted model

pygsti/algorithms/gaugeopt.py

@@ -637,13 +637,13 @@ def _objective_fn(gauge_group_el, oob_check):
                     for opLbl in mdl.operations:
                         wt = item_weights.get(opLbl, opWeight)
                         ret += wt * (1.0 - _tools.entanglement_fidelity(
-                            target_model.operations[opLbl], mdl.operations[opLbl]))**2
+                            target_model.operations[opLbl].to_dense(), mdl.operations[opLbl].to_dense()))**2
 
                 elif gates_metric == "tracedist":
                     for opLbl in mdl.operations:
                         wt = item_weights.get(opLbl, opWeight)
                         ret += opWeight * _tools.jtracedist(
-                            target_model.operations[opLbl], mdl.operations[opLbl])
+                            target_model.operations[opLbl].to_dense(), mdl.operations[opLbl].to_dense())
 
                 else: raise ValueError("Invalid gates_metric: %s" % gates_metric)
 
@@ -666,9 +666,9 @@ def _objective_fn(gauge_group_el, oob_check):
                 elif spam_metric == "fidelity":
                     for preplabel, prep in mdl.preps.items():
                         wt = item_weights.get(preplabel, spamWeight)
-                        rhoMx1 = _tools.vec_to_stdmx(prep, mxBasis)
+                        rhoMx1 = _tools.vec_to_stdmx(prep.to_dense(), mxBasis)
                         rhoMx2 = _tools.vec_to_stdmx(
-                            target_model.preps[preplabel], mxBasis)
+                            target_model.preps[preplabel].to_dense(), mxBasis)
                         ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2
 
                     for povmlabel, povm in mdl.povms.items():
@@ -679,9 +679,9 @@ def _objective_fn(gauge_group_el, oob_check):
                 elif spam_metric == "tracedist":
                     for preplabel, prep in mdl.preps.items():
                         wt = item_weights.get(preplabel, spamWeight)
-                        rhoMx1 = _tools.vec_to_stdmx(prep, mxBasis)
+                        rhoMx1 = _tools.vec_to_stdmx(prep.to_dense(), mxBasis)
                         rhoMx2 = _tools.vec_to_stdmx(
-                            target_model.preps[preplabel], mxBasis)
+                            target_model.preps[preplabel].to_dense(), mxBasis)
                         ret += wt * _tools.tracedist(rhoMx1, rhoMx2)
 
                     for povmlabel, povm in mdl.povms.items():
@@ -768,11 +768,11 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl_pre, mdl_post,
     dS = _np.rollaxis(dS, 2)  # shape (n, d1, d2)
 
     for i, (gl, gate) in enumerate(mdl_post.operations.items()):
-        pre_op = mdl_pre.operations[gl]
+        pre_op = mdl_pre.operations[gl].to_dense()
 
         #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)
+        chi = _tools.fast_jamiolkowski_iso_std(gate.to_dense(), op_basis)
         assert(_np.linalg.norm(chi - chi.T.conjugate()) < 1e-4), \
             "chi should be Hermitian!"
 
@@ -785,7 +785,7 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl_pre, mdl_post,
         # transform).  This shuffle op commutes with the derivative, so that
         # dchi_std/dp := d(M(G'))/dp = M(d(S_inv*G*S)/dp) = M( d(S_inv)*G*S + S_inv*G*dS )
         #              = M( (-S_inv*dS*S_inv)*G*S + S_inv*G*dS ) = M( S_inv*(-dS*S_inv*G*S) + G*dS )
-        left = -1 * _np.dot(dS, gate)  # shape (n,d1,d2)
+        left = -1 * _np.dot(dS, gate.to_dense())  # shape (n,d1,d2)
         right = _np.swapaxes(_np.dot(pre_op, dS), 0, 1)  # shape (d1, n, d2) -> (n,d1,d2)
         result = _np.swapaxes(_np.dot(S_inv, left + right), 0, 1)  # shape (n, d1, d2)
 

pygsti/extras/rpe/rpeconstruction.py

@@ -95,9 +95,9 @@ def create_parameterized_rpe_model(alpha_true, epsilon_true, aux_rot, spam_depol
             effect_labels=ELabels, effect_expressions=EExpressions)
 
         outputModel.operations[loose_axis_gate_label] = \
-            _np.dot(_np.dot(_np.linalg.inv(modelAux1.operations[auxiliary_axis_gate_label]),
-                            outputModel.operations[loose_axis_gate_label]),
-                    modelAux1.operations[auxiliary_axis_gate_label])
+            _np.linalg.inv(modelAux1.operations[auxiliary_axis_gate_label].to_dense())@\
+            outputModel.operations[loose_axis_gate_label].to_dense()@\
+            modelAux1.operations[auxiliary_axis_gate_label].to_dense()
 
     outputModel = outputModel.depolarize(op_noise=gate_depol,
                                          spam_noise=spam_depol)

pygsti/extras/rpe/rpetools.py

@@ -257,7 +257,7 @@ def extract_alpha(model, rpeconfig_inst):
         The value of alpha for the input model.
     """
     op_label = rpeconfig_inst.fixed_axis_gate_label
-    decomp = _decompose_gate_matrix(model.operations[op_label])
+    decomp = _decompose_gate_matrix(model.operations[op_label].to_dense())
     alphaVal = decomp['pi rotations'] * _np.pi
     return alphaVal
 
@@ -284,7 +284,7 @@ def extract_epsilon(model, rpeconfig_inst):
         The value of epsilon for the input model.
     """
     op_label = rpeconfig_inst.loose_axis_gate_label
-    decomp = _decompose_gate_matrix(model.operations[op_label])
+    decomp = _decompose_gate_matrix(model.operations[op_label].to_dense())
 
     epsilonVal = decomp['pi rotations'] * _np.pi
     return epsilonVal
@@ -313,10 +313,10 @@ def extract_theta(model, rpeconfig_inst):
         The value of theta for the input model.
     """
     op_label = rpeconfig_inst.loose_axis_gate_label
-    decomp = _decompose_gate_matrix(model.operations[op_label])
+    decomp = _decompose_gate_matrix(model.operations[op_label].to_dense())
     target_axis = rpeconfig_inst.loose_axis_target
 
-    decomp = _decompose_gate_matrix(model.operations[op_label])
+    decomp = _decompose_gate_matrix(model.operations[op_label].to_dense())
     thetaVal = _np.real_if_close([_np.arccos(
         _np.dot(decomp['axis of rotation'], target_axis))])[0]
     if thetaVal > _np.pi / 2:

pygsti/modelmembers/instruments/instrument.py

@@ -226,7 +226,7 @@ def num_elements(self):
         -------
         int
         """
-        return sum([g.size for g in self.values()])
+        return sum([g._ptr.size for g in self.values()])
 
     @property
     def num_params(self):

pygsti/modelmembers/instruments/tpinstrument.py

@@ -86,11 +86,12 @@ def __init__(self, op_matrices, evotype="default", state_space=None, called_from
 
         if op_matrices is not None:
             if isinstance(op_matrices, dict):
-                matrix_list = [(k, v) for k, v in op_matrices.items()]  # gives definite ordering
+                matrix_list = [(k, v[()]) if (isinstance(v,_np.ndarray) and v.shape == ()) else (k,v) for k, v in op_matrices.items()]  # gives definite ordering
             elif isinstance(op_matrices, list):
-                matrix_list = op_matrices  # assume it's is already an ordered (key,value) list
+                matrix_list = [(k, v[()]) if (isinstance(v,_np.ndarray) and v.shape == ()) else (k,v) for k, v in op_matrices]
             else:
                 raise ValueError("Invalid `op_matrices` arg of type %s" % type(op_matrices))
+            matrix_list = [(k,v.to_dense()) if isinstance(v, _mm.ModelMember) else (k,v) for k,v in matrix_list]
 
             assert(len(matrix_list) > 0 or state_space is not None), \
                 "Must specify `state_space` when there are no instrument members!"
@@ -99,7 +100,7 @@ def __init__(self, op_matrices, evotype="default", state_space=None, called_from
             evotype = _Evotype.cast(evotype, state_space=state_space)
 
             # Create gate objects that are used to parameterize this instrument
-            MT_mx = sum([v for k, v in matrix_list])  # sum-of-instrument-members matrix
+            MT_mx = sum([v for _, v in matrix_list])  # sum-of-instrument-members matrix
             MT = _op.FullTPOp(MT_mx, None, evotype, state_space)
             self.param_ops.append(MT)
 
@@ -114,14 +115,6 @@ def __init__(self, op_matrices, evotype="default", state_space=None, called_from
             items = [(k, _TPInstrumentOp(self.param_ops, i, None))
                      for i, (k, v) in enumerate(matrix_list)]
 
-            #DEBUG
-            #print("POST INIT PARAM GATES:")
-            #for i,v in enumerate(self.param_ops):
-            #    print(i,":\n",v)
-            #
-            #print("POST nINIT ITEMS:")
-            #for k,v in items:
-            #    print(k,":\n",v)
         else:
             assert(state_space is not None), "`state_space` cannot be `None` when there are no members!"
             evotype = _Evotype.cast(evotype, state_space=state_space)
@@ -273,7 +266,7 @@ def num_elements(self):  # same as in Instrument CONSOLIDATE?
         -------
         int
         """
-        return sum([g.size for g in self.values()])
+        return sum([g._ptr.size for g in self.values()])
 
     @property
     def num_params(self):  # same as in Instrument CONSOLIDATE?

pygsti/modelmembers/instruments/tpinstrumentop.py

@@ -122,13 +122,13 @@ def _construct_matrix(self):
         nEls = self.num_instrument_elements
         self._ptr.flags.writeable = True
         if self.index < nEls - 1:
-            self._ptr[:, :] = _np.asarray(self.relevant_param_ops[1]  # i.e. param_ops[self.index + 1]
-                                          + self.relevant_param_ops[0])  # i.e. param_ops[0]
+            self._ptr[:, :] = _np.asarray(self.relevant_param_ops[1].to_dense()  # i.e. param_ops[self.index + 1]
+                                          + self.relevant_param_ops[0].to_dense())  # i.e. param_ops[0]
         else:
             assert(self.index == nEls - 1), \
                 "Invalid index %d > %d" % (self.index, nEls - 1)
-            self._ptr[:, :] = _np.asarray(-sum(self.relevant_param_ops)  # all instrument param_ops == relevant
-                                          - (nEls - 3) * self.relevant_param_ops[0])
+            self._ptr[:, :] = _np.asarray(-sum([op.to_dense() for op in self.relevant_param_ops])  # all instrument param_ops == relevant
+                                          - (nEls - 3) * self.relevant_param_ops[0].to_dense())
 
         assert(self._ptr.shape == (self.dim, self.dim))
         self._ptr.flags.writeable = False

pygsti/modelmembers/operations/__init__.py

@@ -15,7 +15,7 @@
 
 from .composederrorgen import ComposedErrorgen
 from .composedop import ComposedOp
-from .denseop import DenseOperator, DenseOperatorInterface
+from .denseop import DenseOperator
 from .depolarizeop import DepolarizeOp
 from .eigpdenseop import EigenvalueParamDenseOp
 from .embeddederrorgen import EmbeddedErrorgen