Python numpy.einsum() 使用实例

The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

def total_length_selected(ed='empty', coords='empty', ob='empty'):
    '''Returns the total length of all edge segments'''
    if ob == 'empty':
        ob = bpy.context.object
    if coords == 'empty':    
        coords = get_coords(ob)
    if ed == 'empty':    
        ed = get_edge_idx(ob)
    edc = coords[ed]
    e1 = edc[:, 0]
    e2 = edc[:, 1]
    ee1 = e1 - e2
    sel = get_selected_edges(ob)    
    ee = ee1[sel]    
    leng = np.einsum('ij,ij->i', ee, ee)
    return np.sum(np.sqrt(leng)) 

Example 2

def transitions_old(width, height, configs=None, one_per_state=False):
    digit = width * height
    if configs is None:
        configs = generate_configs(digit)
    if one_per_state:
        def pickone(thing):
            index = np.random.randint(0,len(thing))
            return thing[index]
        transitions = np.array([
            generate(
                [c1,pickone(successors(c1,width,height))],width,height)
            for c1 in configs ])
    else:
        transitions = np.array([ generate([c1,c2],width,height)
                                 for c1 in configs for c2 in successors(c1,width,height) ])
    return np.einsum('ab...->ba...',transitions) 

Example 3

def puzzle_plot(p):
    p.setup()
    def name(template):
        return template.format(p.__name__)
    from itertools import islice
    configs = list(islice(p.generate_configs(9), 1000)) # be careful, islice is not immutable!!!
    import numpy.random as random
    random.shuffle(configs)
    configs = configs[:10]
    puzzles = p.generate(configs, 3, 3)
    print(puzzles.shape, "mean", puzzles.mean(), "stdev", np.std(puzzles))
    plot_image(puzzles[-1], name("{}.png"))
    plot_image(np.clip(puzzles[-1]+np.random.normal(0,0.1,puzzles[-1].shape),0,1),name("{}+noise.png"))
    plot_image(np.round(np.clip(puzzles[-1]+np.random.normal(0,0.1,puzzles[-1].shape),0,1)),name("{}+noise+round.png"))
    plot_grid(puzzles, name("{}s.png"))
    _transitions = p.transitions(3,3,configs=configs)
    print(_transitions.shape)
    transitions_for_show = \
        np.einsum('ba...->ab...',_transitions) \
          .reshape((-1,)+_transitions.shape[2:])
    print(transitions_for_show.shape)
    plot_grid(transitions_for_show, name("{}_transitions.png")) 

Example 4

def run(ae,xs):
    zs = ae.encode_binary(xs)
    ys = ae.decode_binary(zs)
    mod_ys = []
    correlations = []
    print(ys.shape)
    print("corrlations:")
    print("bit \ image  {}".format(range(len(xs))))
    for i in range(ae.N):
        mod_zs = np.copy(zs)
        # increase the latent value from 0 to 1 and check the difference
        for j in range(11):
            mod_zs[:,i] = j / 10.0
            mod_ys.append(ae.decode_binary(mod_zs))
        zero_zs,one_zs = np.copy(zs),np.copy(zs)
        zero_zs[:,i] = 0.
        one_zs[:,i] = 1.
        correlation = np.mean(np.square(ae.decode_binary(zero_zs) - ae.decode_binary(one_zs)),
                              axis=(1,2))
        correlations.append(correlation)
        print("{:>5} {}".format(i,correlation))
    plot_grid2(np.einsum("ib...->bi...",np.array(mod_ys)).reshape((-1,)+ys.shape[1:]),
               w=11,path=ae.local("dump_significance.png"))
    return np.einsum("ib->bi",correlations) 

Example 5

def buildFock(self):
        """Routine to build the AO basis Fock matrix"""
        if self.direct:
            if self.incFockRst: # restart incremental fock build?
                self.G = formPT(self.P,np.zeros_like(self.P),self.bfs,
                                self.nbasis,self.screen,self.scrTol)
                self.G = 0.5*(self.G + self.G.T) 
                self.F = self.Core.astype('complex') + self.G
            else:
                self.G = formPT(self.P,self.P_old,self.bfs,self.nbasis,
                                self.screen,self.scrTol)
                self.G = 0.5*(self.G + self.G.T) 
                self.F = self.F_old + self.G

        else:
            self.J = np.einsum('pqrs,sr->pq', self.TwoE.astype('complex'),self.P)
            self.K = np.einsum('psqr,sr->pq', self.TwoE.astype('complex'),self.P)
            self.G = 2.*self.J - self.K
            self.F = self.Core.astype('complex') + self.G 

Example 6

def pointsInRegion(regNum, vor, p, overlap=0.0):
    """
    returns the subset of points p that are inside the regNum region of the voronoi object
    vor. The boundaries of the region are extended by an amount given by 'overlap'.
    """
    reg = vor.regions[vor.point_region[regNum]]  # region associated with the point
    if -1 in reg:
        raise Exception('Open region associated with generator')
    nVerts = len(reg)  # number of verticies in the region
    p0 = vor.points[regNum]

    for i in range(len(reg)):
        vert1, vert2 = vor.vertices[reg[i]], vor.vertices[reg[(i + 1) % len(reg)]]
        dr = vert1 - vert2  # edge
        dr = dr / numpy.linalg.norm(dr)  # normalize
        dn = numpy.array([dr[1], -dr[0]])  # normal to edge
        dn = dn if numpy.dot(dn, vert2 - p0[:2]) > 0 else -dn  # orient so that the normal is outwards
        d1 = numpy.einsum('i,ji', dn, vert2 + dn * overlap - p[:, :2])
        p = p[d1 * numpy.dot(dn, vert2 - p0[:2]) > 0]

    return p 

Example 7

def test_einsum_misc(self):
        # This call used to crash because of a bug in
        # PyArray_AssignZero
        a = np.ones((1, 2))
        b = np.ones((2, 2, 1))
        assert_equal(np.einsum('ij...,j...->i...', a, b), [[[2], [2]]])

        # The iterator had an issue with buffering this reduction
        a = np.ones((5, 12, 4, 2, 3), np.int64)
        b = np.ones((5, 12, 11), np.int64)
        assert_equal(np.einsum('ijklm,ijn,ijn->', a, b, b),
                        np.einsum('ijklm,ijn->', a, b))

        # Issue #2027, was a problem in the contiguous 3-argument
        # inner loop implementation
        a = np.arange(1, 3)
        b = np.arange(1, 5).reshape(2, 2)
        c = np.arange(1, 9).reshape(4, 2)
        assert_equal(np.einsum('x,yx,zx->xzy', a, b, c),
                    [[[1,  3], [3,  9], [5, 15], [7, 21]],
                    [[8, 16], [16, 32], [24, 48], [32, 64]]]) 

Example 8

def test_einsum_all_contig_non_contig_output(self):
        # Issue gh-5907, tests that the all contiguous special case
        # actually checks the contiguity of the output
        x = np.ones((5, 5))
        out = np.ones(10)[::2]
        correct_base = np.ones(10)
        correct_base[::2] = 5
        # Always worked (inner iteration is done with 0-stride):
        np.einsum('mi,mi,mi->m', x, x, x, out=out)
        assert_array_equal(out.base, correct_base)
        # Example 1:
        out = np.ones(10)[::2]
        np.einsum('im,im,im->m', x, x, x, out=out)
        assert_array_equal(out.base, correct_base)
        # Example 2, buffering causes x to be contiguous but
        # special cases do not catch the operation before:
        out = np.ones((2, 2, 2))[..., 0]
        correct_base = np.ones((2, 2, 2))
        correct_base[..., 0] = 2
        x = np.ones((2, 2), np.float32)
        np.einsum('ij,jk->ik', x, x, out=out)
        assert_array_equal(out.base, correct_base) 

Example 9

def dihedral_transform_batch(x):

    g = np.random.randint(low=0, high=8, size=x.shape[0])

    h, w = x.shape[-2:]
    hh = (h - 1) / 2.
    hw = (w - 1) / 2.

    I, J = np.meshgrid(np.linspace(-hh, hh, x.shape[-2]), np.linspace(-hw, hw, x.shape[-1]))
    C = np.r_[[I, J]]
    D4C = np.einsum('...ij,jkl->...ikl', D4, C)
    D4C[:, 0] += hh
    D4C[:, 1] += hw
    D4C = D4C.astype(int)

    x_out = np.empty_like(x)
    for i in range(x.shape[0]):
        I, J = D4C[g[i]]
        x_out[i, :] = x[i][:, J, I]

    return x_out 

Example 10

def get_vol(simplex):
    # Compute the volume via the Cayley-Menger determinant
    # <http://mathworld.wolfram.com/Cayley-MengerDeterminant.html>. One
    # advantage is that it can compute the volume of the simplex indenpendent
    # of the dimension of the space where it's embedded.

    # compute all edge lengths
    edges = numpy.subtract(simplex[:, None], simplex[None, :])
    ei_dot_ej = numpy.einsum('...k,...k->...', edges, edges)

    j = simplex.shape[0] - 1
    a = numpy.empty((j+2, j+2) + ei_dot_ej.shape[2:])
    a[1:, 1:] = ei_dot_ej
    a[0, 1:] = 1.0
    a[1:, 0] = 1.0
    a[0, 0] = 0.0

    a = numpy.moveaxis(a, (0, 1), (-2, -1))
    det = numpy.linalg.det(a)

    vol = numpy.sqrt((-1.0)**(j+1) / 2**j / math.factorial(j)**2 * det)
    return vol 

Example 11

def scalar_product_interval(anchors, indizes_1, indizes_2):
	q = (anchors[1][0]-anchors[0][0])
	
	vector_1 = np.vstack([
		anchors[0][1][indizes_1],     # a_1
		anchors[0][2][indizes_1] * q, # b_1
		anchors[1][1][indizes_1],     # c_1
		anchors[1][2][indizes_1] * q, # d_1
	])
	
	vector_2 = np.vstack([
		anchors[0][1][indizes_2],     # a_2
		anchors[0][2][indizes_2] * q, # b_2
		anchors[1][1][indizes_2],     # c_2
		anchors[1][2][indizes_2] * q, # d_2
	])
	
	return np.einsum(
		vector_1, [0,2],
		sp_matrix, [0,1],
		vector_2, [1,2]
		)*q 

Example 12

def scalar_product_partial(anchors, indizes_1, indizes_2, start):
	q = (anchors[1][0]-anchors[0][0])
	z = (start-anchors[1][0]) / q
	
	vector_1 = np.vstack([
		anchors[0][1][indizes_1],     # a_1
		anchors[0][2][indizes_1] * q, # b_1
		anchors[1][1][indizes_1],     # c_1
		anchors[1][2][indizes_1] * q, # d_1
	])
	
	vector_2 = np.vstack([
		anchors[0][1][indizes_2],     # a_2
		anchors[0][2][indizes_2] * q, # b_2
		anchors[1][1][indizes_2],     # c_2
		anchors[1][2][indizes_2] * q, # d_2
	])
	
	return np.einsum(
		vector_1, [0,2],
		partial_sp_matrix(z), [0,1],
		vector_2, [1,2]
		)*q 

Example 13

def mvl(pha, amp, optimize):
    """Mean Vector Length (Canolty, 2006).

    Parameters
    ----------
    pha : array_like
        Array of phases of shapes (npha, ..., npts)

    amp : array_like
        Array of amplitudes of shapes (namp, ..., npts)

    Returns
    -------
    pac : array_like
        PAC of shape (npha, namp, ...)
    """
    # Number of time points :
    npts = pha.shape[-1]
    return np.abs(np.einsum('i...j, k...j->ik...', amp, np.exp(1j * pha),
                            optimize=optimize)) / npts 

Example 14

def ps(pha, amp, optimize):
    """Phase Synchrony (Penny, 2008; Cohen, 2008).

    Parameters
    ----------
    pha : array_like
        Array of phases of shapes (npha, ..., npts)

    amp : array_like
        Array of amplitudes of shapes (namp, ..., npts)

    Returns
    -------
    pac : array_like
        PAC of shape (npha, namp, ...)
    """
    # Number of time points :
    npts = pha.shape[-1]
    pac = np.einsum('i...j, k...j->ik...', np.exp(-1j * amp), np.exp(1j * pha),
                    optimize=optimize)
    return np.abs(pac) / npts 

Example 15

def half_space(self):
        """Return the half space polytope respresentation of the infinite
        beam."""
        # add half beam width along the normal direction to each of the points
        half = self.normal * self.size / 2
        edges = [Line(self.p1 + half, self.p2 + half),
                 Line(self.p1 - half, self.p2 - half)]

        A = np.ndarray((len(edges), self.dim))
        B = np.ndarray(len(edges))

        for i in range(0, 2):
            A[i, :], B[i] = edges[i].standard

            # test for positive or negative side of line
            if np.einsum('i, i', self.p1._x, A[i, :]) > B[i]:
                A[i, :] = -A[i, :]
                B[i] = -B[i]

        p = pt.Polytope(A, B)
        return p 

Example 16

def forward_prop_random_thru_post_mm(self, model, mx, vx, mu, Su):
        Kuu_noiseless = compute_kernel(
            2 * model.ls, 2 * model.sf, model.zu, model.zu)
        Kuu = Kuu_noiseless + np.diag(jitter * np.ones((self.M, )))
        # TODO: remove inv
        Kuuinv = np.linalg.inv(Kuu)
        A = np.dot(Kuuinv, mu)
        Smm = Su + np.outer(mu, mu)
        B_sto = np.dot(Kuuinv, np.dot(Smm, Kuuinv)) - Kuuinv
        psi0 = np.exp(2.0 * model.sf)
        psi1, psi2 = compute_psi_weave(
            2 * model.ls, 2 * model.sf, mx, vx, model.zu)
        mout = np.einsum('nm,md->nd', psi1, A)
        Bpsi2 = np.einsum('ab,nab->n', B_sto, psi2)[:, np.newaxis]
        vout = psi0 + Bpsi2 - mout**2
        return mout, vout 

Example 17

def _forward_prop_deterministic_thru_post(self, x, return_info=False):
        """Propagate deterministic inputs thru posterior
        
        Args:
            x (float): input values, size K x Din
            return_info (bool, optional): Description
        
        Returns:
            float, size K x Dout: output means
            float, size K x Dout: output variances
        """
        psi0 = np.exp(2 * self.sf)
        psi1 = compute_kernel(2 * self.ls, 2 * self.sf, x, self.zu)
        mout = np.einsum('nm,dm->nd', psi1, self.A)
        Bpsi2 = np.einsum('dab,na,nb->nd', self.B_det, psi1, psi1)
        vout = psi0 + Bpsi2
        if return_info:
            return mout, vout, psi1
        else:
            return mout, vout 

Example 18

def _forward_prop_random_thru_post_mm(self, mx, vx, return_info=False):
        """Propagate uncertain inputs thru posterior, using Moment Matching
        
        Args:
            mx (float): input means, size K x Din
            vx (TYPE): input variances, size K x Din
            return_info (bool, optional): Description
        
        Returns:
            float, size K x Dout: output means
            float, size K x Dout: output variances
        """
        psi0 = np.exp(2.0 * self.sf)
        psi1, psi2 = compute_psi_weave(
            2 * self.ls, 2 * self.sf, mx, vx, self.zu)
        mout = np.einsum('nm,dm->nd', psi1, self.A)
        Bpsi2 = np.einsum('dab,nab->nd', self.B_sto, psi2)
        vout = psi0 + Bpsi2 - mout**2
        if return_info:
            return mout, vout, psi1, psi2
        else:
            return mout, vout 

Example 19

def sample(self, x):
        """Summary
        
        Args:
            x (TYPE): Description
        
        Returns:
            TYPE: Description
        """
        Su = self.Su
        mu = self.mu
        Lu = np.linalg.cholesky(Su)
        epsilon = np.random.randn(self.Dout, self.M)
        u_sample = mu + np.einsum('dab,db->da', Lu, epsilon)

        kff = compute_kernel(2 * self.ls, 2 * self.sf, x, x)
        kff += np.diag(JITTER * np.ones(x.shape[0]))
        kfu = compute_kernel(2 * self.ls, 2 * self.sf, x, self.zu)
        qfu = np.dot(kfu, self.Kuuinv)
        mf = np.einsum('nm,dm->nd', qfu, u_sample)
        vf = kff - np.dot(qfu, kfu.T)
        Lf = np.linalg.cholesky(vf)
        epsilon = np.random.randn(x.shape[0], self.Dout)
        f_sample = mf + np.einsum('ab,bd->ad', Lf, epsilon)
        return f_sample 

Example 20

def _forward_prop_deterministic_thru_cav(self, x):
        """Propagate deterministic inputs thru cavity
        
        Args:
            x (float): input values, size K x Din
        
        Returns:
            float, size K x Dout: output means
            float, size K x Dout: output variances
            float, size K x M: cross covariance matrix
        """
        kff = np.exp(2 * self.sf)
        kfu = compute_kernel(2 * self.ls, 2 * self.sf, x, self.zu)
        mout = np.einsum('nm,dm->nd', kfu, self.Ahat)
        Bkfukuf = np.einsum('dab,na,nb->nd', self.Bhat_det, kfu, kfu)
        vout = kff + Bkfukuf
        return mout, vout, kfu 

Example 21

def _forward_prop_random_thru_cav_mm(self, mx, vx):
        """Propagate uncertain inputs thru cavity, using simple Moment Matching
        
        Args:
            mx (float): input means, size K x Din
            vx (TYPE): input variances, size K x Din
        
        Returns:
            output means and variances, and intermediate info for backprop
        """
        psi0 = np.exp(2 * self.sf)
        psi1, psi2 = compute_psi_weave(
            2 * self.ls, 2 * self.sf, mx, vx, self.zu)
        mout = np.einsum('nm,dm->nd', psi1, self.Ahat)
        Bhatpsi2 = np.einsum('dab,nab->nd', self.Bhat_sto, psi2)
        vout = psi0 + Bhatpsi2 - mout**2
        return mout, vout, psi1, psi2 

Example 22

def psi1compDer(dL_dpsi1, _psi1, variance, lengthscale, Z, mu, S):
    # here are the "statistics" for psi1
    # Produced intermediate results: dL_dparams w.r.t. psi1
    # _dL_dvariance     1
    # _dL_dlengthscale  Q
    # _dL_dZ            MxQ
    # _dL_dgamma        NxQ
    # _dL_dmu           NxQ
    # _dL_dS            NxQ

    lengthscale2 = np.square(lengthscale)

    Lpsi1 = dL_dpsi1 * _psi1
    Zmu = Z[None, :, :] - mu[:, None, :]  # NxMxQ
    denom = 1. / (S + lengthscale2)
    Zmu2_denom = np.square(Zmu) * denom[:, None, :]  # NxMxQ
    _dL_dvar = Lpsi1.sum() / variance
    _dL_dmu = np.einsum('nm,nmq,nq->nq', Lpsi1, Zmu, denom)
    _dL_dS = np.einsum('nm,nmq,nq->nq', Lpsi1, (Zmu2_denom - 1.), denom) / 2.
    _dL_dZ = -np.einsum('nm,nmq,nq->mq', Lpsi1, Zmu, denom)
    _dL_dl = np.einsum('nm,nmq,nq->q', Lpsi1, (Zmu2_denom +
                                               (S / lengthscale2)[:, None, :]), denom * lengthscale)

    return _dL_dvar, _dL_dl, _dL_dZ, _dL_dmu, _dL_dS 

Example 23

def kfucompDer(dL_dkfu, kfu, variance, lengthscale, Z, mu, grad_x):
    # here are the "statistics" for psi1
    # Produced intermediate results: dL_dparams w.r.t. psi1
    # _dL_dvariance     1
    # _dL_dlengthscale  Q
    # _dL_dZ            MxQ

    lengthscale2 = np.square(lengthscale)

    Lpsi1 = dL_dkfu * kfu
    Zmu = Z[None, :, :] - mu[:, None, :]  # NxMxQ
    _dL_dvar = Lpsi1.sum() / variance
    _dL_dZ = -np.einsum('nm,nmq->mq', Lpsi1, Zmu / lengthscale2)
    _dL_dl = np.einsum('nm,nmq->q', Lpsi1, np.square(Zmu) / lengthscale**3)
    if grad_x:
        _dL_dx = np.einsum('nm,nmq->nq', Lpsi1, Zmu / lengthscale2)
        return _dL_dvar, _dL_dl, _dL_dZ, _dL_dx
    else:
        return _dL_dvar, _dL_dl, _dL_dZ 

Example 24

def _forward_prop_deterministic_thru_cav(self, n, x, alpha):
        """Summary

        Args:
            n (TYPE): Description
            x (TYPE): Description
            alpha (TYPE): Description

        Returns:
            TYPE: Description
        """
        muhat, Suhat, SuinvMuhat, Suinvhat = self.compute_cavity(n, alpha)
        Kuuinv = self.Kuuinv
        Ahat = np.einsum('ab,ndb->nda', Kuuinv, muhat)
        Bhat = np.einsum(
            'ab,ndbc->ndac',
            Kuuinv, np.einsum('ndab,bc->ndac', Suhat, Kuuinv)) - Kuuinv
        kff = np.exp(2 * self.sf)
        kfu = compute_kernel(2 * self.ls, 2 * self.sf, x, self.zu)
        mout = np.einsum('nm,ndm->nd', kfu, Ahat)
        Bkfukuf = np.einsum('ndab,na,nb->nd', Bhat, kfu, kfu)
        vout = kff + Bkfukuf
        extra_res = [muhat, Suhat, SuinvMuhat, Suinvhat, kfu, Ahat, Bhat]
        return mout, vout, extra_res 

Example 25

def _forward_prop_deterministic_thru_post(self, x):
        """Summary

        Args:
            x (TYPE): Description

        Returns:
            TYPE: Description
        """
        Kuuinv = self.Kuuinv
        A = np.einsum('ab,db->da', Kuuinv, self.mu)
        B = np.einsum(
            'ab,dbc->dac',
            Kuuinv, np.einsum('dab,bc->dac', self.Su, Kuuinv)) - Kuuinv
        kff = np.exp(2 * self.sf)
        kfu = compute_kernel(2 * self.ls, 2 * self.sf, x, self.zu)
        mout = np.einsum('nm,dm->nd', kfu, A)
        Bpsi2 = np.einsum('dab,na,nb->nd', B, kfu, kfu)
        vout = kff + Bpsi2
        return mout, vout

    # TODO 

Example 26

def _forward_prop_random_thru_post_mm(self, mx, vx):
        """Summary

        Args:
            mx (TYPE): Description
            vx (TYPE): Description

        Returns:
            TYPE: Description
        """
        Kuuinv = self.Kuuinv
        A = np.einsum('ab,db->da', Kuuinv, self.mu)
        Smm = self.Su + np.einsum('da,db->dab', self.mu, self.mu)
        B = np.einsum(
            'ab,dbc->dac',
            Kuuinv, np.einsum('dab,bc->dac', Smm, Kuuinv)) - Kuuinv
        psi0 = np.exp(2.0 * self.sf)
        psi1, psi2 = compute_psi_weave(
            2 * self.ls, 2 * self.sf, mx, vx, self.zu)
        mout = np.einsum('nm,dm->nd', psi1, A)
        Bpsi2 = np.einsum('dab,nab->nd', B, psi2)
        vout = psi0 + Bpsi2 - mout**2
        return mout, vout 

Example 27

def sample(self, x):
        """Summary

        Args:
            x (TYPE): Description

        Returns:
            TYPE: Description
        """
        Su = self.Su
        mu = self.mu
        Lu = np.linalg.cholesky(Su)
        epsilon = np.random.randn(self.Dout, self.M)
        u_sample = mu + np.einsum('dab,db->da', Lu, epsilon)

        kff = compute_kernel(2 * self.ls, 2 * self.sf, x, x)
        kff += np.diag(JITTER * np.ones(x.shape[0]))
        kfu = compute_kernel(2 * self.ls, 2 * self.sf, x, self.zu)
        qfu = np.dot(kfu, self.Kuuinv)
        mf = np.einsum('nm,dm->nd', qfu, u_sample)
        vf = kff - np.dot(qfu, kfu.T)
        Lf = np.linalg.cholesky(vf)
        epsilon = np.random.randn(x.shape[0], self.Dout)
        f_sample = mf + np.einsum('ab,bd->ad', Lf, epsilon)
        return f_sample 

Example 28

def compute_cavity(self, n, alpha=1.0):
        """Summary

        Args:
            n (TYPE): Description
            alpha (float, optional): Description

        Returns:
            TYPE: Description
        """
        # compute the leave one out moments
        t1n = self.t1[n, :, :]
        t2n = self.t2[n, :, :, :]
        Suinvhat = self.Suinv - alpha * t2n
        SuinvMuhat = self.SuinvMu - alpha * t1n
        Suhat = np.linalg.inv(Suinvhat)
        muhat = np.einsum('ndab,ndb->nda', Suhat, SuinvMuhat)
        return muhat, Suhat, SuinvMuhat, Suinvhat 

Example 29

def forward_prop_thru_post(self, x):
        """Summary

        Args:
            x (TYPE): Description

        Returns:
            TYPE: Description
        """
        Kuuinv = self.Kuuinv
        A = np.einsum('ab,db->da', Kuuinv, self.mu)
        B = np.einsum(
            'ab,dbc->dac',
            Kuuinv, np.einsum('dab,bc->dac', self.Su, Kuuinv)) - Kuuinv
        kff = np.exp(2 * self.sf)
        kfu = compute_kernel(2 * self.ls, 2 * self.sf, x, self.zu)
        mout = np.einsum('nm,dm->nd', kfu, A)
        Bpsi2 = np.einsum('dab,na,nb->nd', B, kfu, kfu)
        vout = kff + Bpsi2
        return mout, vout 

Example 30

def update_posterior(self, x_train=None, new_hypers=False):
        """Summary

        Returns:
            TYPE: Description
        """
        # compute the posterior approximation
        if new_hypers and x_train is not None:
            Kfu = compute_kernel(2*self.ls, 2*self.sf, x_train, self.zu)
            KuuinvKuf = np.dot(self.Kuuinv, Kfu.T)
            self.Kfu = Kfu
            self.KuuinvKuf = KuuinvKuf
            self.Kff_diag = compute_kernel_diag(2*self.ls, 2*self.sf, x_train)

        KuuinvKuf_div_var = np.einsum('an,nd->dan', self.KuuinvKuf, 1.0 / self.variances)
        T2u = np.einsum('dan,bn->dab', KuuinvKuf_div_var, self.KuuinvKuf)
        T1u = np.einsum('bn,nd->db', self.KuuinvKuf, self.means / self.variances)
        Vinv = self.Kuuinv + T2u
        self.Suinv = Vinv
        self.Su = np.linalg.inv(Vinv)
        self.mu = np.einsum('dab,db->da', self.Su, T1u)
        self.gamma = np.einsum('ab,db->da', self.Kuuinv, self.mu)
        self.beta = self.Kuuinv - np.einsum('ab,dbc->dac', 
            self.Kuuinv,
            np.einsum('dab,bc->dac', self.Su, self.Kuuinv)) 

Example 31

def compute_cavity(self, idxs, alpha):
        # deletion
        p_i = self.KuuinvKuf[:, idxs].T[:, np.newaxis, :]
        k_i = self.Kfu[idxs, :]
        k_ii = self.Kff_diag[idxs][:, np.newaxis]
        gamma = self.gamma
        beta = self.beta
        h_si = p_i - np.einsum('dab,nb->nda', beta, k_i)
        variance_i = self.variances[idxs, :]
        mean_i = self.means[idxs, :]
        dlogZd_dmi2 = 1.0 / (variance_i/alpha - 
            np.sum(k_i[:, np.newaxis, :] * h_si, axis=2))
        dlogZd_dmi = -dlogZd_dmi2 * (mean_i - 
            np.sum(k_i[:, np.newaxis, :] * gamma, axis=2))
        hd1 = h_si * dlogZd_dmi[:, :, np.newaxis]
        hd2h = np.einsum('nda,ndb->ndab', h_si, h_si) * dlogZd_dmi2[:, :, np.newaxis, np.newaxis]
        gamma_si = gamma + hd1
        beta_si = beta - hd2h

        # projection
        h = p_i - np.einsum('ndab,nb->nda', beta_si, k_i)
        m_si_i = np.einsum('na,nda->nd', k_i, gamma_si)
        v_si_ii = k_ii - np.einsum('na,ndab,nb->nd', k_i, beta_si, k_i)

        return m_si_i, v_si_ii, [h, beta_si, gamma_si] 

Example 32

def project3Dto2D(self, Li, idxs):
        """
        Project 3D point to 2D
        :param Li: joints in normalized 3D
        :param idxs: frames specified by subset
        :return: 2D points, in normalized 2D coordinates
        """

        if not isinstance(idxs, numpy.ndarray):
            idxs = numpy.asarray([idxs])

        # 3D -> 2D projection also shift by M to cropped window
        Li_glob3D = (numpy.reshape(Li, (len(idxs), self.numJoints, 3))*self.Di_scale[idxs][:, None, None]+self.Di_off3D[idxs][:, None, :]).reshape((len(idxs)*self.numJoints, 3))
        Li_glob3D_hom = numpy.concatenate([Li_glob3D, numpy.ones((len(idxs)*self.numJoints, 1), dtype='float32')], axis=1)
        Li_glob2D_hom = numpy.dot(Li_glob3D_hom, self.cam_proj.T)
        Li_glob2D = (Li_glob2D_hom[:, 0:3] / Li_glob2D_hom[:, 3][:, None]).reshape((len(idxs), self.numJoints, 3))
        Li_img2D_hom = numpy.einsum('ijk,ikl->ijl', Li_glob2D, self.Di_trans2D[idxs])
        Li_img2D = (Li_img2D_hom[:, :, 0:2] / Li_img2D_hom[:, :, 2][:, :, None]).reshape((len(idxs), self.numJoints*2))
        Li_img2Dcrop = (Li_img2D - (self.Di.shape[3]/2.)) / (self.Di.shape[3]/2.)
        return Li_img2Dcrop 

Example 33

def compute_dr_wrt(self, wrt):
            if wrt is not self.v:
                return None

            v = self.v.r.reshape(-1, 3)
            blocks = -np.einsum('ij,ik->ijk', v, v) * (self.ss**(-3./2.)).reshape((-1, 1, 1))
            for i in range(3):
                blocks[:, i, i] += self.s_inv

            if True: # pylint: disable=using-constant-test
                data = blocks.ravel()
                indptr = np.arange(0, (self.v.r.size+1)*3, 3)
                indices = col(np.arange(0, self.v.r.size))
                indices = np.hstack([indices, indices, indices])
                indices = indices.reshape((-1, 3, 3))
                indices = indices.transpose((0, 2, 1)).ravel()
                result = sp.csc_matrix((data, indices, indptr), shape=(self.v.r.size, self.v.r.size))
                return result
            else:
                matvec = lambda x: np.einsum('ijk,ik->ij', blocks, x.reshape((blocks.shape[0], 3))).ravel()
                return sp.linalg.LinearOperator((self.v.r.size, self.v.r.size), matvec=matvec) 

Example 34

def test_einsum_misc(self):
        # This call used to crash because of a bug in
        # PyArray_AssignZero
        a = np.ones((1, 2))
        b = np.ones((2, 2, 1))
        assert_equal(np.einsum('ij...,j...->i...', a, b), [[[2], [2]]])

        # The iterator had an issue with buffering this reduction
        a = np.ones((5, 12, 4, 2, 3), np.int64)
        b = np.ones((5, 12, 11), np.int64)
        assert_equal(np.einsum('ijklm,ijn,ijn->', a, b, b),
                        np.einsum('ijklm,ijn->', a, b))

        # Issue #2027, was a problem in the contiguous 3-argument
        # inner loop implementation
        a = np.arange(1, 3)
        b = np.arange(1, 5).reshape(2, 2)
        c = np.arange(1, 9).reshape(4, 2)
        assert_equal(np.einsum('x,yx,zx->xzy', a, b, c),
                    [[[1,  3], [3,  9], [5, 15], [7, 21]],
                    [[8, 16], [16, 32], [24, 48], [32, 64]]]) 

Example 35

def test_einsum_all_contig_non_contig_output(self):
        # Issue gh-5907, tests that the all contiguous special case
        # actually checks the contiguity of the output
        x = np.ones((5, 5))
        out = np.ones(10)[::2]
        correct_base = np.ones(10)
        correct_base[::2] = 5
        # Always worked (inner iteration is done with 0-stride):
        np.einsum('mi,mi,mi->m', x, x, x, out=out)
        assert_array_equal(out.base, correct_base)
        # Example 1:
        out = np.ones(10)[::2]
        np.einsum('im,im,im->m', x, x, x, out=out)
        assert_array_equal(out.base, correct_base)
        # Example 2, buffering causes x to be contiguous but
        # special cases do not catch the operation before:
        out = np.ones((2, 2, 2))[..., 0]
        correct_base = np.ones((2, 2, 2))
        correct_base[..., 0] = 2
        x = np.ones((2, 2), np.float32)
        np.einsum('ij,jk->ik', x, x, out=out)
        assert_array_equal(out.base, correct_base) 

Example 36

def test_exKxz_pairwise(self):
        covall = np.array([self.Xcov, self.Xcovc])
        for k in self.kernels:
            with self.test_context():
                if isinstance(k, ekernels.Linear):
                    continue
                k.compile()
                exKxz = k.compute_exKxz_pairwise(self.Z, self.Xmu, covall)
                Kxz = k.compute_K(self.Xmu[:-1, :], self.Z)  # NxM
                xKxz = np.einsum('nm,nd->nmd', Kxz, self.Xmu[1:, :])
                self.assertTrue(np.allclose(xKxz, exKxz))

#    def test_exKxz(self):
#        for k in self.kernels:
#            with self.test_session():
#                if isinstance(k, ekernels.Linear):
#                    continue
#                k.compile()
#                exKxz = k.compute_exKxz(self.Z, self.Xmu, self.Xcov)
#                Kxz = k.compute_K(self.Xmu, self.Z)  # NxM
#                xKxz = np.einsum('nm,nd->nmd', Kxz, self.Xmu)
#                self.assertTrue(np.allclose(xKxz, exKxz)) 

Example 37

def _accumulate_sufficient_statistics(self, stats, obs, framelogprob,
                                          posteriors, fwdlattice, bwdlattice):
        super(GaussianHMM, self)._accumulate_sufficient_statistics(
            stats, obs, framelogprob, posteriors, fwdlattice, bwdlattice)

        if 'm' in self.params or 'c' in self.params:
            stats['post'] += posteriors.sum(axis=0)
            stats['obs'] += np.dot(posteriors.T, obs)

        if 'c' in self.params:
            if self.covariance_type in ('spherical', 'diag'):
                stats['obs**2'] += np.dot(posteriors.T, obs ** 2)
            elif self.covariance_type in ('tied', 'full'):
                # posteriors: (nt, nc); obs: (nt, nf); obs: (nt, nf)
                # -> (nc, nf, nf)
                stats['obs*obs.T'] += np.einsum(
                    'ij,ik,il->jkl', posteriors, obs, obs) 

Example 38

def __init__(self, matrix, w):
        W = np.sum(w)
        self.w = w
        self.X = matrix
        self.left = None
        self.right = None
        self.mu = np.einsum('ij,i->j', self.X, w)/W
        diff = self.X - np.tile(self.mu, [np.shape(self.X)[0], 1])
        t = np.einsum('ij,i->ij', diff, np.sqrt(w))
        self.cov = (t.T @ t)/W + 1e-5*np.eye(3)
        self.N = self.X.shape[0]
        V, D = np.linalg.eig(self.cov)
        self.lmbda = np.max(np.abs(V))
        self.e = D[np.argmax(np.abs(V))]


# S is measurements vector - dim nxd
# w is weights vector - dim n 

Example 39

def _solve_2D2(self, X, Z):
        """Solves :math:`Z^T N^{-1}X`, where :math:`X`
        and :math:`Z` are 2-d arrays.
        """

        ZNX = np.dot(Z.T / self._nvec, X)
        for slc, jv in zip(self._slices, self._jvec):
            if slc.stop - slc.start > 1:
                Zblock = Z[slc, :]
                Xblock = X[slc, :]
                niblock = 1 / self._nvec[slc]
                beta = 1.0 / (np.einsum('i->', niblock)+1.0/jv)
                zn = np.dot(niblock, Zblock)
                xn = np.dot(niblock, Xblock)
                ZNX -= beta * np.outer(zn.T, xn)
        return ZNX 

Example 40

def ss_framerotate(mjd, planet, x, y, z, dz,
                   offset=None, equatorial=False):
    """
    Rotate planet trajectory given as (n,3) tensor,
    by ecliptic Euler angles x, y, z, and by z rate
    dz. The rate has units of rad/year, and is referred
    to offset 2010/1/1. dates must be given in MJD.
    """
    if equatorial:
        planet = eq2ecl_vec(planet)

    E = euler_vec(z + dz * (mjd - t_offset) / 365.25, y, x,
                  planet.shape[0])

    planet = np.einsum('ijk,ik->ij', E, planet)

    if offset is not None:
        planet = np.array(offset) + planet

    if equatorial:
        planet = ecl2eq_vec(planet)

    return planet 

Example 41

def _log_prod_students_t(self, i, mu, log_prod_var, inv_var, v):
        """
        Return the value of the log of the product of the univariate Student's
        t PDFs at `X[i]`.
        """
        delta = self.X[i, :] - mu
        return (
            self.D * (
                self._cached_gammaln_by_2[v + 1] - self._cached_gammaln_by_2[v]
                - 0.5*self._cached_log_v[v] - 0.5*self._cached_log_pi
                )
            - 0.5*log_prod_var
            - (v + 1.)/2. * (np.log(1. + 1./v * np.square(delta) * inv_var)).sum()
            )


#-----------------------------------------------------------------------------#
#                              UTILITY FUNCTIONS                              #
#-----------------------------------------------------------------------------#

# Below is slightly faster than np.sum, see http://stackoverflow.com/questions/
# 18365073/why-is-numpys-einsum-faster-than-numpys-built-in-functions 

Example 42

def setUp(self):
        self.f = links.Bilinear(
            self.in_shape[0], self.in_shape[1], self.out_size)
        self.f.W.data[...] = _uniform(*self.f.W.data.shape)
        self.f.V1.data[...] = _uniform(*self.f.V1.data.shape)
        self.f.V2.data[...] = _uniform(*self.f.V2.data.shape)
        self.f.b.data[...] = _uniform(*self.f.b.data.shape)
        self.f.zerograds()

        self.W = self.f.W.data.copy()
        self.V1 = self.f.V1.data.copy()
        self.V2 = self.f.V2.data.copy()
        self.b = self.f.b.data.copy()

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

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

Example 43

def inside_triangle(point,triangles):
    v0 = triangles[:,2]-triangles[:,0]
    v1 = triangles[:,1]-triangles[:,0]
    v2 = point-triangles[:,0]

    dot00 = np.einsum('ij,ij->i',v0,v0)
    dot01 = np.einsum('ij,ij->i',v0,v1)
    dot02 = np.einsum('ij,ij->i',v0,v2)
    dot11 = np.einsum('ij,ij->i',v1,v1)
    dot12 = np.einsum('ij,ij->i',v1,v2)
    
    invDenom = 1./(dot00 * dot11-dot01*dot01)
    u = np.float16((dot11 * dot02 - dot01 * dot12)*invDenom)
    v = np.float16((dot00 * dot12 - dot01 * dot02)*invDenom)

    return (u>=0) & (v>=0) & (u+v<=1) 

Example 44

def omgp_model_bound(omgp):
    ''' Calculate the part of the omgp bound which does not depend
    on the response variable.
    '''
    GP_bound = 0.0

    LBs = []
    # Precalculate the bound minus data fit,
    # and LB matrices used for data fit term.
    for i, kern in enumerate(omgp.kern):
        K = kern.K(omgp.X)
        B_inv = np.diag(1. / ((omgp.phi[:, i] + 1e-6) / omgp.variance))
        Bi, LB, LBi, Blogdet = pdinv(K + B_inv)
        LBs.append(LB)

        # Penalty
        GP_bound -= 0.5 * Blogdet

        # Constant
        GP_bound -= 0.5 * omgp.D * np.einsum('j,j->', omgp.phi[:, i], np.log(2 * np.pi * omgp.variance))

    model_bound = GP_bound + omgp.mixing_prop_bound() + omgp.H

    return model_bound, LBs 

Example 45

def _predict_output_derivatives(self, x):
        n = x.shape[0]
        nt = self.nt
        ny = self.training_points[None][0][1].shape[1]
        num = self.num

        dy_dstates = np.empty(n * num['dof'])
        self.rbfc.compute_jac(n, x.flatten(), dy_dstates)
        dy_dstates = dy_dstates.reshape((n, num['dof']))

        dstates_dytl = np.linalg.inv(self.mtx)

        ones = np.ones(self.nt)
        arange = np.arange(self.nt)
        dytl_dyt = csc_matrix((ones, (arange, arange)), shape=(num['dof'], self.nt))

        dy_dyt = (dytl_dyt.T.dot(dstates_dytl.T).dot(dy_dstates.T)).T
        dy_dyt = np.einsum('ij,k->ijk', dy_dyt, np.ones(ny))
        return {None: dy_dyt} 

Example 46

def get_power_spectral_density_matrix(observation, mask=None):
    """
    Calculates the weighted power spectral density matrix.

    This does not yet work with more than one target mask.

    :param observation: Complex observations with shape (bins, sensors, frames)
    :param mask: Masks with shape (bins, frames) or (bins, 1, frames)
    :return: PSD matrix with shape (bins, sensors, sensors)
    """
    bins, sensors, frames = observation.shape

    if mask is None:
        mask = np.ones((bins, frames))
    if mask.ndim == 2:
        mask = mask[:, np.newaxis, :]

    normalization = np.maximum(np.sum(mask, axis=-1, keepdims=True), 1e-6)

    psd = np.einsum('...dt,...et->...de', mask * observation,
                    observation.conj())
    psd /= normalization
    return psd 

Example 47

def get_mvdr_vector(atf_vector, noise_psd_matrix):
    """
    Returns the MVDR beamforming vector.

    :param atf_vector: Acoustic transfer function vector
        with shape (..., bins, sensors)
    :param noise_psd_matrix: Noise PSD matrix
        with shape (bins, sensors, sensors)
    :return: Set of beamforming vectors with shape (..., bins, sensors)
    """

    while atf_vector.ndim > noise_psd_matrix.ndim - 1:
        noise_psd_matrix = np.expand_dims(noise_psd_matrix, axis=0)

    # Make sure matrix is hermitian
    noise_psd_matrix = 0.5 * (
        noise_psd_matrix + np.conj(noise_psd_matrix.swapaxes(-1, -2)))

    numerator = solve(noise_psd_matrix, atf_vector)
    denominator = np.einsum('...d,...d->...', atf_vector.conj(), numerator)
    beamforming_vector = numerator / np.expand_dims(denominator, axis=-1)

    return beamforming_vector 

Example 48

def apply_sdw_mwf(mix, target_psd_matrix, noise_psd_matrix, mu=1, corr=None):
    """
    Apply speech distortion weighted MWF: h = Tpsd * e1 / (Tpsd + mu*Npsd) 
    :param mix: the signal complex FFT
    :param target_psd_matrix (bins, sensors, sensors) 
    :param noise_psd_matrix
    :param mu: the lagrange factor
    :return 
    """
    bins, sensors, frames = mix.shape
    ref_vector = np.zeros((sensors,1), dtype=np.float)    
    if corr is None:
        ref_ch = 0
    else: # choose the channel with highest correlation with the others
        corr=corr.tolist()        
        while len(corr) > sensors:
            corr.remove(np.min(corr))
        ref_ch=np.argmax(corr)
    ref_vector[ref_ch,0]=1 
    
    mwf_vector = solve(target_psd_matrix + mu*noise_psd_matrix, target_psd_matrix[:,:,ref_ch])
    return np.einsum('...a,...at->...t', mwf_vector.conj(), mix) 

Example 49

def test_einsum_misc(self):
        # This call used to crash because of a bug in
        # PyArray_AssignZero
        a = np.ones((1, 2))
        b = np.ones((2, 2, 1))
        assert_equal(np.einsum('ij...,j...->i...', a, b), [[[2], [2]]])

        # The iterator had an issue with buffering this reduction
        a = np.ones((5, 12, 4, 2, 3), np.int64)
        b = np.ones((5, 12, 11), np.int64)
        assert_equal(np.einsum('ijklm,ijn,ijn->', a, b, b),
                        np.einsum('ijklm,ijn->', a, b))

        # Issue #2027, was a problem in the contiguous 3-argument
        # inner loop implementation
        a = np.arange(1, 3)
        b = np.arange(1, 5).reshape(2, 2)
        c = np.arange(1, 9).reshape(4, 2)
        assert_equal(np.einsum('x,yx,zx->xzy', a, b, c),
                    [[[1,  3], [3,  9], [5, 15], [7, 21]],
                    [[8, 16], [16, 32], [24, 48], [32, 64]]]) 

Example 50

def test_against_numpy_einsum(self):
        """ Test against numpy.einsum  """
        a = np.arange(60.).reshape(3,4,5)
        b = np.arange(24.).reshape(4,3,2)
        stream = [a, b]

        from_numpy = np.einsum('ijk,jil->kl', a, b)
        from_stream = last(ieinsum(stream, 'ijk,jil->kl'))

        self.assertSequenceEqual(from_numpy.shape, from_stream.shape)
        self.assertTrue(np.allclose(from_numpy, from_stream)) 
点赞