Python numpy.trace() Examples

The following are 30 code examples for showing how to use numpy.trace(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: OpenFermion-Cirq   Author: quantumlib   File: optimal_givens_decomposition_test.py    License: Apache License 2.0 6 votes vote down vote up
def test_circuit_generation_and_accuracy():
    for dim in range(2, 10):
        qubits = cirq.LineQubit.range(dim)
        u_generator = numpy.random.random(
            (dim, dim)) + 1j * numpy.random.random((dim, dim))
        u_generator = u_generator - numpy.conj(u_generator).T
        assert numpy.allclose(-1 * u_generator, numpy.conj(u_generator).T)

        unitary = scipy.linalg.expm(u_generator)
        circuit = cirq.Circuit()
        circuit.append(optimal_givens_decomposition(qubits, unitary))

        fermion_generator = QubitOperator(()) * 0.0
        for i, j in product(range(dim), repeat=2):
            fermion_generator += jordan_wigner(
                FermionOperator(((i, 1), (j, 0)), u_generator[i, j]))

        true_unitary = scipy.linalg.expm(
            get_sparse_operator(fermion_generator).toarray())
        assert numpy.allclose(true_unitary.conj().T.dot(true_unitary),
                              numpy.eye(2 ** dim, dtype=complex))

        test_unitary = cirq.unitary(circuit)
        assert numpy.isclose(
            abs(numpy.trace(true_unitary.conj().T.dot(test_unitary))), 2 ** dim) 
Example 2
Project: pyscf   Author: pyscf   File: fciqmc.py    License: Apache License 2.0 6 votes vote down vote up
def read_neci_1dms(fciqmcci, norb, nelec, filename='OneRDM.1', directory='.'):
    ''' Read spinned rdms, as they are in the neci output '''

    f = open(os.path.join(directory, filename),'r')

    dm1a = numpy.zeros((norb,norb))
    dm1b = numpy.zeros((norb,norb))
    for line in f.readlines():
        linesp = line.split()
        i, j = int(linesp[0]), int(linesp[1])
        assert((i % 2) == (j % 2))
        if i % 2 == 1:
            # alpha
            assert(all(x<norb for x in (i/2,j/2)))
            dm1a[i/2,j/2] = float(linesp[2])
            dm1a[j/2,i/2] = float(linesp[2])
        else:
            assert(all(x<norb for x in (i/2 - 1,j/2 - 1)))
            dm1b[i/2 - 1,j/2 - 1] = float(linesp[2])
            dm1b[j/2 - 1,i/2 - 1] = float(linesp[2])

    f.close()
    assert(numpy.allclose(dm1a.trace()+dm1b.trace(),sum(nelec)))
    return dm1a, dm1b 
Example 3
Project: pyGSTi   Author: pyGSTio   File: basis.py    License: Apache License 2.0 6 votes vote down vote up
def is_normalized(self):
        '''
        Check if a basis is normalized, meaning that Tr(Bi Bi) = 1.0.

        Available only to bases whose elements are *matrices* for now.

        Returns
        -------
        bool
        '''
        if self.elndim == 2:
            for i, mx in enumerate(self.elements):
                t = _np.trace(_np.dot(mx, mx))
                t = _np.real(t)
                if not _np.isclose(t, 1.0): return False
            return True
        elif self.elndim == 1:
            raise NotImplementedError("TODO: add code so this works for *vector*-valued bases too!")
        else:
            raise ValueError("I don't know what normalized means for elndim == %d!" % self.elndim) 
Example 4
Project: pyGSTi   Author: pyGSTio   File: optools.py    License: Apache License 2.0 6 votes vote down vote up
def tracenorm(A):
    """
    Compute the trace norm of matrix A given by:

      Tr( sqrt{ A^dagger * A } )

    Parameters
    ----------
    A : numpy array
        The matrix to compute the trace norm of.
    """
    if _np.linalg.norm(A - _np.conjugate(A.T)) < 1e-8:
        #Hermitian, so just sum eigenvalue magnitudes
        return _np.sum(_np.abs(_np.linalg.eigvals(A)))
    else:
        #Sum of singular values (positive by construction)
        return _np.sum(_np.linalg.svd(A, compute_uv=False)) 
Example 5
Project: pyGSTi   Author: pyGSTio   File: optools.py    License: Apache License 2.0 6 votes vote down vote up
def jtracedist(A, B, mxBasis='pp'):  # Jamiolkowski trace distance:  Tr(|J(A)-J(B)|)
    """
    Compute the Jamiolkowski trace distance between operation matrices A and B,
    given by:

      D = 0.5 * Tr( sqrt{ (J(A)-J(B))^2 } )

    where J(.) is the Jamiolkowski isomorphism map that maps a operation matrix
    to it's corresponding Choi Matrix.

    Parameters
    ----------
    A, B : numpy array
        The matrices to compute the distance between.

    mxBasis : {'std', 'gm', 'pp', 'qt'} or Basis object
        The source and destination basis, respectively.  Allowed
        values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp),
        and Qutrit (qt) (or a custom basis object).
    """
    JA = _jam.fast_jamiolkowski_iso_std(A, mxBasis)
    JB = _jam.fast_jamiolkowski_iso_std(B, mxBasis)
    return tracedist(JA, JB) 
Example 6
Project: pyGSTi   Author: pyGSTio   File: optools.py    License: Apache License 2.0 6 votes vote down vote up
def povm_jtracedist(model, targetModel, povmlbl):
    """
    Computes the Jamiolkowski trace distance between POVM maps using :func:`jtracedist`.

    Parameters
    ----------
    model, targetModel : Model
        Models containing the two POVMs to compare.

    povmlbl : str
        The POVM label

    Returns
    -------
    float
    """
    povm_mx = get_povm_map(model, povmlbl)
    target_povm_mx = get_povm_map(targetModel, povmlbl)
    return jtracedist(povm_mx, target_povm_mx, targetModel.basis) 
Example 7
Project: rai-python   Author: MarcToussaint   File: transformations.py    License: MIT License 6 votes vote down vote up
def reflection_matrix(point, normal):
    """Return matrix to mirror at plane defined by point and normal vector.
    >>> v0 = numpy.random.random(4) - 0.5
    >>> v0[3] = 1.0
    >>> v1 = numpy.random.random(3) - 0.5
    >>> R = reflection_matrix(v0, v1)
    >>> numpy.allclose(2., numpy.trace(R))
    True
    >>> numpy.allclose(v0, numpy.dot(R, v0))
    True
    >>> v2 = v0.copy()
    >>> v2[:3] += v1
    >>> v3 = v0.copy()
    >>> v2[:3] -= v1
    >>> numpy.allclose(v2, numpy.dot(R, v3))
    True
    """
    normal = unit_vector(normal[:3])
    M = numpy.identity(4)
    M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
    M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
    return M 
Example 8
Project: connecting_the_dots   Author: autonomousvision   File: geometry.py    License: MIT License 6 votes vote down vote up
def axisangle_from_rotm(R):
  # logarithm of rotation matrix
  # R = R.reshape(-1,3,3)
  # tr = np.trace(R, axis1=1, axis2=2)
  # phi = np.arccos(np.clip((tr - 1) / 2, -1, 1))
  # scale = np.zeros_like(phi)
  # div = 2 * np.sin(phi)
  # np.divide(phi, div, out=scale, where=np.abs(div) > 1e-6)
  # A = (R - R.transpose(0,2,1)) * scale.reshape(-1,1,1)
  # aa = np.stack((A[:,2,1], A[:,0,2], A[:,1,0]), axis=1)
  # return aa.squeeze()
  R = R.reshape(-1,3,3)
  omega = np.empty((R.shape[0], 3), dtype=R.dtype)
  omega[:,0] = R[:,2,1] - R[:,1,2]
  omega[:,1] = R[:,0,2] - R[:,2,0]
  omega[:,2] = R[:,1,0] - R[:,0,1]
  r = np.linalg.norm(omega, axis=1).reshape(-1,1)
  t = np.trace(R, axis1=1, axis2=2).reshape(-1,1)
  omega = np.arctan2(r, t-1) * omega
  aa = np.zeros_like(omega)
  np.divide(omega, r, out=aa, where=r != 0)
  return aa.squeeze() 
Example 9
Project: cozmo_driver   Author: OTL   File: transformations.py    License: Apache License 2.0 6 votes vote down vote up
def reflection_matrix(point, normal):
    """Return matrix to mirror at plane defined by point and normal vector.

    >>> v0 = numpy.random.random(4) - 0.5
    >>> v0[3] = 1.0
    >>> v1 = numpy.random.random(3) - 0.5
    >>> R = reflection_matrix(v0, v1)
    >>> numpy.allclose(2., numpy.trace(R))
    True
    >>> numpy.allclose(v0, numpy.dot(R, v0))
    True
    >>> v2 = v0.copy()
    >>> v2[:3] += v1
    >>> v3 = v0.copy()
    >>> v2[:3] -= v1
    >>> numpy.allclose(v2, numpy.dot(R, v3))
    True

    """
    normal = unit_vector(normal[:3])
    M = numpy.identity(4)
    M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
    M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
    return M 
Example 10
Project: SfmLearner-Pytorch   Author: ClementPinard   File: test_pose.py    License: MIT License 6 votes vote down vote up
def compute_pose_error(gt, pred):
    RE = 0
    snippet_length = gt.shape[0]
    scale_factor = np.sum(gt[:,:,-1] * pred[:,:,-1])/np.sum(pred[:,:,-1] ** 2)
    ATE = np.linalg.norm((gt[:,:,-1] - scale_factor * pred[:,:,-1]).reshape(-1))
    for gt_pose, pred_pose in zip(gt, pred):
        # Residual matrix to which we compute angle's sin and cos
        R = gt_pose[:,:3] @ np.linalg.inv(pred_pose[:,:3])
        s = np.linalg.norm([R[0,1]-R[1,0],
                            R[1,2]-R[2,1],
                            R[0,2]-R[2,0]])
        c = np.trace(R) - 1
        # Note: we actually compute double of cos and sin, but arctan2 is invariant to scale
        RE += np.arctan2(s,c)

    return ATE/snippet_length, RE/snippet_length 
Example 11
Project: bop_toolkit   Author: thodan   File: transform.py    License: MIT License 6 votes vote down vote up
def reflection_matrix(point, normal):
  """Return matrix to mirror at plane defined by point and normal vector.

  >>> v0 = numpy.random.random(4) - 0.5
  >>> v0[3] = 1.
  >>> v1 = numpy.random.random(3) - 0.5
  >>> R = reflection_matrix(v0, v1)
  >>> numpy.allclose(2, numpy.trace(R))
  True
  >>> numpy.allclose(v0, numpy.dot(R, v0))
  True
  >>> v2 = v0.copy()
  >>> v2[:3] += v1
  >>> v3 = v0.copy()
  >>> v2[:3] -= v1
  >>> numpy.allclose(v2, numpy.dot(R, v3))
  True

  """
  normal = unit_vector(normal[:3])
  M = numpy.identity(4)
  M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
  M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
  return M 
Example 12
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    License: MIT License 6 votes vote down vote up
def test_trace():
    rng = numpy.random.RandomState(utt.fetch_seed())
    x = theano.tensor.matrix()
    g = trace(x)
    f = theano.function([x], g)

    for shp in [(2, 3), (3, 2), (3, 3)]:
        m = rng.rand(*shp).astype(config.floatX)
        v = numpy.trace(m)
        assert v == f(m)

    xx = theano.tensor.vector()
    ok = False
    try:
        trace(xx)
    except TypeError:
        ok = True
    assert ok 
Example 13
Project: OpenFermion-Cirq   Author: quantumlib   File: higham_test.py    License: Apache License 2.0 5 votes vote down vote up
def test_reconstruction():
    dim = 20
    np.random.seed(42)
    mat = np.random.random((dim, dim))
    mat = 0.5 * (mat + mat.T)
    test_mat = np.zeros_like(mat)
    w, v = np.linalg.eigh(mat)
    for i in range(w.shape[0]):
        test_mat += w[i] * v[:, [i]].dot(v[:, [i]].T)
    assert np.allclose(test_mat - mat, 0.0)

    test_mat = fixed_trace_positive_projection(mat, np.trace(mat))
    assert np.isclose(np.trace(test_mat), np.trace(mat))
    w, v = np.linalg.eigh(test_mat)
    assert np.all(w >= -(float(4.0E-15)))

    mat = np.arange(16).reshape((4, 4))
    mat = 0.5 * (mat + mat.T)
    mat_tensor = map_to_tensor(mat)
    trace_mat = np.trace(mat)
    true_mat = fixed_trace_positive_projection(mat, trace_mat)
    test_mat = map_to_matrix(fixed_trace_positive_projection(mat_tensor,
                                                             trace_mat))
    assert np.allclose(true_mat, test_mat)

    w, v = np.linalg.eigh(true_mat)
    assert np.allclose(true_mat, fixed_trace_positive_projection(true_mat,
                                                                 trace_mat)) 
Example 14
Project: pyscf   Author: pyscf   File: test_kpoint.py    License: Apache License 2.0 5 votes vote down vote up
def test_rdm1(self):
        cell = pbcgto.Cell()
        cell.atom = '''Al 0 0 0'''
        cell.basis = 'gth-szv'
        cell.pseudo = 'gth-pade'
        cell.a = '''
        2.47332919 0 1.42797728
        0.82444306 2.33187713 1.42797728
        0, 0, 2.85595455
        '''
        cell.unit = 'angstrom'
        cell.build()
        cell.verbose = 4
        cell.incore_anyway = True

        abs_kpts = cell.make_kpts((2, 1, 1), scaled_center=(.1, .2, .3))
        kmf = pbcscf.KRHF(cell, abs_kpts)
        kmf.conv_tol = 1e-12
        kmf.kernel()
        mp = pyscf.pbc.mp.kmp2.KMP2(kmf)
        mp.kernel(with_t2=True)
        self.assertAlmostEqual(mp.e_corr, -0.00162057921874043)
        dm = mp.make_rdm1()
        np.testing.assert_allclose(np.trace(dm[0]) + np.trace(dm[1]), 6)
        for kdm in dm:
            np.testing.assert_allclose(kdm, kdm.conj().T) 
Example 15
Project: pyscf   Author: pyscf   File: test_uhf.py    License: Apache License 2.0 5 votes vote down vote up
def test_det_ovlp(self):
        s, x = mf.det_ovlp(mf.mo_coeff, mf.mo_coeff, mf.mo_occ, mf.mo_occ)
        self.assertAlmostEqual(s, 1.000000000, 9)
        self.assertAlmostEqual(numpy.trace(x[0]), mol.nelec[0]*1.000000000, 9)
        self.assertAlmostEqual(numpy.trace(x[0]), mol.nelec[1]*1.000000000, 9) 
Example 16
Project: pyscf   Author: pyscf   File: tdscf.py    License: Apache License 2.0 5 votes vote down vote up
def trdot(A,B):
    C = np.trace(np.dot(A,B))
    return C 
Example 17
Project: pyscf   Author: pyscf   File: tdscf.py    License: Apache License 2.0 5 votes vote down vote up
def energy(self,dm_lao,fmat,jmat,kmat):
        """
        Args:
            dm_lao: complex
                Density in LAO basis.
            fmat: complex
                Fock matrix in Lowdin AO basis
            jmat: complex
                Coulomb matrix in AO basis
            kmat: complex
                Exact Exchange in AO basis
        Returns:
            e_tot: float
                Total Energy of a system
        """
        if (self.params["Model"] == "TDHF"):
            hlao = transmat(self.h,self.x)
            e_tot = (self.enuc+np.trace(np.dot(dm_lao,hlao+fmat))).real
            return e_tot
        elif self.params["Model"] == "TDDFT":
            dm = transmat(dm_lao,self.x,-1)
            exc = self.exc
            if(self.hyb > 0.01):
                exc -= 0.5 * self.hyb * trdot(dm,kmat)
            # if not using auxmol
            eh = trdot(dm,2*self.h)
            ej = trdot(dm,jmat)
            e_tot = (eh + ej + exc + self.enuc).real
            return e_tot 
Example 18
Project: pyscf   Author: pyscf   File: tdscf.py    License: Apache License 2.0 5 votes vote down vote up
def loginstant(self, rho, c_am, v_lm, fmat, jmat, kmat, tnow, it):
        """
        time is logged in atomic units.
        Args:
            rho: complex
                MO density matrix.
            c_am: complex
                Transformation Matrix |AO><MO|
            v_lm: complex
                Transformation Matrix |LAO><MO|
            fmat: complex
                Fock matrix in Lowdin AO basis
            jmat: complex
                Coulomb matrix in AO basis
            kmat: complex
                Exact Exchange in AO basis
            tnow: float
                Current time in propagation in A.U.
            it: int
                Number of iteration of propagation
        Returns:
            tore: str
                |t, dipole(x,y,z), energy|

        """
        np.set_printoptions(precision = 7)
        tore = str(tnow)+" "+str(self.dipole(rho, c_am).real).rstrip("]").lstrip("[")+\
         " " +str(self.energy(transmat(rho,v_lm,-1),fmat, jmat, kmat))

        if it%self.params["StatusEvery"] ==0 or it == self.params["MaxIter"]-1:
            logger.log(self, "t: %f fs    Energy: %f a.u.   Total Density: %f",\
            tnow*FSPERAU,self.energy(transmat(rho,v_lm,-1),fmat, jmat, kmat), \
            2*np.trace(rho))
            logger.log(self, "Dipole moment(X, Y, Z, au): %8.5f, %8.5f, %8.5f",\
             self.dipole(rho, c_am).real[0],self.dipole(rho, c_am).real[1],\
             self.dipole(rho, c_am).real[2])
        return tore 
Example 19
Project: pyGSTi   Author: pyGSTio   File: spamvec.py    License: Apache License 2.0 5 votes vote down vote up
def __init__(self, vec, basis, truncate=False):
        """
        Initialize a CPTPSPAMVec object.

        Parameters
        ----------
        vec : array_like or SPAMVec
            a 1D numpy array representing the SPAM operation.  The
            shape of this array sets the dimension of the SPAM op.

        basis : {"std", "gm", "pp", "qt"} or Basis
            The basis `vec` is in.  Needed because this parameterization
            requires we construct the density matrix corresponding to
            the Lioville vector `vec`.

        trunctate : bool, optional
            Whether or not a non-positive, trace=1 `vec` should
            be truncated to force a successful construction.
        """
        vector = SPAMVec.convert_to_vector(vec)
        basis = _Basis.cast(basis, len(vector))

        self.basis = basis
        self.basis_mxs = basis.elements  # shape (len(vec), dmDim, dmDim)
        self.basis_mxs = _np.rollaxis(self.basis_mxs, 0, 3)  # shape (dmDim, dmDim, len(vec))
        assert(self.basis_mxs.shape[-1] == len(vector))

        # set self.params and self.dmDim
        self._set_params_from_vector(vector, truncate)

        #scratch space
        self.Lmx = _np.zeros((self.dmDim, self.dmDim), 'complex')

        DenseSPAMVec.__init__(self, vector, "densitymx", "prep") 
Example 20
Project: pyGSTi   Author: pyGSTio   File: spamvec.py    License: Apache License 2.0 5 votes vote down vote up
def _construct_vector(self):
        dmDim = self.dmDim

        #  params is an array of length dmDim^2-1 that
        #  encodes a lower-triangular matrix "Lmx" via:
        #  Lmx[i,i] = params[i*dmDim + i] / param-norm  # i = 0...dmDim-2
        #     *last diagonal el is given by sqrt(1.0 - sum(L[i,j]**2))
        #  Lmx[i,j] = params[i*dmDim + j] + 1j*params[j*dmDim+i] (i > j)

        param2Sum = _np.vdot(self.params, self.params)  # or "dot" would work, since params are real
        paramNorm = _np.sqrt(param2Sum)  # also the norm of *all* Lmx els

        for i in range(dmDim):
            self.Lmx[i, i] = self.params[i * dmDim + i] / paramNorm
            for j in range(i):
                self.Lmx[i, j] = (self.params[i * dmDim + j] + 1j * self.params[j * dmDim + i]) / paramNorm

        Lmx_norm = _np.trace(_np.dot(self.Lmx.T.conjugate(), self.Lmx))  # sum of magnitude^2 of all els
        assert(_np.isclose(Lmx_norm, 1.0)), "Violated trace=1 condition!"

        #The (complex, Hermitian) density matrix is build by
        # assuming Lmx is its Cholesky decomp, which makes
        # the density matrix is pos-def.
        density_mx = _np.dot(self.Lmx, self.Lmx.T.conjugate())
        assert(_np.isclose(_np.trace(density_mx), 1.0)), "density matrix must be trace == 1"

        # write density matrix in given basis: = sum_i alpha_i B_i
        # ASSUME that basis is orthogonal, i.e. Tr(Bi^dag*Bj) = delta_ij
        basis_mxs = _np.rollaxis(self.basis_mxs, 2)  # shape (dmDim, dmDim, len(vec))
        vec = _np.array([_np.trace(_np.dot(M.T.conjugate(), density_mx)) for M in basis_mxs])

        #for now, assume Liouville vector should always be real (TODO: add 'real' flag later?)
        assert(_np.linalg.norm(_np.imag(vec)) < IMAG_TOL)
        vec = _np.real(vec)

        self.base1D.flags.writeable = True
        self.base1D[:] = vec[:]  # so shape is (dim,1) - the convention for spam vectors
        self.base1D.flags.writeable = False 
Example 21
Project: pyGSTi   Author: pyGSTio   File: reportables.py    License: Apache License 2.0 5 votes vote down vote up
def choi_trace(gate, mxBasis):
    """Trace of the Choi matrix"""
    choi = _tools.jamiolkowski_iso(gate, mxBasis, mxBasis)
    return _np.trace(choi) 
Example 22
Project: pyGSTi   Author: pyGSTio   File: reportables.py    License: Apache License 2.0 5 votes vote down vote up
def evaluate_nearby(self, nearby_model):
            """Evaluate at a nearby model"""
            mxBasis = nearby_model.basis
            JAstd = self.d * _tools.fast_jamiolkowski_iso_std(
                nearby_model.product(self.circuit), mxBasis)
            JBstd = self.d * _tools.fast_jamiolkowski_iso_std(self.B, mxBasis)
            Jt = (JBstd - JAstd).T
            return 0.5 * _np.trace(_np.dot(Jt.real, self.W.real) + _np.dot(Jt.imag, self.W.imag))

    #def circuit_half_diamond_norm(modelA, modelB, circuit):
    #    A = modelA.product(circuit) # "gate"
    #    B = modelB.product(circuit) # "target gate"
    #    return half_diamond_norm(A, B, modelB.basis)
    #Circuit_half_diamond_norm = _modf.modelfn_factory(circuit_half_diamond_norm)
    #  # init args == (modelA, modelB, circuit) 
Example 23
Project: pyGSTi   Author: pyGSTio   File: reportables.py    License: Apache License 2.0 5 votes vote down vote up
def povm_jt_diff(modelA, modelB, povmlbl):
    """
    POVM Jamiolkowski trace distance between `modelA` and `modelB`, equal to
    `Jamiolkowski_trace_distance(POVM_MAP)` where `POVM_MAP` is the extension
    of the POVM from the classical space of k-outcomes to the space of
    (diagonal) k by k density matrices.
    """
    return _tools.povm_jtracedist(modelA, modelB, povmlbl) 
Example 24
Project: pyGSTi   Author: pyGSTio   File: reportables.py    License: Apache License 2.0 5 votes vote down vote up
def maximum_trace_dist(gate, mxBasis):
    """ Jamiolkowski trace distance between `gate` and its closest unitary"""
    closestUOpMx = _alg.find_closest_unitary_opmx(gate)
    #closestUJMx = _tools.jamiolkowski_iso(closestUOpMx, mxBasis, mxBasis)
    _tools.jamiolkowski_iso(closestUOpMx, mxBasis, mxBasis)
    return _tools.jtracedist(gate, closestUOpMx) 
Example 25
Project: pyGSTi   Author: pyGSTio   File: reportables.py    License: Apache License 2.0 5 votes vote down vote up
def jt_diff(A, B, mxBasis):  # assume vary model1, model2 fixed
    """ Jamiolkowski trace distance between A and B"""
    return _tools.jtracedist(A, B, mxBasis) 
Example 26
Project: pyGSTi   Author: pyGSTio   File: reportables.py    License: Apache License 2.0 5 votes vote down vote up
def evaluate_nearby(self, nearby_model):
            """Evaluates at a nearby model"""
            gl = self.oplabel; mxBasis = nearby_model.basis
            JAstd = self.d * _tools.fast_jamiolkowski_iso_std(
                nearby_model.operations[gl].todense(), mxBasis)
            JBstd = self.d * _tools.fast_jamiolkowski_iso_std(self.B, mxBasis)
            Jt = (JBstd - JAstd).T
            return 0.5 * _np.trace(_np.dot(Jt.real, self.W.real) + _np.dot(Jt.imag, self.W.imag)) 
Example 27
Project: pyGSTi   Author: pyGSTio   File: optools.py    License: Apache License 2.0 5 votes vote down vote up
def tracedist(A, B):
    """
    Compute the trace distance between matrices A and B,
    given by:

      D = 0.5 * Tr( sqrt{ (A-B)^dagger * (A-B) } )

    Parameters
    ----------
    A, B : numpy array
        The matrices to compute the distance between.
    """
    return 0.5 * tracenorm(A - B) 
Example 28
Project: pyGSTi   Author: pyGSTio   File: optools.py    License: Apache License 2.0 5 votes vote down vote up
def entanglement_fidelity(A, B, mxBasis='pp'):
    """
    Returns the "entanglement" process fidelity between gate
    matrices A and B given by :

      F = Tr( sqrt{ sqrt(J(A)) * J(B) * sqrt(J(A)) } )^2

    where J(.) is the Jamiolkowski isomorphism map that maps a operation matrix
    to it's corresponding Choi Matrix.

    Parameters
    ----------
    A, B : numpy array
        The matrices to compute the fidelity between.

    mxBasis : {'std', 'gm', 'pp', 'qt'} or Basis object
        The basis of the matrices.  Allowed values are Matrix-unit (std),
        Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt)
        (or a custom basis object).
    """
    d2 = A.shape[0]
    def isTP(x): return _np.isclose(x[0, 0], 1.0) and all(
        [_np.isclose(x[0, i], 0) for i in range(d2)])

    def isUnitary(x): return _np.allclose(_np.identity(d2, 'd'), _np.dot(x, x.conjugate().T))

    if isTP(A) and isTP(B) and isUnitary(B):  # then assume TP-like gates & use simpler formula
        TrLambda = _np.trace(_np.dot(A, B.conjugate().T))  # same as using _np.linalg.inv(B)
        d2 = A.shape[0]
        return TrLambda / d2

    JA = _jam.jamiolkowski_iso(A, mxBasis, mxBasis)
    JB = _jam.jamiolkowski_iso(B, mxBasis, mxBasis)
    return fidelity(JA, JB) 
Example 29
Project: pyGSTi   Author: pyGSTio   File: germselection.py    License: Apache License 2.0 5 votes vote down vote up
def _SuperOpForPerfectTwirl(wrt, eps):
    """Return super operator for doing a perfect twirl with respect to wrt.
    """
    assert wrt.shape[0] == wrt.shape[1]  # only square matrices allowed
    dim = wrt.shape[0]
    SuperOp = _np.zeros((dim**2, dim**2), 'complex')

    # Get spectrum and eigenvectors of wrt
    wrtEvals, wrtEvecs = _np.linalg.eig(wrt)
    wrtEvecsInv = _np.linalg.inv(wrtEvecs)

    # We want to project  X -> M * (Proj_i * (Minv * X * M) * Proj_i) * Minv,
    # where M = wrtEvecs. So A = B = M * Proj_i * Minv and so
    # superop = A tensor B^T == A tensor A^T
    # NOTE: this == (A^T tensor A)^T while *Maple* germ functions seem to just
    # use A^T tensor A -> ^T difference
    for i in range(dim):
        # Create projector onto i-th eigenspace (spanned by i-th eigenvector
        # and other degenerate eigenvectors)
        Proj_i = _np.diag([(1 if (abs(wrtEvals[i] - wrtEvals[j]) <= eps)
                            else 0) for j in range(dim)])
        A = _np.dot(wrtEvecs, _np.dot(Proj_i, wrtEvecsInv))
        #if _np.linalg.norm(A.imag) > 1e-6:
        #    print("DB: imag = ",_np.linalg.norm(A.imag))
        #assert(_np.linalg.norm(A.imag) < 1e-6)
        #A = _np.real(A)
        # Need to normalize, because we are overcounting projectors onto
        # subspaces of dimension d > 1, giving us d * Proj_i tensor Proj_i^T.
        # We can fix this with a division by tr(Proj_i) = d.
        SuperOp += _np.kron(A, A.T) / _np.trace(Proj_i)
        # SuperOp += _np.kron(A.T,A) # Mimic Maple version (but I think this is
        # wrong... or it doesn't matter?)
    return SuperOp  # a op_dim^2 x op_dim^2 matrix 
Example 30
Project: pyGSTi   Author: pyGSTio   File: test_jamiolkowski.py    License: Apache License 2.0 5 votes vote down vote up
def checkBasis(self, cmb):
        #Op with Jamio map on gate in std and gm bases
        Jmx1 = j.jamiolkowski_iso(self.testGate, opMxBasis=self.stdSmall,
                                  choiMxBasis=cmb)
        Jmx2 = j.jamiolkowski_iso(self.testGateGM_mx, opMxBasis=self.gmSmall,
                                  choiMxBasis=cmb)
        print("Jmx1.shape = ", Jmx1.shape)

        #Make sure these yield the same trace == 1 matrix
        self.assertArraysAlmostEqual(Jmx1, Jmx2)
        self.assertAlmostEqual(np.trace(Jmx1), 1.0)

        #Op on expanded gate in std and gm bases
        JmxExp1 = j.jamiolkowski_iso(self.expTestGate_mx, opMxBasis=self.std, choiMxBasis=cmb)
        JmxExp2 = j.jamiolkowski_iso(self.expTestGateGM_mx, opMxBasis=self.gm, choiMxBasis=cmb)
        print("JmxExp1.shape = ", JmxExp1.shape)

        #Make sure these are the same as operating on the contracted basis
        self.assertArraysAlmostEqual(Jmx1, JmxExp1)
        self.assertArraysAlmostEqual(Jmx1, JmxExp2)

        #Reverse transform should yield back the operation matrix
        revTestGate_mx = j.jamiolkowski_iso_inv(Jmx1, choiMxBasis=cmb,
                                                opMxBasis=self.gmSmall)
        self.assertArraysAlmostEqual(revTestGate_mx, self.testGateGM_mx)

        #Reverse transform without specifying stateSpaceDims, then contraction, should yield same result
        revExpTestGate_mx = j.jamiolkowski_iso_inv(Jmx1, choiMxBasis=cmb, opMxBasis=self.std)
        self.assertArraysAlmostEqual(bt.resize_std_mx(revExpTestGate_mx, 'contract', self.std, self.stdSmall),
                                     self.testGate)