diff --git a/pygsti/modelmembers/operations/eigpdenseop.py b/pygsti/modelmembers/operations/eigpdenseop.py index 6173681fe..fd1cadf5b 100644 --- a/pygsti/modelmembers/operations/eigpdenseop.py +++ b/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 diff --git a/pygsti/modelmembers/operations/linearop.py b/pygsti/modelmembers/operations/linearop.py index b464498a0..6eb9bbdd6 100644 --- a/pygsti/modelmembers/operations/linearop.py +++ b/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): """ diff --git a/pygsti/modelmembers/povms/effect.py b/pygsti/modelmembers/povms/effect.py index d583cdc2e..c7bd387d6 100644 --- a/pygsti/modelmembers/povms/effect.py +++ b/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): """ diff --git a/pygsti/modelmembers/states/state.py b/pygsti/modelmembers/states/state.py index dd7d8225b..d072f4364 100644 --- a/pygsti/modelmembers/states/state.py +++ b/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): """ diff --git a/pygsti/report/workspaceplots.py b/pygsti/report/workspaceplots.py index 1e6a105f7..595357c4b 100644 --- a/pygsti/report/workspaceplots.py +++ b/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" diff --git a/pygsti/tools/jamiolkowski.py b/pygsti/tools/jamiolkowski.py index e204c5206..cafabb2d9 100644 --- a/pygsti/tools/jamiolkowski.py +++ b/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) diff --git a/pygsti/tools/lindbladtools.py b/pygsti/tools/lindbladtools.py index 9b24a9688..fcefc41ca 100644 --- a/pygsti/tools/lindbladtools.py +++ b/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 diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py index 0e176ca2e..9a32613a3 100644 --- a/pygsti/tools/matrixtools.py +++ b/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 diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 6b894e939..8d68a73ff 100644 --- a/pygsti/tools/optools.py +++ b/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 diff --git a/test/unit/construction/test_modelconstruction.py b/test/unit/construction/test_modelconstruction.py index 6e22cb481..77f2717e8 100644 --- a/test/unit/construction/test_modelconstruction.py +++ b/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) diff --git a/test/unit/tools/test_lindbladtools.py b/test/unit/tools/test_lindbladtools.py index 0870714e8..1be0230ab 100644 --- a/test/unit/tools/test_lindbladtools.py +++ b/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')))