Python numpy.tensordot函数代码示例

    def backward_cpu(self, inputs, grad_outputs):
        x, W = inputs[:2]
        b = inputs[2] if len(inputs) == 3 else None

        if not type_check.same_types(*inputs):
            if b is not None:
                raise ValueError('numpy and cupy must not be used together\n'
                                 'type(W): {0}, type(x): {1}, type(b): {2}'
                                 .format(type(W), type(x), type(b)))
            else:
                raise ValueError('numpy and cupy must not be used together\n'
                                 'type(W): {0}, type(x): {1}'
                                 .format(type(W), type(x)))

        gy = grad_outputs[0]
        h, w = x.shape[2:]

        gW = numpy.tensordot(
            gy, self.col, ((0, 2, 3), (0, 4, 5))).astype(W.dtype, copy=False)

        if not self.requires_x_grad:
            gx = None
        else:
            gcol = numpy.tensordot(W, gy, (0, 1)).astype(x.dtype, copy=False)
            gcol = numpy.rollaxis(gcol, 3)
            gx = conv.col2im_cpu(gcol, self.sy, self.sx, self.ph, self.pw,
                                 h, w)

        if b is None:
            return gx, gW
        else:
            gb = gy.sum(axis=(0, 2, 3))
            return gx, gW, gb

def innerProductOBC(mpsA,mpsB):
    """ Inner product <A|B> using transfer matrices
        where A and B are MPS representations of }A> and }B>
        with open boundary conditions (OBC).
    """
    # Take adjoint of |A> to get <A|
    A = []
    for a in mpsA:
        A.append(np.conj(a))
    
    B = mpsB
    N = len(A) # Number of qubits
    d = A[1].shape[1] # d = 2 for qubits

    # Construct list of transfer matrices by contracting pairs of
    # tensors from A and B.
    transfer = []
    t = np.tensordot(A[0],B[0],axes=(0,0))
    t = np.reshape(t,A[0].shape[1]*B[0].shape[1])
    transfer.append(t)
    for i in xrange(1,N-1):
        t = np.tensordot(A[i],B[i],axes=(1,1))
        t = np.transpose(t,axes=(0,2,1,3))
        t = np.reshape(t,(A[i].shape[0]*B[i].shape[0],
                          A[i].shape[2]*B[i].shape[2]))
        transfer.append(t)
    t = np.tensordot(A[N-1],B[N-1],axes=(1,1))
    t = np.reshape(t,A[N-1].shape[0]*B[N-1].shape[0])
    transfer.append(t)
    # Contract the transfer matrices.
    prod = transfer[0]
    for i in xrange(1,N-1):
        prod = np.dot(prod,transfer[i])
    prod = np.dot(prod,transfer[N-1])
    return prod

    def r_log_spiral(self,phi):
	"""
	return distance from center for angle phi of logarithmic spiral

	Parameters
	----------
	phi: scalar or np.array with polar angle values

	Returns
	-------
	r(phi) = rx * exp(b * phi) as np.array

	Notes
	-----
	see http://en.wikipedia.org/wiki/Logarithmic_spiral
	"""
	if np.isscalar(phi):
	    phi = np.array([phi])
	ones = np.ones(phi.shape[0])

	# self.rx.shape = 8
	# phi.shape = p
	# then result is given as (8,p)-dim array, each row stands for one rx

	result = np.tensordot(self.rx , np.exp((phi - 3.*pi*ones) / np.tan(pi/2. - self.idisk)),axes = 0)
	result = np.vstack((result, np.tensordot(self.rx , np.exp((phi - pi*ones) / np.tan(pi/2. - self.idisk)),axes = 0) ))
	result = np.vstack((result, np.tensordot(self.rx , np.exp((phi + pi*ones) / np.tan(pi/2. - self.idisk)),axes = 0) ))
	return np.vstack((result, np.tensordot(self.rx , np.exp((phi + 3.*pi*ones) / np.tan(pi/2. - self.idisk)),axes = 0) ))

def tensordot_adjoint_1(A, G, axes, A_ndim, B_ndim):
    # The adjoint of the operator
    # B |--> np.tensordot(A, B, axes)
    if A_ndim == 0:
        return G * A

    G_axes = onp.arange(onp.ndim(G))
    if type(axes) is int:
        axes = max(axes, 0)
        A_axes = onp.arange(A_ndim)
        return onp.tensordot(A, G, [A_axes[:A_ndim-axes], G_axes[:A_ndim-axes]])
    elif type(axes[0]) is int:
        axes = [axes[0] % A_ndim, axes[1] % B_ndim]
        A_axes = onp.arange(A_ndim)
        return onp.tensordot(A, G, [onp.delete(A_axes, axes[0]), G_axes[:A_ndim-1]])
    else:
        A_axes = onp.arange(A_ndim)
        B_axes = onp.arange(B_ndim)
        summed_axes = [onp.asarray(axes[0]) % A_ndim,
                       onp.asarray(axes[1]) % B_ndim]
        other_axes  = [onp.delete(A_axes, summed_axes[0]),
                       onp.delete(B_axes, summed_axes[1])]
        out = onp.tensordot(A, G, [other_axes[0], G_axes[:len(other_axes[0])]])
        perm = onp.argsort(onp.concatenate(
            (summed_axes[1][onp.argsort(summed_axes[0])], other_axes[1])))
        return onp.transpose(out, perm)

 def energy(self,R1,M0):
     R0=np.tensordot(M0,self.MPO[0],axes=(0,0))
     R0=np.tensordot(R0,np.conj(M0),axes=(2,0))
     R0=np.transpose(R0,[0,2,4,1,3,5])
     energy1=np.squeeze(np.tensordot(R0,R1,axes=([3,4,5],[0,1,2])))
     norm=np.tensordot(M0,np.conj(M0),axes=([0,1,2],[0,1,2]))
     return energy1/norm

def concprof(n0, T, r, phi):

    kb = 8.61733034e-5
    beta = 1.0 / (kb * T)

    ux = strain(r, phi)

    st = np.zeros((3, 3))
    for i in xrange(3):
        for j in xrange(3):
            for k in xrange(3):
                for l in xrange(3):
                    st[i, j] = Cijkl[i, j, k, l] * ux[k, l] + st[i, j]

    F1 = -float(np.tensordot(dpl1, ux))
    F2 = -float(np.tensordot(dpl2, ux))

    chi1 = n0 / (1 - n0) * np.exp(-F1 * beta)
    chi2 = n0 / (1 - n0) * np.exp(-F2 * beta)

    n1 = chi1 / (1 + chi1)
    n2 = chi2 / (1 + chi2)

    n1 = n0 / (np.exp(F1 * beta) * (1 - n0) + n0)
    n2 = n0 / (np.exp(F2 * beta) * (1 - n0) + n0)

    dn = (n1 + n2) - n0

    return dn

    def learn(self, Xd):

        # Initialize a data particle
        X = [Xd]+[self.X[l]*0+self.O[l]
                           for l in range(1, len(self.X))]

        # Alternate gibbs sampler on data and free particles
        for l in (range(1, len(self.X), 2)+range(2, len(self.X), 2))*5:
            self.gibbs(X, l)
        for l in (range(1, len(self.X), 2)+range(0, len(self.X), 2))*1:
            self.gibbs(self.X, l)

        # Parameter update
        self.W[0] += lr*(np.tensordot(X[0]-self.O[0], X[1]-self.O[1], axes=([0],[0])) -
                        np.tensordot(self.X[0]-self.O[0], self.X[1]-self.O[1], axes=([0],[0])))/len(Xd)
        for i in range(1, len(self.W)):
            self.W[i] += lr*(numpy.dot((X[i]-self.O[i]).T,     X[i+1]-self.O[i+1]) -
                             numpy.dot((self.X[i]-self.O[i]).T, self.X[i+1]-self.O[i+1]))/len(Xd)
        for i in range(0, len(self.B)):
            self.B[i] += lr*(X[i]-self.X[i]).mean(axis=0)

        # Reparameterization
        for l in range(0, len(self.B)):
            self.reparamB(X, l)
        for l in range(0, len(self.O)):
            self.reparamO(X, l)

def princ2(cauchyStress,C):
    """
    Calculate the principal values using
    Sig = C:T:s form

    Arguments
    ---------
    cauchyStress - the linearly transformed stress Sigma
    C (6x6) tensor for linear transformation

    Returns
    -------
    S1
    S2
    S3
    """
    sqrt=np.sqrt
    s = cauchyStress.copy()
    T = cpb_lib.returnT()
    Sig = np.tensordot(C,np.tensordot(T,s,axes=(1,0)),axes=(1,0))

    Sxx,Syy,Szz,Sxy = Sig[0], Sig[1], Sig[2], Sig[5]

    S1 = 0.5*(Sxx + Syy + sqrt( (Sxx-Syy)**2 + 4*Sxy**2))
    S2 = 0.5*(Sxx + Syy - sqrt( (Sxx-Syy)**2 + 4*Sxy**2))
    S3 = Szz
    return S1, S2, S3

def discretizedGaussian(amp, mu, cov, grid):
    '''Convenience method for discretized Gaussian evaluation'''
    #eigenvalue decomposition of precision matrix
    P = la.inv(cov)     #precision matrix
    evl, M = la.eig(P)
    
    #assert np.allclose(np.diag(evl), iM.dot(cov).dot(M))
    #check if covariance is positive definite
    if np.any(evl < 0):
        raise ValueError('Covariance matrix should be positive definite')
    
    
    #make column vector for arithmetic
    mu = np.array(mu, ndmin=grid.ndim, dtype=float).T
    evl = np.array(evl, ndmin=grid.ndim, dtype=float).T

    xm = grid - mu #(2,...) shape
    
    #return M, xm
    
    f = np.sqrt(2 / evl)
    pf = np.sqrt(np.pi / evl)
    td0 = np.tensordot(M, xm + 0.5, 1)
    td1 = np.tensordot(M, xm - 0.5, 1)
    w = pf*(erf(f*td0) - erf(f*td1))
    
    return amp * np.prod(w, axis=0)

 def get_chisquareds(self, n):
     dd = self.get_datum(n)
     if n == 0:
         foo = np.sum(np.sum(np.tensordot(dd, self.get_ivarts(), axes=(1,1)) * dd[:,None,:],
                             axis=2), axis=0)
     return np.sum(np.sum(np.tensordot(dd, self.get_ivarts(), axes=(1,1)) * dd[:,None,:],
                          axis=2), axis=0) # should be length T

 def update(self):
     check_inputs_synchronized(self)
     
     y_dot = self.input.y_dot
     y = self.input.y
     u = self.input.u
     T = self.T
     
     check_multiple([
         ('shape(x),(array[HxW]|array[N])', y),
         ('shape(x),(array[HxW]|array[N])', y_dot),
         ('array[K]', u),
         ('array[Kx2xHxW]|array[Kx1xN]', T), # TODO: use references
     ])
     K = u.size
     
     # update covariance of gradients
     gy = generalized_gradient(y)
     
     Tgy = np.tensordot([1, 1], T * gy , axes=(0, 1)) 
             
     y_dot_pred = np.tensordot(u, Tgy, axes=(0, 0))
      
     assert y_dot_pred.ndim == 2
     
     error = np.maximum(0, -y_dot_pred * y_dot)
     # error = np.abs(np.sign(y_dot_pred) - np.sign(y_dot))
     
     self.output.y_dot_pred = y_dot_pred
     self.output.error = error
     self.output.error_sensel = np.sum(error, axis=0)

def dev_trapz(ds):
	d = ds.variables['development'][:]/1e4
	t = ds.variables['time'][:]
	t_floor = np.floor(t).astype(int)
	#now set up a loop to solve the index for each whole integer or iterate up
	#we can figure out if we have a leap year or not also by looking at the max day
	#i.e. np.max(t_floor) == 365 , means this is a leap year dataset

	max_days = np.max(t_floor) + 1
	z,m,n = d.shape
	daily_dev = np.zeros((max_days, m,n))
	start = 0 #first integration point
	k_array = np.zeros(max_days)
	for k in range(max_days):
		ones = np.ones((m,n))
		if k == max_days-1:
			tm = np.tensordot(t[start:], ones, axes=0)
			daily_dev[k] = np.trapz(d[start:],tm,axis=0) 
		else:
			end = np.argmax(t_floor>k)
			tm = np.tensordot(t[start:end+1], ones, axes=0) #plus one to integrate the end point
			daily_dev[k] = np.trapz(d[start:end+1],tm,axis=0) 
			#define the new start
			start = end
		k_array[k] = k+1 #counter for each day
	return(daily_dev, k_array)

def four_index_transform_cholesky(ao_integrals, orb0, orb1=None, method='tensordot'):
    """Perform four index transformation on a Cholesky-decomposed four-index object.

    Parameters
    ----------
    oa_integrals : np.ndarray, shape=(nvec, nbasis, nbasis)
        Cholesky decomposition of four-index object in the AO basis.
    orb0
        A Orbitals object with molecular orbitals.
    orb1
        Can be provided to transform the second index differently.
    method
        Either ``einsum`` or ``tensordot`` (default).
    """
    if orb1 is None:
        orb1 = orb0
    result = np.zeros(ao_integrals.shape)
    if method == 'einsum':
        result = np.einsum('ai,kac->kic', orb0.coeffs, ao_integrals)
        result = np.einsum('cj,kic->kij', orb1.coeffs, result)
    elif method == 'tensordot':
        result = np.tensordot(ao_integrals, orb0.coeffs, axes=([1],[0]))
        result = np.tensordot(result, orb1.coeffs, axes=([1],[0]))
    else:
        raise ValueError('The method must either be \'einsum\' or \'tensordot\'.')
    return result

 def d2y_SS(z_i):
     DF,df,Hf = self.DF(z_i),self.df(z_i),self.Hf(z_i)
     d = self.get_d(z_i)
     DFi = DF[n:]
     return np.tensordot(np.linalg.inv(DFi[:,y] + DFi[:,e].dot(df) + DFi[:,v].dot(Ivy)),
         -self.HFhat[S,S](z_i) - np.tensordot(DFi[:,e],quadratic_dot(Hf,d[y,S],d[y,S]),1)
             , axes=1)

def test_solve_fredholm_reconstr_ac():
    """
        here we see that the reconstruction quality is independent of the integration weights

        differences occur when checking validity of the interpolated time continuous Fredholm equation
    """
    _WC_ = 2
    def lac(t):
        return np.exp(- np.abs(t) - 1j*_WC_*t)
    t_max = 10
    tol = 2e-10
    for ng in range(11,500,30):
        t, w = sp.method_kle.get_mid_point_weights_times(t_max, ng)
        r = lac(t.reshape(-1,1)-t.reshape(1,-1))
        _eig_val, _eig_vec = sp.method_kle.solve_hom_fredholm(r, w)
        _eig_vec_ast = np.conj(_eig_vec)  # (N_gp, N_ev)
        tmp = _eig_val.reshape(1, -1) * _eig_vec  # (N_gp, N_ev)
        recs_bcf = np.tensordot(tmp, _eig_vec_ast, axes=([1], [1]))
        rd = np.max(np.abs(recs_bcf - r) / np.abs(r))
        assert rd < tol, "rd={} >= {}".format(rd, tol)

        t, w = sp.method_kle.get_simpson_weights_times(t_max, ng)
        r = lac(t.reshape(-1, 1) - t.reshape(1, -1))
        _eig_val, _eig_vec = sp.method_kle.solve_hom_fredholm(r, w)
        _eig_vec_ast = np.conj(_eig_vec)  # (N_gp, N_ev)
        tmp = _eig_val.reshape(1, -1) * _eig_vec  # (N_gp, N_ev)
        recs_bcf = np.tensordot(tmp, _eig_vec_ast, axes=([1], [1]))
        rd = np.max(np.abs(recs_bcf - r) / np.abs(r))
        assert rd < tol, "rd={} >= {}".format(rd, tol)