Python scipy.kron函数代码示例

    def K(self):
        if self.dim > _MAX_DIM:
            raise TooExpensiveOperationError(msg_too_expensive_dim(my_name(),
                                                                   _MAX_DIM))

        rv = sp.kron(self.Cg.K(), self.R) + sp.kron(self.Cn.K(), sp.eye(self.dim_r))
        return rv

    def _LMLgrad_s(self,hyperparams,debugging=False):
        """
        evaluate gradients with respect to covariance matrix Sigma
        """
        try:
            KV = self.get_covariances(hyperparams, debugging=debugging)
        except LA.LinAlgError:
            LG.error('linalg exception in _LMLgrad_x_sigma')
            return {'X_s':SP.zeros(hyperparams['X_s'].shape)}

        Si = 1./KV['Stilde_os']
        SS = SP.dot(unravel(Si,self.n,self.t).T,KV['Stilde_o'])
        USU = SP.dot(KV['USi_c'],KV['Utilde_s'])        
        Yhat = unravel(Si * ravel(KV['UYtildeU_os']),self.n,self.t)
        RV = {}
        
        if 'X_s' in hyperparams:
            USUY = SP.dot(USU,Yhat.T)
            USUYSYUSU = SP.dot(USUY,(KV['Stilde_o']*USUY).T)
        
            LMLgrad = SP.zeros((self.t,self.covar_s.n_dimensions))
            LMLgrad_det = SP.zeros((self.t,self.covar_s.n_dimensions))
            LMLgrad_quad = SP.zeros((self.t,self.covar_s.n_dimensions))
        
            for d in xrange(self.covar_s.n_dimensions):
                Kd_grad = self.covar_s.Kgrad_x(hyperparams['covar_s'],d)
                UsU = SP.dot(Kd_grad.T,USU)*USU
                LMLgrad_det[:,d] = SP.dot(UsU,SS.T)
                # calculate gradient of squared form
                LMLgrad_quad[:,d] = -(USUYSYUSU*Kd_grad).sum(0)
            LMLgrad = LMLgrad_det + LMLgrad_quad
            RV['X_s'] = LMLgrad
            
            if debugging:
                _LMLgrad = SP.zeros((self.t,self.covar_s.n_dimensions))
                for t in xrange(self.t):
                    for d in xrange(self.covar_s.n_dimensions):
                        Kgrad_x = self.covar_s.Kgrad_x(hyperparams['covar_s'],d,t)
                        Kgrad_x = SP.kron(Kgrad_x,KV['K_o'])
                        _LMLgrad[t,d] = 0.5*(KV['W']*Kgrad_x).sum()

                assert SP.allclose(LMLgrad,_LMLgrad,rtol=1E-3,atol=1E-2), 'ouch, something is wrong: %.2f'%LA.norm(LMLgrad-_LMLgrad)

        if 'covar_s' in hyperparams:
            theta = SP.zeros(len(hyperparams['covar_s']))
            for i in range(len(theta)):
                Kgrad_s = self.covar_s.Kgrad_theta(hyperparams['covar_s'],i)
                UdKU = SP.dot(USU.T, SP.dot(Kgrad_s, USU))
                SYUdKU = SP.dot(UdKU,KV['Stilde_o'] * Yhat.T)
                LMLgrad_det = SP.sum(Si*SP.kron(SP.diag(UdKU),KV['Stilde_o']))
                LMLgrad_quad = -(Yhat.T*SYUdKU).sum()
                LMLgrad = 0.5*(LMLgrad_det + LMLgrad_quad)
                theta[i] = LMLgrad

                if debugging:
                    Kd = SP.kron(Kgrad_s, KV['K_o'])
                    _LMLgrad = 0.5 * (KV['W']*Kd).sum()
                    assert SP.allclose(LMLgrad,_LMLgrad,rtol=1E-3,atol=1E-2), 'ouch, something is wrong: %.2f'%LA.norm(LMLgrad-_LMLgrad)
            RV['covar_s'] = theta
        return RV

    def _LMLgrad_covar_debug(self,covar):

        assert self.N*self.P<2000, 'gp2kronSum:: N*P>=2000'

        y  = SP.reshape(self.Y,(self.N*self.P), order='F') 

        K  = SP.kron(self.Cg.K(),self.XX)
        K += SP.kron(self.Cn.K()+self.offset*SP.eye(self.P),SP.eye(self.N))

        cholK = LA.cholesky(K).T
        Ki  = LA.cho_solve((cholK,True),SP.eye(y.shape[0]))
        Kiy   = LA.cho_solve((cholK,True),y)

        if covar=='Cr':     n_params = self.Cr.getNumberParams()
        elif covar=='Cg':   n_params = self.Cg.getNumberParams()
        elif covar=='Cn':   n_params = self.Cn.getNumberParams()

        RV = SP.zeros(n_params)

        for i in range(n_params):
            #0. calc grad_i
            if covar=='Cg':
                C   = self.Cg.Kgrad_param(i)
                Kgrad  = SP.kron(C,self.XX)
            elif covar=='Cn':
                C   = self.Cn.Kgrad_param(i)
                Kgrad  = SP.kron(C,SP.eye(self.N))

            #1. der of log det
            RV[i]  = 0.5*(Ki*Kgrad).sum()
            
            #2. der of quad form
            RV[i] -= 0.5*(Kiy*SP.dot(Kgrad,Kiy)).sum()

        return RV

	def _update_indicator(self,K,L):
		""" update the indicator """
		_update = {'term': self.n_terms*SP.ones((K,L)).T.ravel(),
					'row': SP.kron(SP.arange(K)[:,SP.newaxis],SP.ones((1,L))).T.ravel(),
					'col': SP.kron(SP.ones((K,1)),SP.arange(L)[SP.newaxis,:]).T.ravel()} 
		for key in _update.keys():
			self.indicator[key] = SP.concatenate([self.indicator[key],_update[key]])

 def varMLE(self):
     """ calculate inverse of fisher information """
     self._update_cache()
     Sr = {}
     Sr['Cg'] = self.cache['Srstar']
     Sr['Cn'] = SP.ones(self.N)
     n_params = self.Cg.getNumberParams()+self.Cn.getNumberParams()
     fisher = SP.zeros((n_params,n_params))
     header = SP.zeros((n_params,n_params),dtype='|S10')
     C1  = SP.zeros((self.P,self.P))
     C2  = SP.zeros((self.P,self.P))
     idx1 = 0
     for key1 in ['Cg','Cn']:
         for key1_p1 in range(self.P):
             for key1_p2 in range(key1_p1,self.P):
                 C1[key1_p1,key1_p2] = C1[key1_p2,key1_p1] = 1
                 LCL1 = SP.dot(self.cache['Lc'],SP.dot(C1,self.cache['Lc'].T))
                 CSr1 = SP.kron(Sr[key1][:,SP.newaxis,SP.newaxis],LCL1[SP.newaxis,:])
                 DCSr1 = self.cache['D'][:,:,SP.newaxis]*CSr1
                 idx2 = 0
                 for key2 in ['Cg','Cn']:
                     for key2_p1 in range(self.P):
                         for key2_p2 in range(key2_p1,self.P):
                             C2[key2_p1,key2_p2] = C2[key2_p2,key2_p1] = 1
                             LCL2 = SP.dot(self.cache['Lc'],SP.dot(C2,self.cache['Lc'].T))
                             CSr2 = SP.kron(Sr[key2][:,SP.newaxis,SP.newaxis],LCL2[SP.newaxis,:])
                             DCSr2 = self.cache['D'][:,:,SP.newaxis]*CSr2
                             fisher[idx1,idx2] = 0.5*(DCSr1*DCSr2).sum()
                             header[idx1,idx2] = '%s%d%d_%s%d%d'%(key1,key1_p1,key1_p1,key2,key2_p1,key2_p2)
                             C2[key2_p1,key2_p2] = C2[key2_p2,key2_p1] = 0
                             idx2+=1
                 C1[key1_p1,key1_p2] = C1[key1_p2,key1_p1] = 0
                 idx1+=1
     RV = LA.inv(fisher)
     return RV,header

    def LMLdebug(self):
        """
        LML function for debug
        """
        assert self.N*self.P<5000, 'gp2kronSum:: N*P>=5000'

        y = SP.reshape(self.Y,(self.N*self.P), order='F') 
        V = SP.kron(SP.eye(self.P),self.F)

        XX = SP.dot(self.Xr,self.Xr.T)
        K  = SP.kron(self.Cr.K(),XX)
        K += SP.kron(self.Cn.K()+self.offset*SP.eye(self.P),SP.eye(self.N))

        # inverse of K
        cholK = LA.cholesky(K)
        Ki = LA.cho_solve((cholK,False),SP.eye(self.N*self.P))

        # Areml and inverse
        Areml = SP.dot(V.T,SP.dot(Ki,V))
        cholAreml = LA.cholesky(Areml)
        Areml_i = LA.cho_solve((cholAreml,False),SP.eye(self.K*self.P))

        # effect sizes and z
        b = SP.dot(Areml_i,SP.dot(V.T,SP.dot(Ki,y)))
        z = y-SP.dot(V,b)
        Kiz = SP.dot(Ki,z)

        # lml
        lml  = y.shape[0]*SP.log(2*SP.pi)
        lml += 2*SP.log(SP.diag(cholK)).sum()
        lml += 2*SP.log(SP.diag(cholAreml)).sum()
        lml += SP.dot(z,Kiz)
        lml *= 0.5

        return lml

 def diag_Ctilde_o_Sr(self, i):
     if i < self.Cg.getNumberParams():
         r = sp.kron(sp.diag(self.LcGradCgLc(i)), self.Sr())
     else:
         _i = i - self.Cg.getNumberParams()
         r = sp.kron(sp.diag(self.LcGradCnLc(_i)), sp.ones(self.dim_r))
     return r

def sqrtm3(X):
    M = sp.copy(X)
    m, fb, fe = block_structure(M)
    n = M.shape[0]
    for i in range(0,m):
        M[fb[i]:fe[i],fb[i]:fe[i]] = twobytworoot(M[fb[i]:fe[i],fb[i]:fe[i]])
        #print M

    for j in range(1,m):
        for i in range(0,m-j):
            #print M[fb[i]:fe[i],fb[JJ]:fe[JJ]]
            JJ = i+j
            Tnoto = M[fb[i]:fe[i],fb[JJ]:fe[JJ]] #dopo togliere il copy
            #print "Tnot: "
            #print Tnoto
            for k in range(i+1,JJ):
                Tnoto -= (M[fb[i]:fe[i],fb[k]:fe[k]]).dot(M[fb[k]:fe[k],fb[JJ]:fe[JJ]])
                #print M[fb[i]:fe[i],fb[k]:fe[k]]
                #print M[fb[k]:fe[k],fb[JJ]:fe[JJ]]

            if((M[fb[i]:fe[i],fb[JJ]:fe[JJ]]).shape==(1,1)):
                #print "forma 1"
                #print M[fb[i]:fe[i],fb[JJ]:fe[JJ]]           #  Uij
                #print M[fb[i]:fe[i],fb[i]:fe[i]]               #  Uii
                #print M[fb[JJ]:fe[JJ],fb[JJ]:fe[JJ]]       #  Ujj
                M[fb[i]:fe[i],fb[JJ]:fe[JJ]] = Tnoto/(M[fb[i]:fe[i],fb[i]:fe[i]] + M[fb[JJ]:fe[JJ],fb[JJ]:fe[JJ]])

            else:
                Uii = M[fb[i]:fe[i],fb[i]:fe[i]]
                Ujj = M[fb[JJ]:fe[JJ],fb[JJ]:fe[JJ]]
                shapeUii = Uii.shape[0]
                shapeUjj = Ujj.shape[0]
                """
                print "------------"
                print Tnoto
                print Tnoto.shape
                print sp.kron(sp.eye(shapeUjj),Uii)
                print sp.kron(Ujj.T,sp.eye(shapeUii))
                print Tnoto
                """
                #M[fb[i]:fe[i],fb[JJ]:fe[JJ]] = sp.linalg.solve_sylvester(Uii, Ujj, Tnoto)

                """
                x, scale, info = dtrsyl(Uii, Ujj, Tnoto

                if (scale==1.0):
                     = x

                else:
                    M[fb[i]:fe[i],fb[JJ]:fe[JJ]] = x*scale
                    print "scale!=0"
                """
                Tnoto = Tnoto.reshape((shapeUii*shapeUjj),1,order="F")
                M[fb[i]:fe[i],fb[JJ]:fe[JJ]] = \
                linalg.solve(sp.kron(sp.eye(shapeUjj),Uii) +
                sp.kron(Ujj.T,sp.eye(shapeUii)),
                Tnoto).reshape(shapeUii,shapeUjj,order="F")


    return M

def AlphaBetaCoeffs_old(n, a, b):
    " Construct the alpha and beta coefficient matrices. "
    Z = sp.matrix([[1,0],[0,-1]])
    I = sp.identity(2)
    alpha = sp.zeros((2**n,2**n))
    beta = sp.zeros((2**n,2**n))
    m1 = []
    m2 = []
    for i in range(0,n):
        for m in range(0,n-1): m1.append(I)
        m1.insert(i, Z)
        temp1 = m1[0]
        m1.pop(0)

        while (len(m1) != 0):
            temp1 = sp.kron(temp1, m1[0])
            m1.pop(0)
        alpha += temp1*a[i]
        for j in range(i+1, n):
            for m in range(0, n-2): m2.append(I)
            m2.insert(i, Z)
            m2.insert(j, Z)
            temp2 = m2[0]
            m2.pop(0)
            while (len(m2) != 0):
                temp2 = sp.kron(temp2, m2[0])
                m2.pop(0)
            beta += (temp2)*b[i,j]
    return [alpha, beta]

    def _LMLgrad_covar(self,hyperparams,debugging=False):
        """
        evaluates the gradient of the log marginal likelihood with respect to the
        hyperparameters of the covariance function
        """
        try:
            KV = self.get_covariances(hyperparams,debugging=debugging)
        except LA.LinAlgError:
            LG.error('linalg exception in _LMLgrad_covar')
            return {'covar_r':SP.zeros(len(hyperparams['covar_r'])),'covar_c':SP.zeros(len(hyperparams['covar_c'])),'covar_r':SP.zeros(len(hyperparams['covar_r']))}
        except ValueError:
            LG.error('value error in _LMLgrad_covar')
            return {'covar_r':SP.zeros(len(hyperparams['covar_r'])),'covar_c':SP.zeros(len(hyperparams['covar_c'])),'covar_r':SP.zeros(len(hyperparams['covar_r']))}
 
        RV = {}
        Si = unravel(1./KV['S'],self.n,self.t)

        if 'covar_r' in hyperparams:
            theta = SP.zeros(len(hyperparams['covar_r']))
            for i in range(len(theta)):
                Kgrad_r = self.covar_r.Kgrad_theta(hyperparams['covar_r'],i)
                d=(KV['U_r']*SP.dot(Kgrad_r,KV['U_r'])).sum(0)
                LMLgrad_det = SP.dot(d,SP.dot(Si,KV['S_c']))
                UdKU = SP.dot(KV['U_r'].T,SP.dot(Kgrad_r,KV['U_r']))
                SYUdKU = SP.dot(UdKU,(KV['Ytilde']*SP.tile(KV['S_c'][SP.newaxis,:],(self.n,1))))
                LMLgrad_quad = - (KV['Ytilde']*SYUdKU).sum()
                LMLgrad = 0.5*(LMLgrad_det + LMLgrad_quad)
                theta[i] = LMLgrad

                if debugging:
                    Kd = SP.kron(KV['K_c'], Kgrad_r)
                    _LMLgrad = 0.5 * (KV['W']*Kd).sum()
                    assert SP.allclose(LMLgrad,_LMLgrad), 'ouch, gradient is wrong for covar_r'
                    
            RV['covar_r'] = theta

        if 'covar_c' in hyperparams:
            theta = SP.zeros(len(hyperparams['covar_c']))
            for i in range(len(theta)):
                Kgrad_c = self.covar_c.Kgrad_theta(hyperparams['covar_c'],i)

                d=(KV['U_c']*SP.dot(Kgrad_c,KV['U_c'])).sum(0)
                LMLgrad_det = SP.dot(KV['S_r'],SP.dot(Si,d))

                UdKU = SP.dot(KV['U_c'].T,SP.dot(Kgrad_c,KV['U_c']))
                SYUdKU = SP.dot((KV['Ytilde']*SP.tile(KV['S_r'][:,SP.newaxis],(1,self.t))),UdKU.T)
                LMLgrad_quad = -SP.sum(KV['Ytilde']*SYUdKU)
                LMLgrad = 0.5*(LMLgrad_det + LMLgrad_quad)
                theta[i] = LMLgrad
            
                if debugging:
                    Kd = SP.kron(Kgrad_c, KV['K_r'])
                    _LMLgrad = 0.5 * (KV['W']*Kd).sum()
                    assert SP.allclose(LMLgrad,_LMLgrad), 'ouch, gradient is wrong for covar_c'
                    
                RV['covar_c'] = theta

        return RV

def get_linds(N, eps):
    #Lindblad operators must have same range as Hamiltonian terms. In this case they are nearest-neighbour.
    Sp1 = (sp.kron(Sp, sp.eye(2))).reshape(2, 2, 2, 2)
    Sm2 = (sp.kron(sp.eye(2), Sm)).reshape(2, 2, 2, 2)
    
    L1 = (1, sp.sqrt(eps) * Sp1)
    L2 = (N-1, sp.sqrt(eps) * Sm2)
    
    return [L1, L2]

 def _fit(self, type, vc=False):
     #2. init
     if type=='null':
         self.gp[type].covar.Cg.setCovariance(0.5*self.covY)
         self.gp[type].covar.Cn.setCovariance(0.5*self.covY)
     elif type=='full':
         Cg0_K = self.gp['null'].covar.Cg.K()
         Cn0_K = self.gp['null'].covar.Cn.K()
         self.gp[type].covar.Cr.setCovariance((Cn0_K+Cg0_K)/3.)
         self.gp[type].covar.Cg.setCovariance(2.*Cg0_K/3.)
         self.gp[type].covar.Cn.setCovariance(2.*Cn0_K/3.)
     elif type=='block':
         Crf_K = self.gp['full'].covar.Cr.K()
         Cnf_K = self.gp['full'].covar.Cn.K()
         self.gp[type].covar.Cr.scale = sp.mean(Crf_K)
         self.gp[type].covar.Cn.setCovariance(Cnf_K)
     elif type=='rank1':
         Crf_K = self.gp['full'].covar.Cr.K()
         Cnf_K = self.gp['full'].covar.Cn.K()
         self.gp[type].covar.Cr.setCovariance(Crf_K)
         self.gp[type].covar.Cn.setCovariance(Cnf_K)
     else:
         print('poppo')
     self.gp[type].optimize(factr=self.factr, verbose=False)
     RV = {'Cg': self.gp[type].covar.Cg.K(),
             'Cn': self.gp[type].covar.Cn.K(),
             'LML': sp.array([self.gp[type].LML()]),
             'LMLgrad': sp.array([sp.mean((self.gp[type].LML_grad()['covar'])**2)])}
     if type=='null':        RV['Cr'] = sp.zeros(RV['Cg'].shape) 
     else:                   RV['Cr'] = self.gp[type].covar.Cr.K() 
     if vc:
         # tr(P CoR) = tr(C)tr(R) - tr(Ones C) tr(Ones R) / float(NP)
         #           = tr(C)tr(R) - C.sum() * R.sum() / float(NP)
         trRr = (self.Xr**2).sum()
         var_r = sp.trace(RV['Cr'])*trRr / float(self.Y.size-1)
         var_g = sp.trace(RV['Cg'])*self.trRg / float(self.Y.size-1)
         var_n = sp.trace(RV['Cn'])*self.Y.shape[0] 
         var_n-= RV['Cn'].sum() / float(RV['Cn'].shape[0])
         var_n/= float(self.Y.size-1) 
         RV['var'] = sp.array([var_r, var_g, var_n])
         if 0 and self.Y.size<5000:
             pdb.set_trace()
             Kr = sp.kron(RV['Cr'], sp.dot(self.Xr, self.Xr.T))
             Kn = sp.kron(RV['Cn'], sp.eye(self.Y.shape[0]))
             _var_r = sp.trace(Kr-Kr.mean(0)) / float(self.Y.size-1)
             _var_n = sp.trace(Kn-Kn.mean(0)) / float(self.Y.size-1)
             _var = sp.array([_var_r, var_g, _var_n])
             print(((_var-RV['var'])**2).mean())
         if type=='full':
             # calculate within region vcs 
             Cr_block = sp.mean(RV['Cr']) * sp.ones(RV['Cr'].shape)
             Cr_rank1 = lowrank_approx(RV['Cr'], rank=1)
             var_block = sp.trace(Cr_block)*trRr / float(self.Y.size-1)
             var_rank1 = sp.trace(Cr_rank1)*trRr / float(self.Y.size-1)
             RV['var_r'] = sp.array([var_block, var_rank1-var_block, var_r-var_rank1])
     return RV

def distanz(x, y=None):
    r"""
    Calculate the distanz between two colon vectors

    Parameters
    ----------
    x : ndarray
        First colon vector
    y : ndarray
        Second colon vector

    Returns
    -------
    d : ndarray
        Distance between x and y

    Examples
    --------
    >>> import numpy as np
    >>> from pygsp import utils
    >>> x = np.random.rand(16)
    >>> y = np.random.rand(16)
    >>> distanz = utils.distanz(x, y)

    """
    try:
        x.shape[1]
    except IndexError:
        x = x.reshape(1, x.shape[0])

    if y is None:
        y = x

    else:
        try:
            y.shape[1]
        except IndexError:
            y = y.reshape(1, y.shape[0])

    rx, cx = x.shape
    ry, cy = y.shape

    # Size verification
    if rx != ry:
        raise("The sizes of x and y do not fit")

    xx = (x*x).sum(axis=0)
    yy = (y*y).sum(axis=0)
    xy = np.dot(x.T, y)

    d = abs(sp.kron(sp.ones((cy, 1)), xx).T +
            sp.kron(sp.ones((cx, 1)), yy) - 2*xy)

    return np.sqrt(d)

    def test_inv(self):
        C = self.C

        L = sp.kron(C.Lc(), sp.eye(C.dim_r))
        W = sp.kron(C.Wc(), C.Wr())
        WdW = sp.dot(W.T, C.d()[:, sp.newaxis] * W)
        I_WdW = sp.eye(C.dim_c * C.dim_r) - WdW
        inv1 = sp.dot(L.T, sp.dot(I_WdW, L))

        inv2 = C.inv()
        np.testing.assert_array_almost_equal(inv1, inv2)

 def diag_Ctilde_o_Sr(self, i):
     np_r = self.Cr.getNumberParams()
     np_g = self.Cg.getNumberParams()
     if i < np_r:
         r = sp.kron(sp.diag(self.LcGradCrLc(i)), self.diagWrWr())
     elif i < (np_r + np_g):
         _i = i - np_r
         r = sp.kron(sp.diag(self.LcGradCgLc(_i)), self.Sr())
     else:
         _i = i - np_r - np_g
         r = sp.kron(sp.diag(self.LcGradCnLc(_i)), sp.ones(self.dim_r))
     return r