fix bug in previous change to eigendenseop. Finish going through the …

…codebase.

pygsti/modelmembers/operations/eigpdenseop.py

@@ -432,7 +432,7 @@ def deriv_wrt_params(self, wrt_filter=None):
             dMx = _np.zeros((self.dim, self.dim), 'complex')
             for prefactor, (i, j) in pdesc:
                 dMx[i, j] = prefactor
-            tmp = self.B @ (dMx, self.Bi)
+            tmp = self.B @ (dMx @ self.Bi)
             if _np.linalg.norm(tmp.imag) >= IMAG_TOL:  # just a warning until we figure this out.
                 print("EigenvalueParamDenseOp deriv_wrt_params WARNING:"
                       " Imag part = ", _np.linalg.norm(tmp.imag), " pdesc = ", pdesc)  # pragma: no cover

pygsti/modelmembers/operations/linearop.py

@@ -483,12 +483,14 @@ def residuals(self, other_op, transform=None, inv_transform=None):
         numpy.ndarray
             A 1D-array of size equal to that of the flattened operation matrix.
         """
-        if transform is None and inv_transform is None:
-            return _ot.residuals(self.to_dense(on_space='minimal'), other_op.to_dense(on_space='minimal'))
+        dense_self = self.to_dense(on_space='minimal')
+        if transform is not None:
+            assert inv_transform is not None
+            dense_self = inv_transform @ (dense_self @ transform)
         else:
-            return _ot.residuals(_np.dot(
-                inv_transform, _np.dot(self.to_dense(on_space='minimal'), transform)),
-                other_op.to_dense(on_space='minimal'))
+            assert inv_transform is None
+        return (dense_self - other_op.to_dense(on_space='minimal')).ravel()
+
 
     def jtracedist(self, other_op, transform=None, inv_transform=None):
         """

pygsti/modelmembers/povms/effect.py

@@ -170,11 +170,9 @@ def residuals(self, other_spam_vec, transform=None, inv_transform=None):
         float
         """
         vec = self.to_dense()
-        if transform is None:
-            return _ot.residuals(vec, other_spam_vec.to_dense())
-        else:
-            return _ot.residuals(_np.dot(_np.transpose(transform),
-                                         vec), other_spam_vec.to_dense())
+        if transform is not None:
+            vec = transform.T @ vec
+        return (vec - other_spam_vec.to_dense()).ravel()
 
     def transform_inplace(self, s):
         """

pygsti/modelmembers/states/state.py

@@ -350,11 +350,10 @@ def residuals(self, other_spam_vec, transform=None, inv_transform=None):
         float
         """
         vec = self.to_dense(on_space='minimal')
-        if inv_transform is None:
-            return _ot.residuals(vec, other_spam_vec.to_dense(on_space='minimal'))
-        else:
-            return _ot.residuals(_np.dot(inv_transform, vec),
-                                 other_spam_vec.to_dense(on_space='minimal'))
+        if inv_transform is not None:
+            vec = inv_transform @ vec
+        return (vec - other_spam_vec.to_dense(on_space='minimal')).ravel()
+
 
     def transform_inplace(self, s):
         """

pygsti/report/workspaceplots.py

@@ -2588,7 +2588,7 @@ def _create(self, evals_list, colors, labels, scale, amp, center_text):
                 amp_evals = evals**amp
                 trace = go.Scatterpolar(
                     r=list(_np.absolute(amp_evals).flat),
-                    theta=list(_np.angle(amp_evals).flatten() * (180.0 / _np.pi)),
+                    theta=list(_np.angle(amp_evals).ravel() * (180.0 / _np.pi)),
                     showlegend=False,
                     mode='markers',
                     marker=dict(
@@ -2899,15 +2899,16 @@ def __init__(self, ws, evals, errbars=None, scale=1.0):
 
     def _create(self, evals, errbars, scale):
 
+        flat_errbars = errbars.ravel()
         HOVER_PREC = 7
         xs = list(range(evals.size))
-        ys = []; colors = []; texts = []
-        for i, ev in enumerate(evals.flatten()):
+        ys, colors, texts = [], [], []
+        for i, ev in enumerate(evals.ravel()):
             ys.append(abs(ev.real))
             colors.append('rgb(200,200,200)' if ev.real > 0 else 'red')
             if errbars is not None:
                 texts.append("%g +/- %g" % (round(ev.real, HOVER_PREC),
-                                            round(errbars.flatten()[i].real, HOVER_PREC)))
+                                            round(flat_errbars[i].real, HOVER_PREC)))
             else:
                 texts.append("%g" % round(ev.real, HOVER_PREC))
 
@@ -3033,13 +3034,13 @@ def _create(self, dataset, target, maxlen, fixed_lists, scale):
 
         xs = list(range(svals.size))
         trace1 = go.Bar(
-            x=xs, y=list(svals.flatten()),
+            x=xs, y=list(svals.flat),
             marker=dict(color="blue"),
             hoverinfo='y',
             name="from Data"
         )
         trace2 = go.Bar(
-            x=xs, y=list(target_svals.flatten()),
+            x=xs, y=list(target_svals.flat),
             marker=dict(color="black"),
             hoverinfo='y',
             name="from Target"

pygsti/tools/jamiolkowski.py

@@ -225,7 +225,7 @@ def fast_jamiolkowski_iso_std(operation_mx, op_mx_basis):
     N2 = opMxInStdBasis.shape[0]; N = int(_np.sqrt(N2))
     assert(N * N == N2)  # make sure N2 is a perfect square
     Jmx = opMxInStdBasis.reshape((N, N, N, N))
-    Jmx = _np.swapaxes(Jmx, 1, 2).flatten()
+    Jmx = _np.swapaxes(Jmx, 1, 2).ravel()
     Jmx = Jmx.reshape((N2, N2))
 
     # This construction results in a Jmx with trace == dim(H) = sqrt(gateMxInPauliBasis.shape[0])
@@ -261,7 +261,7 @@ def fast_jamiolkowski_iso_std_inv(choi_mx, op_mx_basis):
     N2 = choi_mx.shape[0]; N = int(_np.sqrt(N2))
     assert(N * N == N2)  # make sure N2 is a perfect square
     opMxInStdBasis = choi_mx.reshape((N, N, N, N)) * N
-    opMxInStdBasis = _np.swapaxes(opMxInStdBasis, 1, 2).flatten()
+    opMxInStdBasis = _np.swapaxes(opMxInStdBasis, 1, 2).ravel()
     opMxInStdBasis = opMxInStdBasis.reshape((N2, N2))
     op_mx_basis = _bt.create_basis_for_matrix(opMxInStdBasis, op_mx_basis)
 

pygsti/tools/lindbladtools.py

@@ -83,18 +83,21 @@ def create_elementary_errorgen_dual(typ, p, q=None, sparse=False, normalization_
             rho1 = (p @ rho0 @ qdag + q @ rho0 @ pdag)  # 1 / (2 * d2) *
         elif typ == 'A':
             rho1 = 1j * (p @ rho0 @ qdag - q @ rho0 @ pdag)  # 1j / (2 * d2)
-        elem_errgen[:, i] = rho1.flatten()[:, None] if sparse else rho1.flatten()
+        elem_errgen[:, i] = rho1.ravel()
+        # ^ That line used to branch depending on the value of "sparse", but it
+        #   turns out that both codepaths produced the same result.
 
     return_normalization = bool(normalization_factor == 'auto_return')
     if normalization_factor in ('auto', 'auto_return'):
         primal = create_elementary_errorgen(typ, p, q, sparse)
         if sparse:
-            normalization_factor = _np.vdot(elem_errgen.toarray().flatten(), primal.toarray().flatten())
+            normalization_factor = _np.vdot(elem_errgen.toarray(), primal.toarray())
         else:
-            normalization_factor = _np.vdot(elem_errgen.flatten(), primal.flatten())
+            normalization_factor = _np.vdot(elem_errgen, primal)
     elem_errgen *= _np.real_if_close(1 / normalization_factor).item()  # item() -> scalar
 
-    if sparse: elem_errgen = elem_errgen.tocsr()
+    if sparse:
+        elem_errgen = elem_errgen.tocsr()
     return (elem_errgen, normalization_factor) if return_normalization else elem_errgen
 
 
@@ -133,7 +136,8 @@ def create_elementary_errorgen(typ, p, q=None, sparse=False):
     -------
     ndarray or Scipy CSR matrix
     """
-    d = p.shape[0]; d2 = d**2
+    d = p.shape[0]
+    d2 = d**2
     if sparse:
         elem_errgen = _sps.lil_matrix((d2, d2), dtype=p.dtype)
     else:
@@ -162,10 +166,12 @@ def create_elementary_errorgen(typ, p, q=None, sparse=False):
             rho1 = p @ rho0 @ qdag + q @ rho0 @ pdag - 0.5 * (pq_plus_qp @ rho0 + rho0 @ pq_plus_qp)
         elif typ == 'A':
             rho1 = 1j * (p @ rho0 @ qdag - q @ rho0 @ pdag + 0.5 * (pq_minus_qp @ rho0 + rho0 @ pq_minus_qp))
+        elem_errgen[:, i] = rho1.ravel()
+        # ^ That line used to branch depending on the value of sparse, but both
+        #   branches had the same effect.
 
-        elem_errgen[:, i] = rho1.flatten()[:, None] if sparse else rho1.flatten()
-
-    if sparse: elem_errgen = elem_errgen.tocsr()
+    if sparse:
+        elem_errgen = elem_errgen.tocsr()
     return elem_errgen
 
 def create_lindbladian_term_errorgen(typ, Lm, Ln=None, sparse=False):  # noqa N803
@@ -225,7 +231,10 @@ def create_lindbladian_term_errorgen(typ, Lm, Ln=None, sparse=False):  # noqa N8
         elif typ == 'O':
             rho1 = Ln @ rho0 @ Lm_dag - 0.5 * (Lmdag_Ln @ rho0 + rho0 @ Lmdag_Ln)
         else: raise ValueError("Invalid lindblad term errogen type!")
-        lind_errgen[:, i] = rho1.flatten()[:, None] if sparse else rho1.flatten()
+        lind_errgen[:, i] = rho1.ravel()
+        # ^ That line used to branch based on the value of sparse, but both branches
+        #   produced the same result.
 
-    if sparse: lind_errgen = lind_errgen.tocsr()
+    if sparse:
+        lind_errgen = lind_errgen.tocsr()
     return lind_errgen

pygsti/tools/matrixtools.py

@@ -682,11 +682,16 @@ def approximate_matrix_log(m, target_logm, target_weight=10.0, tol=1e-6):
     assert(_np.linalg.norm(m.imag) < 1e-8), "Argument `m` must be a *real* matrix!"
     mx_shape = m.shape
 
+    #
+    #   Riley note: I'd like to remove all commented-out code in this function.
+    #   @Corey or @Stefan -- you okay with that?
+    #
+
     def _objective(flat_logm):
         logM = flat_logm.reshape(mx_shape)
         testM = _spl.expm(logM)
         ret = target_weight * _np.linalg.norm(logM - target_logm)**2 + \
-            _np.linalg.norm(testM.flatten() - m.flatten(), 1)
+            _np.linalg.norm(testM.ravel() - m.ravel(), 1)
         #print("DEBUG: ",ret)
         return ret
 
@@ -703,7 +708,7 @@ def _objective(flat_logm):
     print_obj_func = None
 
     logM = _np.real(real_matrix_log(m, action_if_imaginary="ignore"))  # just drop any imaginary part
-    initial_flat_logM = logM.flatten()  # + 0.1*target_logm.flatten()
+    initial_flat_logM = logM.ravel()  # + 0.1*target_logm.flatten()
     # Note: adding some of target_logm doesn't seem to help; and hurts in easy cases
 
     if _objective(initial_flat_logM) > 1e-16:  # otherwise initial logM is fine!
@@ -1274,9 +1279,9 @@ def _fas(a, inds, rhs, add=False):
                 indx_tups = list(_itertools.product(*b))
                 inds = tuple(zip(*indx_tups))  # un-zips to one list per dim
                 if add:
-                    a[inds] += rhs.flatten()
+                    a[inds] += rhs.ravel()
                 else:
-                    a[inds] = rhs.flatten()
+                    a[inds] = rhs.ravel()
 
             #OLD DEBUG: just a reference for building the C-implementation (this is very slow in python!)
             ##Alt: C-able impl

pygsti/tools/optools.py

@@ -34,12 +34,14 @@
 def _flat_mut_blks(i, j, block_dims):
     # like _mut(i,j,dim).flatten() but works with basis *blocks*
     N = sum(block_dims)
-    mx = _np.zeros((N, N), 'd'); mx[i, j] = 1.0
+    mx = _np.zeros((N, N), 'd')
+    mx[i, j] = 1.0
     ret = _np.zeros(sum([d**2 for d in block_dims]), 'd')
     i = 0; off = 0
     for d in block_dims:
-        ret[i:i + d**2] = mx[off:off + d, off:off + d].flatten()
-        i += d**2; off += d
+        ret[i:i + d**2] = mx[off:off + d, off:off + d].ravel()
+        i += d**2
+        off += d
     return ret
 
 
@@ -221,26 +223,6 @@ def frobeniusdist_squared(a, b):
     return frobeniusdist(a, b)**2
 
 
-def residuals(a, b):
-    """
-    Calculate residuals between the elements of two matrices
-
-    Parameters
-    ----------
-    a : numpy array
-        First matrix.
-
-    b : numpy array
-        Second matrix.
-
-    Returns
-    -------
-    np.array
-        residuals
-    """
-    return (a - b).flatten()
-
-
 def tracenorm(a):
     """
     Compute the trace norm of matrix `a` given by:
@@ -1181,8 +1163,8 @@ def state_to_dmvec(psi):
         The vectorized density matrix.
     """
     psi = psi.reshape((psi.size, 1))  # convert to (N,1) shape if necessary
-    dm = _np.dot(psi, _np.conjugate(psi.T))
-    return dm.flatten()
+    dm = psi @ psi.conj().T
+    return dm.ravel()
 
 
 def dmvec_to_state(dmvec, tol=1e-6):
@@ -1661,7 +1643,7 @@ def extract_elementary_errorgen_coefficients(errorgen, elementary_errorgen_label
                                          errorgen_basis.create_simple_equivalent('std'))
     else:
         errorgen_std = _bt.change_basis(errorgen, errorgen_basis, "std")
-    flat_errorgen_std = errorgen_std.toarray().flatten() if _sps.issparse(errorgen_std) else errorgen_std.flatten()
+    flat_errorgen_std = errorgen_std.toarray().ravel() if _sps.issparse(errorgen_std) else errorgen_std.ravel()
 
     d2 = errorgen_std.shape[0]
     d = int(_np.sqrt(d2))
@@ -1677,7 +1659,7 @@ def extract_elementary_errorgen_coefficients(errorgen, elementary_errorgen_label
         bel_lbls = key.basis_element_labels
         bmx0 = elementary_errorgen_basis[bel_lbls[0]]
         bmx1 = elementary_errorgen_basis[bel_lbls[1]] if (len(bel_lbls) > 1) else None
-        flat_projector = _lt.create_elementary_errorgen_dual(key.errorgen_type, bmx0, bmx1, sparse=False).flatten()
+        flat_projector = _lt.create_elementary_errorgen_dual(key.errorgen_type, bmx0, bmx1, sparse=False).ravel()
         projections[key] = _np.real_if_close(_np.vdot(flat_projector, flat_errorgen_std), tol=1000)
         if return_projected_errorgen:
             space_projector[:, i] = flat_projector
@@ -1757,7 +1739,7 @@ def project_errorgen(errorgen, elementary_errorgen_type, elementary_errorgen_bas
                                          errorgen_basis.create_simple_equivalent('std'))
     else:
         errorgen_std = _bt.change_basis(errorgen, errorgen_basis, "std")
-    flat_errorgen_std = errorgen_std.toarray().flatten() if _sps.issparse(errorgen_std) else errorgen_std.flatten()
+    flat_errorgen_std = errorgen_std.toarray().ravel() if _sps.issparse(errorgen_std) else errorgen_std.ravel()
 
     d2 = errorgen_std.shape[0]
     d = int(_np.sqrt(d2))
@@ -1771,7 +1753,7 @@ def project_errorgen(errorgen, elementary_errorgen_type, elementary_errorgen_bas
         space_projector = _np.empty((d2 * d2, len(projectors)), complex)
 
     for i, (lbl, projector) in enumerate(projectors.items()):
-        flat_projector = projector.flatten()
+        flat_projector = projector.ravel()
         proj = _np.real_if_close(_np.vdot(flat_projector, flat_errorgen_std), tol=1000)
         if return_projected_errorgen:
             space_projector[:, i] = flat_projector

test/unit/construction/test_modelconstruction.py

@@ -43,8 +43,8 @@ def multikron(args):
         mdl1Q = smq1Q_XYI.target_model()
         mdlE0, mdlE1 = mdl1Q.povms['Mdefault']['0'].to_dense(), mdl1Q.povms['Mdefault']['1'].to_dense()
         ss1Q = pygsti.baseobjs.statespace.QubitSpace(1)
-        E0 = mc.create_spam_vector(0, ss1Q, 'pp').flatten()
-        E1 = mc.create_spam_vector(1, ss1Q, 'pp').flatten()
+        E0 = mc.create_spam_vector(0, ss1Q, 'pp').ravel()
+        E1 = mc.create_spam_vector(1, ss1Q, 'pp').ravel()
 
         self.assertArraysAlmostEqual(mdlE0, E0)
         self.assertArraysAlmostEqual(mdlE1, E1)
@@ -60,9 +60,10 @@ def multikron(args):
         for i in range(2**nQubits):
             bin_i = '{{0:0{}b}}'.format(nQubits).format(i)  # first .format creates format str, e.g. '{0:04b}'
             print("Testing state %d (%s)" % (i, bin_i))
-            created = mc.create_spam_vector(i, ssNQ, 'pp').flatten()
+            created = mc.create_spam_vector(i, ssNQ, 'pp').ravel()
             krond = multikron([E[digit] for digit in bin_i])
-            v = np.zeros(2**nQubits, complex); v[i] = 1.0
+            v = np.zeros(2**nQubits, complex)
+            v[i] = 1.0
             alt_created = st.create_from_pure_vector(v, 'static', 'pp', 'default', state_space=ssNQ).to_dense()
             self.assertArraysAlmostEqual(created, krond)
             self.assertArraysAlmostEqual(created, alt_created)

test/unit/tools/test_lindbladtools.py

@@ -85,6 +85,6 @@ def test_elementary_errorgen_bases(self):
             dot_mx = np.empty((len(duals), len(primals)), complex)
             for i, dual in enumerate(duals):
                 for j, primal in enumerate(primals):
-                    dot_mx[i,j] = np.vdot(dual.flatten(), primal.flatten())
+                    dot_mx[i,j] = np.vdot(dual, primal)
 
             self.assertTrue(np.allclose(dot_mx, np.identity(len(lbls), 'd')))