python类linalg()的实例源码

def logdet(C,eps=1e-6,safe=0):
    '''
    Logarithm of the determinant of a matrix
    Works only with real-valued positive definite matrices
    '''
    try:
        return 2.0*np.sum(np.log(np.diag(chol(C))))
    except numpy.linalg.linalg.LinAlgError:
        if safe: C = check_covmat(C,eps=eps)
        w = np.linalg.eigh(C)[0]
        w = np.real(w)
        w[w<eps]=eps
        det = np.sum(np.log(w))
        return det

def nufft_T(N, J, K, alpha, beta):
    '''
     The Equation (29) and (26) in Fessler and Sutton 2003.
     Create the overlapping matrix CSSC (diagonal dominent matrix)
     of J points and find out the pseudo-inverse of CSSC '''

#     import scipy.linalg
    L = numpy.size(alpha) - 1
#     print('L = ', L, 'J = ',J, 'a b', alpha,beta )
    cssc = numpy.zeros((J, J))
    [j1, j2] = numpy.mgrid[1:J + 1, 1:J + 1]
    overlapping_mat = j2 - j1
    for l1 in range(-L, L + 1):
        for l2 in range(-L, L + 1):
            alf1 = alpha[abs(l1)]
#             if l1 < 0: alf1 = numpy.conj(alf1)
            alf2 = alpha[abs(l2)]
#             if l2 < 0: alf2 = numpy.conj(alf2)
            tmp = overlapping_mat + beta * (l1 - l2)

            tmp = dirichlet(1.0 * tmp / (1.0 * K / N))
            cssc = cssc + alf1 * alf2 * tmp

    return mat_inv(cssc)

def update_per_row(self, y_i, phi_i, J, mu0, c, v, r_prev_i, u_prev_i, phi_r_i, phi_u):
        # Params:
        #   J - column indices

        nnz_i = len(J)
        residual_i = y_i - mu0 - c[J]
        prior_Phi = np.diag(np.concatenate(([phi_r_i], phi_u)))
        v_T = np.hstack((np.ones((nnz_i, 1)), v[J, :]))
        post_Phi_i = prior_Phi + \
                     np.dot(v_T.T,
                            np.tile(phi_i[:, np.newaxis], (1, 1 + self.num_factor)) * v_T)  # Weighted sum of v_j * v_j.T
        post_mean_i = np.squeeze(np.dot(phi_i * residual_i, v_T))
        C, lower = scipy.linalg.cho_factor(post_Phi_i)
        post_mean_i = scipy.linalg.cho_solve((C, lower), post_mean_i)
        # Generate Gaussian, recycling the Cholesky factorization from the posterior mean computation.
        ru_i = math.sqrt(1 - self.relaxation ** 2) * scipy.linalg.solve_triangular(C, np.random.randn(len(post_mean_i)),
                                                                                   lower=lower)
        ru_i += post_mean_i + self.relaxation * (post_mean_i - np.concatenate(([r_prev_i], u_prev_i)))
        r_i = ru_i[0]
        u_i = ru_i[1:]

        return r_i, u_i

def update_per_col(self, y_j, phi_j, I, mu0, r, u, c_prev_j, v_prev_j, phi_c_j, phi_v):

        prior_Phi = np.diag(np.concatenate(([phi_c_j], phi_v)))
        nnz_j = len(I)
        residual_j = y_j - mu0 - r[I]
        u_T = np.hstack((np.ones((nnz_j, 1)), u[I, :]))
        post_Phi_j = prior_Phi + \
                     np.dot(u_T.T,
                            np.tile(phi_j[:, np.newaxis], (1, 1 + self.num_factor)) * u_T)  # Weighted sum of u_i * u_i.T
        post_mean_j = np.squeeze(np.dot(phi_j * residual_j, u_T))
        C, lower = scipy.linalg.cho_factor(post_Phi_j)
        post_mean_j = scipy.linalg.cho_solve((C, lower), post_mean_j)
        # Generate Gaussian, recycling the Cholesky factorization from the posterior mean computation.
        cv_j = math.sqrt(1 - self.relaxation ** 2) * scipy.linalg.solve_triangular(C, np.random.randn(len(post_mean_j)),
                                                                              lower=lower)
        cv_j += post_mean_j + self.relaxation * (post_mean_j - np.concatenate(([c_prev_j], v_prev_j)))
        c_j = cv_j[0]
        v_j = cv_j[1:]

        return c_j, v_j

def test_scipy(pyi_builder):
    pyi_builder.test_source(
        """
        from distutils.version import LooseVersion

        # Test top-level SciPy importability.
        import scipy
        from scipy import *

        # Test hooked SciPy modules.
        import scipy.io.matlab
        import scipy.sparse.csgraph

        # Test problematic SciPy modules.
        import scipy.linalg
        import scipy.signal

        # SciPy >= 0.16 privatized the previously public "scipy.lib" package as
        # "scipy._lib". Since this package is problematic, test its
        # importability regardless of SciPy version.
        if LooseVersion(scipy.__version__) >= LooseVersion('0.16.0'):
            import scipy._lib
        else:
            import scipy.lib
        """)

def matrixMult(A, B):
    try:
        linalg.blas
    except AttributeError:
        return np.dot(A, B)

    if not A.flags['F_CONTIGUOUS']:
        AA = A.T
        transA = True
    else:
        AA = A
        transA = False

    if not B.flags['F_CONTIGUOUS']:
        BB = B.T
        transB = True
    else:
        BB = B
        transB = False

    return linalg.blas.dgemm(alpha=1., a=AA, b=BB, trans_a=transA, trans_b=transB)

def factor(X, rho):
    """
    computes cholesky factorization of the kernel K = 1/rho*XX^T + I

    Input:
    X design matrix: n_s x n_f (we assume n_s << n_f)
    rho: regularizaer

    Output:
    L  lower triangular matrix
    U  upper triangular matrix
    """
    n_s, n_f = X.shape
    K = 1 / rho * scipy.dot(X, X.T) + scipy.eye(n_s)
    U = linalg.cholesky(K)
    return U

def FeCalc(self, Aoptimality=True):
        '''
        Compute estimation efficiency.

        :param Aoptimality: Kind of optimality to optimize, A- or D-optimality
        :type Aoptimality: boolean
        '''
        try:
            invM = scipy.linalg.inv(self.X)
        except scipy.linalg.LinAlgError:
            try:
                invM = scipy.linalg.pinv(self.X)
            except np.linalg.linalg.LinAlgError:
                invM = np.nan
        sys.exc_clear()
        invM = np.array(invM)
        st1 = np.dot(self.CX, invM)
        CMC = np.dot(st1, t(self.CX))
        if Aoptimality == True:
            self.Fe = float(self.CX.shape[0] / np.matrix.trace(CMC))
        else:
            self.Fe = float(np.linalg.det(CMC)**(-1 / len(self.C)))
        self.Fe = self.Fe / self.experiment.FeMax
        return self

def FdCalc(self, Aoptimality=True):
        '''
        Compute detection power.

        :param Aoptimality: Kind of optimality to optimize: A- or D-optimality
        :type Aoptimality: boolean
        '''
        try:
            invM = scipy.linalg.inv(self.Z)
        except scipy.linalg.LinAlgError:
            try:
                invM = scipy.linalg.pinv(self.Z)
            except np.linalg.linalg.LinAlgError:
                invM = np.nan
        sys.exc_clear()
        invM = np.array(invM)
        CMC = np.matrix(self.C) * invM * np.matrix(t(self.C))
        if Aoptimality == True:
            self.Fd = float(len(self.C) / np.matrix.trace(CMC))
        else:
            self.Fd = float(np.linalg.det(CMC)**(-1 / len(self.C)))
        self.Fd = self.Fd / self.experiment.FdMax
        return self

def decompose(P):
    M = P[:3, :3]
    T = P[:3, 3]

    K, R = scipy.linalg.rq(M)

    for i in range(2):
        if K[i,i] < 0:
            K[:, i] *= -1
            R[i, :] *= -1

    if K[2,2] > 0:
        K[:, 2] *= -1
        R[2, :] *= -1

    if det(R) < 0:
        R *= -1

    T = linalg.inv(dot(K, -R)).dot(T.reshape((3, 1)))
    K /= -K[2,2]

    return K, R, T

def _gauss_from_coefficients_numpy(alpha, beta):
    assert isinstance(alpha, numpy.ndarray)
    assert isinstance(beta, numpy.ndarray)

    # eigh_tridiagonal is only available from scipy 1.0.0
    try:
        from scipy.linalg import eigh_tridiagonal
    except ImportError:
        # Use eig_banded
        x, V = eig_banded(numpy.vstack((numpy.sqrt(beta), alpha)), lower=False)
        w = beta[0]*scipy.real(scipy.power(V[0, :], 2))
    else:
        x, V = eigh_tridiagonal(alpha, numpy.sqrt(beta[1:]))
        w = beta[0] * V[0, :]**2

    return x, w

def test_scipy(pyi_builder):
    pyi_builder.test_source(
        """
        from distutils.version import LooseVersion

        # Test top-level SciPy importability.
        import scipy
        from scipy import *

        # Test hooked SciPy modules.
        import scipy.io.matlab
        import scipy.sparse.csgraph

        # Test problematic SciPy modules.
        import scipy.linalg
        import scipy.signal

        # SciPy >= 0.16 privatized the previously public "scipy.lib" package as
        # "scipy._lib". Since this package is problematic, test its
        # importability regardless of SciPy version.
        if LooseVersion(scipy.__version__) >= LooseVersion('0.16.0'):
            import scipy._lib
        else:
            import scipy.lib
        """)

def MVG_DKL(M0,P0,M1,P1):
    '''
    KL divergence between two Gaussians

    Example
    -------
        M = randn(10)
        Q = randn(N,N)
        P = Q.dot(Q.T)
        MGV_DKL(M,P,M,P)
    '''
    MVG_check(M0,P0)
    MVG_check(M1,P1)
    N = len(M0)
    M1M0 = M1-M0
    return 0.5*(np.sum(P1*pinv(P0))+logdet(P0)-logdet(P1)-N+M1M0.T.dot(P1).dot(M1M0))
    #return 0.5*(np.sum(np.diag(P1.dot(np.linalg.pinv(P0))))+logdet(P0)-logdet(P1)-N+M1M0.T.dot(P1).dot(M1M0))

def gaussian1DblurOperator(n,sigma):
    '''
    Returns a 1D Gaussan blur operator of size n
    '''
    x   = np.linspace(0,n-1,n); # 1D domain
    tau = 1.0/sigma**2;       # precision
    k   = np.exp(-tau*x**2);    # compute (un-normalized) 1D kernel
    op  = scipy.linalg.special_matrices.toeplitz(k,k);     # convert to an operator from n -> n
    # normalize rows so density is conserved
    op /= np.sum(op)
    # truncate small entries
    big = np.max(op)
    toosmall = 1e-4*big
    op[op<toosmall] = 0
    # (re) normalize rows so density is conserved
    op /= np.sum(op)
    return op

def gaussian2DblurOperator(n,sigma):
    '''
    Returns a 2D Gaussan blur operator for a n x n sized domain
    Constructed as a tensor product of two 1d blurs of size n
    '''
    x   = np.linspace(0,n-1,n) # 1D domain
    tau = 1.0/sigma**2       # precision
    k   = np.exp(-tau*x**2)    # compute (un-normalized) 1D kernel
    tp  = scipy.linalg.special_matrices.toeplitz(k,k)     # convert to an operator from n -> n
    op  = scipy.linalg.special_matrices.kron(tp,tp)       # take the tensor product to get 2D operator
    # normalize rows so density is conserved
    op /= np.sum(op,axis=1)
    # truncate small entries
    big = np.max(op)
    toosmall = 1e-4*big
    op[op<toosmall] = 0
    # (re) normalize rows so density is conserved
    op /= np.sum(op,axis=1)
    return op

def _getAplus(A):
    '''
    Please see the documentation for nearPDHigham
    '''
    eigval, eigvec = np.linalg.eig(A)
    Q = np.matrix(eigvec)
    xdiag = np.matrix(np.diag(np.maximum(eigval, 0)))
    return Q*xdiag*Q.T

def measure_background(image, Fibers, width=30, niter=3, order=3):
    t = []
    a,b = image.shape
    ygrid,xgrid = np.indices(image.shape)
    ygrid = 1. * ygrid.ravel() / a
    xgrid = 1. * xgrid.ravel() / b
    image = image.ravel()
    s = np.arange(a*b)
    for fiber in Fibers:
        t.append(fiber.D*fiber.yind + fiber.xind)
    t = np.hstack(t)
    t = np.array(t, dtype=int)
    ind = np.setdiff1d(s,t)
    mask = np.zeros((a*b))
    mask[ind] = 1.
    mask[ind] = 1.-is_outlier(image[ind])
    sel = np.where(mask==1.)[0]
    for i in xrange(niter):
        V = polyvander2d(xgrid[sel],ygrid[sel],[order,order])
        sol = np.linalg.lstsq(V, image[sel])[0]
        vals = np.dot(V,sol) - image[sel]
        sel = sel[~is_outlier(vals)]
    V = polyvander2d(xgrid,ygrid,[order,order])
    back = np.dot(V, sol).reshape(a,b)    
    return back

def matrixMult(A, B):
    try:
        linalg.blas
    except AttributeError:
        return np.dot(A, B)

    if not A.flags['F_CONTIGUOUS']:
        AA = A.T
        transA = True
    else:
        AA = A
        transA = False

    if not B.flags['F_CONTIGUOUS']:
        BB = B.T
        transB = True
    else:
        BB = B
        transB = False

    return linalg.blas.dgemm(alpha=1., a=AA, b=BB, trans_a=transA, trans_b=transB)

def mat_inv(A):
#     I = numpy.eye(A.shape[0], A.shape[1])
    B = scipy.linalg.pinv2(A)
    return B

def _batch_incremental_pca(x, G, S):
    r = G.shape[1]
    b = x.shape[0]

    xh = G.T.dot(x)
    H  = x - G.dot(xh)
    J, W = scipy.linalg.qr(H, overwrite_a=True, mode='full', check_finite=False)

    Q = np.bmat( [[np.diag(S), xh], [np.zeros((b,r), dtype=np.float32), W]] )

    G_new, St_new, Vtoss = scipy.linalg.svd(Q, full_matrices=False, check_finite=False)
    St_new=St_new[:r]
    G_new= np.asarray(np.bmat([G, J]).dot( G_new[:,:r] ))

    return G_new, St_new

def test_recoverG(self):
    '''
    Test GCCA implementation by seeing if it can recover G.
    '''

    eps = 1.e-10

    Vs = [self.V1, self.V2, self.V3]

    wgcca = WGCCA.WeightedGCCA(3, [self.F1, self.F2, self.F3],
                               self.k, [eps, eps, eps], verbose=True)
    wgcca.learn(Vs)
    U1 = wgcca.U[0]
    U2 = wgcca.U[1]
    U3 = wgcca.U[2]

    Gprime   = wgcca.G

    # Rotate G to minimize norm of difference between G and G'
    R, B = scipy.linalg.orthogonal_procrustes(self.G, Gprime)
    normDiff = scipy.linalg.norm(self.G.dot(R) - Gprime)

    print ('Recovered G up to rotation; difference in norm:', normDiff)

    self.assertTrue( normDiff < 1.e-6 )
    self.assertTrue( np.allclose(self.G.dot(R), Gprime) )

def for_loop_update_row_param_blockwise(self, y_csr, phi_csr, mu0, c, v, r_prev, u_prev):

        nrow = y_csr.shape[0]
        num_factor = v.shape[1]
        prior_Phi = np.diag(np.hstack((self.prior_param['row_bias_scale'] ** -2,
                                       np.tile(self.prior_param['factor_scale'] ** -2, num_factor))))

        # Pre-allocate
        r = np.zeros(nrow)
        u = np.zeros((nrow, num_factor))

        # NOTE: The loop through 'i' is completely parallelizable.
        for i in range(nrow):
            j = y_csr[i, :].indices
            nnz_i = len(j)
            residual_i = y_csr[i, :].data - mu0 - c[j]
            phi_i = phi_csr[i, :].data.copy()

            v_T = np.hstack((np.ones((nnz_i, 1)), v[j, :]))
            post_Phi_i = prior_Phi + \
                         np.dot(v_T.T,
                                np.tile(phi_i[:, np.newaxis], (1, 1 + num_factor)) * v_T)  # Weighted sum of v_j * v_j.T
            post_mean_i = np.squeeze(np.dot(phi_i * residual_i, v_T))

            C, lower = scipy.linalg.cho_factor(post_Phi_i)
            post_mean_i = scipy.linalg.cho_solve((C, lower), post_mean_i)
            # Generate Gaussian, recycling the Cholesky factorization from the posterior mean computation.
            ru_i = math.sqrt(1 - self.relaxation ** 2) * scipy.linalg.solve_triangular(C, np.random.randn(len(post_mean_i)),
                                                                                       lower=lower)
            ru_i += post_mean_i + self.relaxation * (post_mean_i - np.concatenate(([r_prev[i]], u_prev[i, :])))
            r[i] = ru_i[0]
            u[i, :] = ru_i[1:]

        return r, u

def for_loop_update_col_param_blockwise(self, y_csc, phi_csc, mu0, r, u, c_prev, v_prev):

        ncol = y_csc.shape[1]
        num_factor = u.shape[1]
        prior_Phi = np.diag(np.hstack((self.prior_param['col_bias_scale'] ** -2,
                                       np.tile(self.prior_param['factor_scale'] ** -2, num_factor))))

        # Pre-allocate
        c = np.zeros(ncol)
        v = np.zeros((ncol, num_factor))

        # NOTE: The loop through 'j' is completely parallelizable.
        for j in range(ncol):
            i = y_csc[:, j].indices
            nnz_j = len(i)
            residual_j = y_csc[:, j].data - mu0 - r[i]
            phi_j = phi_csc[:, j].data

            u_T = np.hstack((np.ones((nnz_j, 1)), u[i, :]))
            post_Phi_j = prior_Phi + \
                         np.dot(u_T.T,
                                np.tile(phi_j[:, np.newaxis], (1, 1 + num_factor)) * u_T)  # Weighted sum of u_i * u_i.T
            post_mean_j = np.squeeze(np.dot(phi_j * residual_j, u_T))

            C, lower = scipy.linalg.cho_factor(post_Phi_j)
            post_mean_j = scipy.linalg.cho_solve((C, lower), post_mean_j)
            # Generate Gaussian, recycling the Cholesky factorization from the posterior mean computation.
            cv_j = math.sqrt(1 - self.relaxation ** 2) * scipy.linalg.solve_triangular(C, np.random.randn(len(post_mean_j)),
                                                                                       lower=lower)
            cv_j += post_mean_j + self.relaxation * (post_mean_j - np.concatenate(([c_prev[j]], v_prev[j, :])))
            c[j] = cv_j[0]
            v[j, :] = cv_j[1:]

        return c, v

def check_gradient(self, X, T, eps=1e-10):
        thetas = self.flatten()

        grad1 = self.numerical_gradient(thetas, X, T)
        _, grad2 = self.compute_cost(thetas, X, T)

        diff = linalg.norm(grad1 - grad2) / linalg.norm(grad1 + grad2)
        print(np.c_[grad1, grad2, np.abs(grad1 - grad2)])
        print('diff = {0}'.format(diff))

        return diff < eps

def ComputeCCA(X, Y):
  assert X.shape[0] == Y.shape[0], (X.shape, Y.shape, "Unequal number of rows")
  assert X.shape[0] > 1, (X.shape, "Must have more than 1 row")

  X = NormCenterMatrix(X)
  Y = NormCenterMatrix(Y)
  X_q, _, _ = decomp_qr.qr(X, overwrite_a=True, mode='economic', pivoting=True)
  Y_q, _, _ = decomp_qr.qr(Y, overwrite_a=True, mode='economic', pivoting=True)
  C = np.dot(X_q.T, Y_q)
  r = linalg.svd(C, full_matrices=False, compute_uv=False)
  d = min(X.shape[1], Y.shape[1])
  r = r[:d]
  r = np.minimum(np.maximum(r, 0.0), 1.0)  # remove roundoff errs
  return r.mean()

def compute_phi_prior(self):
        logZ_prior = 0
        for i in range(self.no_layers):
            Dout_i = self.layer_sizes[i+1]
            if not self.zu_tied:
                for d in range(Dout_i):
                    (sign, logdet) = np.linalg.slogdet(self.Kuu[i][d])
                    logZ_prior += 0.5 * logdet
            else:
                (sign, logdet) = np.linalg.slogdet(self.Kuu[i])
                logZ_prior += Dout_i * 0.5 * logdet

        return logZ_prior

def compute_phi_posterior(self):
        logZ_posterior = 0
        for i in range(self.no_layers):

            Dout_i = self.layer_sizes[i+1]
            for d in range(Dout_i):
                mud_val = self.mu[i][d]
                Sud_val = self.Su[i][d]
                (sign, logdet) = np.linalg.slogdet(Sud_val)
                # print 'phi_poste: ', 0.5 * logdet, 0.5 * np.sum(mud_val * spalg.solve(Sud_val, mud_val.T))
                logZ_posterior += 0.5 * logdet
                logZ_posterior += 0.5 * np.sum(mud_val * spalg.solve(Sud_val, mud_val.T))
        return logZ_posterior

def compute_phi_cavity(self):
        phi_cavity = 0
        for i in range(self.no_layers):
            Dout_i = self.layer_sizes[i+1]
            for d in range(Dout_i):
                muhatd_val = self.muhat[i][d]
                Suhatd_val = self.Suhat[i][d]
                (sign, logdet) = np.linalg.slogdet(Suhatd_val)
                phi_cavity += 0.5 * logdet
                phi_cavity += 0.5 * np.sum(muhatd_val * spalg.solve(Suhatd_val, muhatd_val.T))
        return phi_cavity

def CreateLmComp(self):
        '''
        This function generates components for the linear model: hrf, whitening matrix, autocorrelation matrix and CX
        '''

        # hrf
        self.canonical(0.1)

        # contrasts
        # expand contrasts to resolution of HRF (?)
        prec = int(self.laghrf/(self.hrf_precision/self.resolution))
        self.CX = np.kron(self.C, np.eye(prec))

        # drift
        self.S = self.drift(np.arange(0, self.n_scans))  # [tp x 1]
        self.S = np.matrix(self.S)

        # square of the whitening matrix
        base = [1 + self.rho**2, -1 * self.rho] + [0] * (self.n_scans - 2)
        self.V2 = scipy.linalg.toeplitz(base)
        self.V2[0, 0] = 1
        self.V2 = np.matrix(self.V2)
        self.V2[self.n_scans - 1, self.n_scans - 1] = 1

        self.white = self.V2 - self.V2 * \
            t(self.S) * np.linalg.pinv(self.S *
                                       self.V2 * t(self.S)) * self.S * self.V2

        return self

def find_frame_of_reference(target, source):
    target = copy(target)
    source = copy(source[:, :target.shape[1]])

    target /= target[3,:]
    source /= source[3,:]

    return dot(linalg.pinv(source.T), target.T).T