numpy.einsum

Here are the examples of the python api numpy.einsum taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

200 Examples 7

Example 1

Project: pyscf
Source File: ccsd_slow.py
View license
def update_amp_t1(cc, t1, t2, eris, foo, fov, fvv):
    nocc, nvir = t1.shape
    fock = eris.fock

    g2 = 2 * eris.ovov - eris.oovv.transpose(1,2,0,3)
    t1new = fock[:nocc,nocc:] \
          + numpy.einsum('ib,ab->ia', t1, fvv) \
          - numpy.einsum('ja,ji->ia', t1, foo) \
          + numpy.einsum('jb,iajb->ia', t1, g2)
    theta = t2 * 2 - t2.transpose(2,1,0,3)
    t1new += numpy.einsum('jb,iajb->ia', fov, theta)
    t1new += numpy.einsum('ibjc,jcba->ia', theta, eris.ovvv)
    t1new -= numpy.einsum('jbki,jbka->ia', eris.ovoo, theta)
    mo_e = fock.diagonal()
    eia = mo_e[:nocc,None] - mo_e[None,nocc:]
    t1new /= eia
    return t1new

Example 2

Project: pyscf
Source File: ccsd_slow.py
View license
def energy(cc, t1, t2, eris):
    nocc, nvir = t1.shape
    tau = t2 + numpy.einsum('ia,jb->iajb', t1, t1)
    theta = tau * 2 - tau.transpose(2,1,0,3)
    fock = eris.fock
    e = numpy.einsum('ia,ia', fock[:nocc,nocc:], t1) * 2 \
      + numpy.einsum('iajb,iajb', eris.ovov, theta)
    return e

Example 3

Project: pyscf
Source File: 21-concatenate_molecules.py
View license
def get_vhf(mol, dm, *args, **kwargs):
    naux = eri2c.shape[0]
    rho = numpy.einsum('ijp,ij->p', eri3c, dm)
    rho = numpy.linalg.solve(eri2c, rho)
    jmat = numpy.einsum('p,ijp->ij', rho, eri3c)
    kpj = numpy.einsum('ijp,jk->ikp', eri3c, dm)
    pik = numpy.linalg.solve(eri2c, kpj.reshape(-1,naux).T).reshape(-1,nao,nao)
    kmat = numpy.einsum('pik,kjp->ij', pik, eri3c)
    return jmat - kmat * .5

Example 4

Project: pyscf
Source File: icmpspt.py
View license
def makeheff(int1, int2popo, int2ppoo, E1, ncore, nvirt, frozen):
        nc = int1.shape[0]-nvirt

        int1_eff = 1.*int1 + 2.0*numpy.einsum('mnii->mn', int2ppoo[:, :, frozen:ncore, frozen:ncore])-numpy.einsum('mini->mn', int2popo[:,frozen:ncore,:, frozen:ncore])

        int1_eff[:ncore, :ncore] += numpy.einsum('lmjk,jk->lm',int2ppoo[:ncore,:ncore,ncore:nc,ncore:nc], E1) - 0.5*numpy.einsum('ljmk,jk->lm',int2popo[:ncore,ncore:nc,:ncore,ncore:nc], E1)
        int1_eff[nc:, nc:] += numpy.einsum('lmjk,jk->lm',int2ppoo[nc:,nc:,ncore:nc,ncore:nc], E1) - 0.5*numpy.einsum('ljmk,jk->lm',int2popo[nc:,ncore:nc,nc:,ncore:nc], E1)
        #int1_eff[nc:, nc:] += numpy.einsum('ljmk,jk->lm',int2[nc:,ncore:nc,nc:,ncore:nc], E1) - 0.5*numpy.einsum('ljkm,jk->lm',int2[nc:,ncore:nc,ncore:nc,nc:], E1)
        return int1_eff

Example 5

Project: pyscf
Source File: nevpt2.py
View license
def make_a17(h1e,h2e,dm2,dm3):
    h1e = h1e - numpy.einsum('mjjn->mn',h2e)

    a17 = -numpy.einsum('pi,cabi->abcp',h1e,dm2)\
          -numpy.einsum('kpij,cabjki->abcp',h2e,dm3)
    return a17

Example 6

Project: pyscf
Source File: nevpt2.py
View license
def make_a19(h1e,h2e,dm1,dm2):
    h1e = h1e - numpy.einsum('mjjn->mn',h2e)

    a19 = -numpy.einsum('pi,ai->ap',h1e,dm1)\
          -numpy.einsum('kpij,ajki->ap',h2e,dm2)
    return a19

Example 7

Project: pyscf
Source File: nevpt2.py
View license
def make_a23(h1e,h2e,dm1,dm2,dm3):
    a23 = -numpy.einsum('ip,caib->abcp',h1e,dm2)\
          -numpy.einsum('pijk,cajbik->abcp',h2e,dm3)\
          +2.0*numpy.einsum('bp,ca->abcp',h1e,dm1)\
          +2.0*numpy.einsum('pibk,caik->abcp',h2e,dm2)

    return a23

Example 8

Project: pyscf
Source File: nevpt2.py
View license
def make_a25(h1e,h2e,dm1,dm2):

    a25 = -numpy.einsum('pi,ai->ap',h1e,dm1)\
          -numpy.einsum('pijk,jaik->ap',h2e,dm2)\
          +2.0*numpy.einsum('ap->pa',h1e)\
          +2.0*numpy.einsum('piaj,ij->ap',h2e,dm1)

    return a25

Example 9

Project: pyscf
Source File: nevpt2.py
View license
def make_a9(h1e,h2e,hdm1,hdm2,hdm3):
    a9 =  numpy.einsum('ib,pqai->pqab',h1e,hdm2)
    a9 += numpy.einsum('ijib,pqaj->pqab',h2e,hdm2)*2.0
    a9 -= numpy.einsum('ijjb,pqai->pqab',h2e,hdm2)
    a9 -= numpy.einsum('ijkb,pkqaij->pqab',h2e,hdm3)
    a9 += numpy.einsum('ia,pqib->pqab',h1e,hdm2)
    a9 -= numpy.einsum('ijja,pqib->pqab',h2e,hdm2)
    a9 -= numpy.einsum('ijba,pqji->pqab',h2e,hdm2)
    a9 += numpy.einsum('ijia,pqjb->pqab',h2e,hdm2)*2.0
    a9 -= numpy.einsum('ijka,pqkjbi->pqab',h2e,hdm3)
    return a9

Example 10

Project: pyscf
Source File: nevpt2.py
View license
def make_a12(h1e,h2e,dm1,dm2,dm3):
    a12 = numpy.einsum('ia,qpib->pqab',h1e,dm2)\
        - numpy.einsum('bi,qpai->pqab',h1e,dm2)\
        + numpy.einsum('ijka,qpjbik->pqab',h2e,dm3)\
        - numpy.einsum('kbij,qpajki->pqab',h2e,dm3)\
        - numpy.einsum('bjka,qpjk->pqab',h2e,dm2)\
        + numpy.einsum('jbij,qpai->pqab',h2e,dm2)
    return a12

Example 11

Project: chainer
Source File: test_bilinear.py
View license
    def setUp(self):
        super(TestBilinear2, self).setUp()

        assert self.in_shape[1] % 2 == 0
        self.e1 = _uniform(self.batch_size, 1, self.in_shape[0])
        self.e2 = _uniform(self.batch_size, self.in_shape[1] // 2, 2)
        self.gy = _uniform(self.batch_size, self.out_size)

        e1 = array.as_mat(self.e1)
        e2 = array.as_mat(self.e2)

        self.y = (
            numpy.einsum('ij,ik,jkl->il', e1, e2, self.W) +
            e1.dot(self.V1) + e2.dot(self.V2) + self.b)

Example 12

Project: chainer
Source File: test_bilinear.py
View license
    def setUp(self):
        super(TestBilinearWOBias2, self).setUp()

        assert self.in_shape[1] % 2 == 0
        self.e1 = _uniform(self.batch_size, 1, self.in_shape[0])
        self.e2 = _uniform(self.batch_size, 2, self.in_shape[1] // 2)
        self.gy = _uniform(self.batch_size, self.out_size)

        e1 = array.as_mat(self.e1)
        e2 = array.as_mat(self.e2)

        self.y = numpy.einsum('ij,ik,jkl->il', e1, e2, self.W)

Example 13

Project: attention-lvcsr
Source File: test_blocksparse.py
View license
    @staticmethod
    def gemv_numpy3(o, W, h, iIdx, oIdx):
        """
        Other implementation
        """
        from numpy import ix_
        for b in range(o.shape[0]):
            w = W[ix_(iIdx[b], oIdx[b])]
            # The next three lines do the same operation. The last one is the
            # fastest
            # o[b] += (h[b][:, None, :, None] * w).sum(axis=(0, 2))
            # o[b] += numpy.tensordot(h[b], w, [(0,1),(0,2)])
            o[b] += numpy.einsum('ik,ijkl', h[b], w)
        return o

Example 14

Project: scikit-learn
Source File: extmath.py
View license
def row_norms(X, squared=False):
    """Row-wise (squared) Euclidean norm of X.

    Equivalent to np.sqrt((X * X).sum(axis=1)), but also supports sparse
    matrices and does not create an X.shape-sized temporary.

    Performs no input validation.
    """
    if issparse(X):
        if not isinstance(X, csr_matrix):
            X = csr_matrix(X)
        norms = csr_row_norms(X)
    else:
        norms = np.einsum('ij,ij->i', X, X)

    if not squared:
        np.sqrt(norms, norms)
    return norms

Example 15

Project: scot
Source File: connectivity.py
View license
    @memoize
    def G(self):
        """Inverse cross-spectral density.

        .. math:: \mathbf{G}(f) = \mathbf{A}(f) \mathbf{C}^{-1} \mathbf{A}'(f)
        """
        if self.c is None:
            raise RuntimeError('Inverse cross spectral density requires '
                               'invertible noise covariance matrix c.')
        A = self.A()
        # TODO: can we do that more efficiently?
        G = np.einsum('ji..., jk... ->ik...', A.conj(), self.Cinv())
        G = np.einsum('ij..., jk... ->ik...', G, A)
        return G

Example 16

Project: scot
Source File: connectivity.py
View license
    @memoize
    def COH(self):
        """Coherence.

        .. math:: \mathrm{COH}_{ij}(f) = \\frac{S_{ij}(f)}
                                               {\sqrt{S_{ii}(f) S_{jj}(f)}}

        References
        ----------
        P. L. Nunez, R. Srinivasan, A. F. Westdorp, R. S. Wijesinghe,
        D. M. Tucker, R. B. Silverstein, P. J. Cadusch. EEG coherency I:
        statistics, reference electrode, volume conduction, Laplacians,
        cortical imaging, and interpretation at multiple scales. Electroenceph.
        Clin. Neurophysiol. 103(5): 499-515, 1997.
        """
        S = self.S()
        # TODO: can we do that more efficiently?
        return S / np.sqrt(np.einsum('ii..., jj... ->ij...', S, S.conj()))

Example 17

Project: scot
Source File: connectivity.py
View license
    @memoize
    def pCOH(self):
        """Partial coherence.

        .. math:: \mathrm{pCOH}_{ij}(f) = \\frac{G_{ij}(f)}
                                                {\sqrt{G_{ii}(f) G_{jj}(f)}}

        References
        ----------
        P. J. Franaszczuk, K. J. Blinowska, M. Kowalczyk. The application of
        parametric multichannel spectral estimates in the study of electrical
        brain activity. Biol. Cybernetics 51(4): 239-247, 1985.
        """
        G = self.G()
        # TODO: can we do that more efficiently?
        return G / np.sqrt(np.einsum('ii..., jj... ->ij...', G, G))

Example 18

Project: scot
Source File: connectivity.py
View license
    @memoize
    def PDCF(self):
        """Partial directed coherence factor.

        .. math:: \mathrm{PDCF}_{ij}(f) =
        \\frac{A_{ij}(f)}{\sqrt{A_{:j}'(f) \mathbf{C}^{-1} A_{:j}(f)}}

        References
        ----------
        L. A. Baccalá, K. Sameshima. Partial directed coherence: a new concept
        in neural structure determination. Biol. Cybernetics 84(6): 463-474,
        2001.
        """
        A = self.A()
        # TODO: can we do that more efficiently?
        return np.abs(A / np.sqrt(np.einsum('aj..., ab..., bj... ->j...',
                                            A.conj(), self.Cinv(), A)))

Example 19

Project: scot
Source File: connectivity.py
View license
    @memoize
    def GDTF(self):
        """Generalized directed transfer function.

        .. math:: \mathrm{GPDC}_{ij}(f) = \\frac{\sigma_j |H_{ij}(f)|}
            {\sqrt{H_{i:}(f) \mathrm{diag}(\mathbf{C}) H_{i:}'(f)}}

        References
        ----------
        L. Faes, S. Erla, G. Nollo. Measuring connectivity in linear
        multivariate processes: definitions, interpretation, and practical
        analysis. Comput. Math. Meth. Med. 2012: 140513, 2012.
        """
        H = self.H()
        tmp = H / np.sqrt(np.einsum('ia..., aa..., ia..., j... ->ij...',
                                    H.conj(), self.c, H,
                                    1 / self.c.diagonal()))
        return np.abs(tmp)

Example 20

Project: sfepy
Source File: shell10x.py
View license
def create_strain_matrix(bfgm, dxidx, dsg):
    """
    Create the strain operator matrix.
    """
    tg = create_strain_transform(dxidx)

    mtx_b = nm.concatenate((bfgm[..., 0:1, 0, :],
                            bfgm[..., 1:2, 1, :],
                            bfgm[..., 2:3, 2, :],
                            bfgm[..., 0:1, 1, :] + bfgm[..., 1:2, 0, :],
                            nm.einsum('ciaj,cijk->ciak', tg[..., 4:5, 4:], dsg),
                            nm.einsum('ciaj,cijk->ciak', tg[..., 5:6, 4:], dsg)),
                           2)

    return mtx_b

Example 21

Project: sfepy
Source File: shell10x.py
View license
def rotate_elastic_tensor(mtx_d, bfu, ebs):
    """
    Rotate the elastic tensor into the local coordinate system of each cell.
    The local coordinate system results from interpolation of `ebs` with the
    bilinear basis.
    """
    # ! mtx_ts is not orthonormal for non-planar cells!
    mtx_ts = nm.einsum('qi,cijk->cqjk', bfu[:, 0, :], ebs)
    # print nm.einsum('cqji,cqjk->cqik', mtx_ts, mtx_ts)

    # ? double entries here?
    tt = create_strain_transform(mtx_ts)

    mtx_dr = ddot(ddot(tt, mtx_d, 'ATB'), tt, 'AB')
    return mtx_dr

Example 22

Project: GPy
Source File: basis_funcs.py
View license
    def update_gradients_full(self, dL_dK, X, X2=None):
        if self.ARD:
            phi1 = self.phi(X)
            if X2 is None or X is X2:
                self.variance.gradient = np.einsum('ij,iq,jq->q', dL_dK, phi1, phi1)
            else:
                phi2 = self.phi(X2)
                self.variance.gradient = np.einsum('ij,iq,jq->q', dL_dK, phi1, phi2)
        else:
            self.variance.gradient = np.einsum('ij,ij', dL_dK, self._K(X, X2)) * self.beta

Example 23

Project: GPy
Source File: basis_funcs.py
View license
    def update_gradients_diag(self, dL_dKdiag, X):
        if self.ARD:
            phi1 = self.phi(X)
            self.variance.gradient = np.einsum('i,iq,iq->q', dL_dKdiag, phi1, phi1)
        else:
            self.variance.gradient = np.einsum('i,i', dL_dKdiag, self.Kdiag(X)) * self.beta

Example 24

Project: GPy
Source File: splitKern.py
View license
    def update_gradients_full(self, dL_dK, X, X2=None):
        if X2 is None:
            X2 = X
                        
        k1 = self.kern.K(X,self.Xp)
        k2 = self.kern.K(self.Xp,X2)
        k3 = self.kern.K(self.Xp,self.Xp)
        dL_dk1 = np.einsum('ij,j->i',dL_dK,k2[0])/k3[0,0]
        dL_dk2 = np.einsum('ij,i->j',dL_dK,k1[:,0])/k3[0,0]
        dL_dk3 = np.einsum('ij,ij->',dL_dK,-np.dot(k1,k2)/(k3[0,0]*k3[0,0]))

        self.kern.update_gradients_full(dL_dk1[:,None],X,self.Xp)
        grad = self.kern.gradient.copy()
        self.kern.update_gradients_full(dL_dk2[None,:],self.Xp,X2)
        grad += self.kern.gradient.copy()
        self.kern.update_gradients_full(np.array([[dL_dk3]]),self.Xp,self.Xp)
        grad += self.kern.gradient.copy()
        
        self.kern.gradient = grad

Example 25

Project: GPy
Source File: trunclinear.py
View license
    @Cache_this(limit=3)
    def _product(self, X, X2=None):
        if X2 is None:
            X2 = X
        XX = np.einsum('nq,mq->nmq',X-self.delta,X2-self.delta)
        XX[XX<0] = 0
        return XX

Example 26

Project: GPy
Source File: trunclinear.py
View license
    def update_gradients_full(self, dL_dK, X, X2=None):
        dK_dvar = self._product(X, X2)
        if X2 is None:
            X2=X
        dK_ddelta = self.variances*(2*self.delta-X[:,None,:]-X2[None,:,:])*(dK_dvar>0)
        if self.ARD:
            self.variances.gradient[:] = np.einsum('nmq,nm->q',dK_dvar,dL_dK)
            self.delta.gradient[:] = np.einsum('nmq,nm->q',dK_ddelta,dL_dK)
        else:
            self.variances.gradient[:] = np.einsum('nmq,nm->',dK_dvar,dL_dK)
            self.delta.gradient[:] = np.einsum('nmq,nm->',dK_ddelta,dL_dK)

Example 27

Project: GPy
Source File: trunclinear.py
View license
    def gradients_X(self, dL_dK, X, X2=None):
        XX = self._product(X, X2)
        if X2 is None:
            Xp = (self.variances*(X-self.delta))*(XX>0)
        else:
            Xp = (self.variances*(X2-self.delta))*(XX>0)
        if X2 is None:
            return np.einsum('nmq,nm->nq',Xp,dL_dK)+np.einsum('mnq,nm->mq',Xp,dL_dK)
        else:
            return np.einsum('nmq,nm->nq',Xp,dL_dK)

Example 28

Project: GPy
Source File: trunclinear.py
View license
    def update_gradients_full(self, dL_dK, X, X2=None):
        dK_dvar = self._product(X, X2)
        if self.ARD:
            self.variances.gradient[:] = np.einsum('nmq,nm->q',dK_dvar,dL_dK)
        else:
            self.variances.gradient[:] = np.einsum('nmq,nm->',dK_dvar,dL_dK)

Example 29

Project: GPy
Source File: trunclinear.py
View license
    def gradients_X(self, dL_dK, X, X2=None):
        XX = self._product(X, X2)
        if X2 is None:
            Xp = (self.variances*(X-self.delta))*(XX>0)
        else:
            Xp = (self.variances*(X2-self.delta))*(XX>0)
        if X2 is None:
            return np.einsum('nmq,nm->nq',Xp,dL_dK)+np.einsum('mnq,nm->mq',Xp,dL_dK)
        else:
            return np.einsum('nmq,nm->nq',Xp,dL_dK)

Example 30

Project: GPy
Source File: pca.py
View license
    def _dual_eig(self, X):
        dual_eigvals, dual_eigvects = numpy.linalg.eigh(numpy.einsum('ij,kj->ik',X,X))
        relevant_dimensions = numpy.argsort(numpy.abs(dual_eigvals))[-X.shape[1]:]
        eigvals = dual_eigvals[relevant_dimensions]
        eigvects = dual_eigvects[:, relevant_dimensions]
        eigvects = (1. / numpy.sqrt(X.shape[0] * numpy.abs(eigvals))) * X.T.dot(eigvects)
        eigvects /= numpy.sqrt(numpy.diag(eigvects.T.dot(eigvects)))
        return eigvals, eigvects

Example 31

Project: python-skyfield
Source File: positionlib.py
View license
def ITRF_to_GCRS(t, rITRF):  # todo: velocity

    # Todo: wobble

    spin = rot_z(t.gast * tau / 24.0)
    position = einsum('ij...,j...->i...', spin, array(rITRF))
    return einsum('ij...,j...->i...', t.MT, position)

Example 32

Project: python-skyfield
Source File: relativity.py
View license
def light_time_difference(position, observer_position):
    """Returns the difference in light-time, for a star,
      between the barycenter of the solar system and the observer (or
      the geocenter).

    """
    # From 'pos1', form unit vector 'u1' in direction of star or light
    # source.

    dis = length_of(position)
    u1 = position / dis

    # Light-time returned is the projection of vector 'pos_obs' onto the
    # unit vector 'u1' (formed from 'pos1'), divided by the speed of light.

    diflt = einsum('a...,a...', u1, observer_position) / C_AUDAY
    return diflt

Example 33

Project: python-skyfield
Source File: toposlib.py
View license
    def icrf_vector_at(self, t):
        """Return the GCRS position, velocity of this Topos at `t`."""
        topos = self.target
        pos, vel = terra(topos.latitude.radians, topos.longitude.radians,
                         topos.elevation.au, t.gast)
        pos = einsum('ij...,j...->i...', t.MT, pos)
        vel = einsum('ij...,j...->i...', t.MT, vel)
        if topos.x:
            R = rot_y(topos.x * ASEC2RAD)
            pos = einsum('ij...,j...->i...', R, pos)
        if topos.y:
            R = rot_x(topos.y * ASEC2RAD)
            pos = einsum('ij...,j...->i...', R, pos)
        # TODO: also rotate velocity
        return pos, vel

Example 34

Project: bayespy
Source File: concat_gaussian.py
View license
    def _compute_message_to_parent(self, i, m, *u_nodes):
        r = self.slices

        # Pick the proper parts from the message array
        m0 = m[0][...,r[i]:r[i+1]]
        m1 = m[1][...,r[i]:r[i+1],r[i]:r[i+1]]

        # Handle cross-covariance terms by using the mean of the covariate node
        for (j, u) in enumerate(u_nodes):
            if j != i:
                m0 = m0 + 2 * np.einsum(
                    '...ij,...j->...i',
                    m[1][...,r[i]:r[i+1],r[j]:r[j+1]],
                    u[0]
                )

        return [m0, m1]

Example 35

Project: bayespy
Source File: gamma.py
View license
    @staticmethod
    def _compute_message_to_parent(index, m_children, *u_parents):

        # Take the diagonal
        m0 = np.einsum('...ii->...i', m_children[0])
        m1 = np.reshape(m_children[1], np.shape(m_children[1]) + (1,))

        return [m0, m1]

Example 36

Project: bayespy
Source File: linalg.py
View license
def m_dot(A,b):
    raise DeprecationWarning()
    # Compute matrix-vector product over the last two axes of A and
    # the last axes of b.  Other axes are broadcasted. If A has shape
    # (..., M, N) and b has shape (..., N), then the result has shape
    # (..., M)
    
    #b = reshape(b, shape(b)[:-1] + (1,) + shape(b)[-1:])
    #return np.dot(A, b)
    return np.einsum('...ik,...k->...i', A, b)

Example 37

Project: bayespy
Source File: misc.py
View license
def m_dot(A,b):
    # Compute matrix-vector product over the last two axes of A and
    # the last axes of b.  Other axes are broadcasted. If A has shape
    # (..., M, N) and b has shape (..., N), then the result has shape
    # (..., M)
    
    #b = reshape(b, shape(b)[:-1] + (1,) + shape(b)[-1:])
    #return np.dot(A, b)
    return np.einsum('...ik,...k->...i', A, b)

Example 38

Project: pyscf
Source File: ccsd.py
View license
def energy(cc, t1, t2, eris):
    nocc, nvir = t1.shape
    fock = eris.fock
    e = numpy.einsum('ia,ia', fock[:nocc,nocc:], t1) * 2
    tau = numpy.empty((1,nocc,nvir,nvir))
    for p0 in range(nocc):
        p1 = p0 + 1
        make_tau(t2[p0:p1], t1[p0:p1], t1, 1, out=tau)
        theta = tau*2 - tau.transpose(0,1,3,2)
        e += numpy.einsum('ijab,ijab', theta,
                          eris.ovov[p0:p1].transpose(0,2,1,3))
    return e

Example 39

Project: pyscf
Source File: ccsd_incore.py
View license
def energy(cc, t1, t2, eris):
    nocc, nvir = t1.shape
    fock = eris.fock
    e = numpy.einsum('ia,ia', fock[:nocc,nocc:], t1) * 2
    tau = numpy.empty((1,nocc,nvir,nvir))
    for p0 in range(nocc):
        p1 = p0 + 1
        make_tau(t2[p0:p1], t1[p0:p1], t1, 1, out=tau)
        theta = tau*2 - tau.transpose(0,1,3,2)
        e += numpy.einsum('ijab,ijab', theta,
                          eris.ovov[p0:p1].transpose(0,2,1,3))
    return e

Example 40

Project: pyscf
Source File: uks.py
View license
def energy_elec(ks, dm, h1e=None, vhf=None):
    if h1e is None:
        h1e = ks.get_hcore()
    e1 = numpy.einsum('ij,ji', h1e, dm[0]) + numpy.einsum('ij,ji', h1e, dm[1])
    tot_e = e1.real + ks._ecoul + ks._exc
    logger.debug(ks, 'Ecoul = %s  Exc = %s', ks._ecoul, ks._exc)
    return tot_e, ks._ecoul+ks._exc

Example 41

Project: pyscf
Source File: spin_op.py
View license
def spin_square0(fcivec, norb, nelec):
    '''Spin square for RHF-FCI CI wfn only (obtained from spin-degenerated
    Hamiltonian)'''
    ci1 = contract_ss(fcivec, norb, nelec)
    ss = numpy.einsum('ij,ij->', fcivec.reshape(ci1.shape), ci1)
    s = numpy.sqrt(ss+.25) - .5
    multip = s*2+1
    return ss, multip

Example 42

Project: pyscf
Source File: fciqmc.py
View license
def one_from_two_pdm(two_pdm, nelec):
    '''Return a 1-rdm, given a 2-rdm to contract.

    Args:
        two_pdm : ndarray
            A (spin-free) 2-particle reduced density matrix.
        nelec: int
            The number of electrons contributing to the RDMs.

    Returns:
        one_pdm : ndarray
            The (spin-free) 1-particle reduced density matrix.
    '''

    # Last two indices refer to middle two second quantized operators in the 2RDM
    one_pdm = numpy.einsum('ikjj->ik', two_pdm)
    one_pdm /= (numpy.sum(nelec)-1)
    return one_pdm

Example 43

Project: pyscf
Source File: boys.py
View license
def dipole_integral(mol, mo_coeff):
    # The gauge origin has no effects for maximization |<r>|^2
    # Set to charge center for physical significance of <r>
    charge_center = numpy.einsum('z,zx->x', mol.atom_charges(), mol.atom_coords())
    mol.set_common_orig(charge_center)
    dip = numpy.asarray([reduce(numpy.dot, (mo_coeff.T, x, mo_coeff))
                         for x in mol.intor_symmetric('cint1e_r_sph', comp=3)])
    return dip

Example 44

Project: pyscf
Source File: test_mdf_jk.py
View license
    def test_hcore(self):
        mf = pscf.RHF(cell)
        odf = mdf.MDF(cell)
        odf.auxbasis = 'weigend'
        odf.gs = (10,)*3
        dm = mf.get_init_guess()
        self.assertAlmostEqual(numpy.einsum('ij,ji->', odf.get_nuc(), dm),
                               -150.61096872500249, 9)

Example 45

Project: pyscf
Source File: kuks.py
View license
    def energy_elec(self, dm_kpts=None, h1e_kpts=None, vhf=None):
        if h1e_kpts is None: h1e_kpts = self.get_hcore(self.cell, self.kpts)
        if dm_kpts is None: dm_kpts = self.make_rdm1()

        nkpts = len(h1e_kpts)
        e1 = 1./nkpts *(np.einsum('kij,kji', h1e_kpts, dm_kpts[0]) +
                        np.einsum('kij,kji', h1e_kpts, dm_kpts[1])).real

        tot_e = e1 + self._ecoul + self._exc
        logger.debug(self, 'E1 = %s  Ecoul = %s  Exc = %s', e1, self._ecoul, self._exc)
        return tot_e, self._ecoul + self._exc

Example 46

Project: pyscf
Source File: khf.py
View license
def energy_elec(mf, dm_kpts=None, h1e_kpts=None, vhf_kpts=None):
    '''Following pyscf.scf.hf.energy_elec()
    '''
    if dm_kpts is None: dm_kpts = mf.make_rdm1()
    if h1e_kpts is None: h1e_kpts = mf.get_hcore()
    if vhf_kpts is None: vhf_kpts = mf.get_veff(mf.cell, dm_kpts)

    nkpts = len(dm_kpts)
    e1 = 1./nkpts * np.einsum('kij,kji', dm_kpts, h1e_kpts)
    e_coul = 1./nkpts * np.einsum('kij,kji', dm_kpts, vhf_kpts) * 0.5
    if abs(e_coul.imag > 1.e-7):
        raise RuntimeError("Coulomb energy has imaginary part, "
                           "something is wrong!", e_coul.imag)
    e1 = e1.real
    e_coul = e_coul.real
    logger.debug(mf, 'E_coul = %.15g', e_coul)
    return e1+e_coul, e_coul

Example 47

Project: pyscf
Source File: kuhf.py
View license
def energy_elec(mf, dm_kpts=None, h1e_kpts=None, vhf_kpts=None):
    '''Following pyscf.scf.hf.energy_elec()
    '''
    if dm_kpts is None: dm_kpts = mf.make_rdm1()
    if h1e_kpts is None: h1e_kpts = mf.get_hcore()
    if vhf_kpts is None: vhf_kpts = mf.get_veff(mf.cell, dm_kpts)

    nkpts = len(h1e_kpts)
    e1 = 1./nkpts * np.einsum('kij,kji', dm_kpts[0], h1e_kpts)
    e1+= 1./nkpts * np.einsum('kij,kji', dm_kpts[1], h1e_kpts)
    e_coul = 1./nkpts * np.einsum('kij,kji', dm_kpts[0], vhf_kpts[0]) * 0.5
    e_coul+= 1./nkpts * np.einsum('kij,kji', dm_kpts[1], vhf_kpts[1]) * 0.5
    if abs(e_coul.imag > 1.e-10):
        raise RuntimeError("Coulomb energy has imaginary part, "
                           "something is wrong!", e_coul.imag)
    e1 = e1.real
    e_coul = e_coul.real
    logger.debug(mf, 'E_coul = %.15g', e_coul)
    return e1+e_coul, e_coul

Example 48

Project: pyscf
Source File: geom.py
View license
def closest_axes(axes, ref):
    xcomp, ycomp, zcomp = numpy.einsum('ix,jx->ji', axes, ref)
    z_id = numpy.argmax(abs(zcomp))
    xcomp[z_id] = ycomp[z_id] = 0       # remove z
    x_id = numpy.argmax(abs(xcomp))
    ycomp[x_id] = 0                     # remove x
    y_id = numpy.argmax(abs(ycomp))
    return x_id, y_id, z_id

Example 49

Project: PyFR
Source File: elements.py
View license
    def _gen_pnorm_fpts(self):
        smats = self.smat_at_np('fpts').transpose(1, 3, 0, 2)

        # We need to compute |J|*[(J^{-1})^{T}.N] where J is the
        # Jacobian and N is the normal for each fpt.  Using
        # J^{-1} = S/|J| where S are the smats, we have S^{T}.N.
        pnorm_fpts = np.einsum('ijlk,il->ijk', smats, self.basis.norm_fpts)

        # Compute the magnitudes of these flux point normals
        mag_pnorm_fpts = np.einsum('...i,...i', pnorm_fpts, pnorm_fpts)
        mag_pnorm_fpts = np.sqrt(mag_pnorm_fpts)

        # Check that none of these magnitudes are zero
        if np.any(mag_pnorm_fpts < 1e-10):
            raise RuntimeError('Zero face normals detected')

        # Normalize the physical normals at the flux points
        self._norm_pnorm_fpts = pnorm_fpts / mag_pnorm_fpts[..., None]
        self._mag_pnorm_fpts = mag_pnorm_fpts

Example 50

Project: DESMAN
Source File: Eta_Sampler.py
View license
    def logLikelihoodGene(self,variants,cTau,cGamma,eta,cEpsilon):
        """Computes data log likelihood given parameter states"""
        logLL = 0.0
        V = variants.shape[0]
        if eta.sum() > 0:
            gammaR = self.maskGamma(cGamma,eta)
        
            probVS = np.einsum('ijk,lj,km->ilm',cTau,gammaR,cEpsilon)
        
            #loop each variant

            for v in range(V):
            #loop each site
                for s in range(self.S):                    
                    logLL += log_multinomial_pdf(variants[v,s,:], probVS[v,s,:])
        else:
            logLL = float(V)*ETA_PENALTY
        return logLL