Python scipy.eye() Examples

def test_complex_lu(self):
        """Getting factors of complex matrix"""
        umfpack = um.UmfpackContext("zi")

        for A in self.complex_matrices:
            umfpack.numeric(A)

            (L,U,P,Q,R,do_recip) = umfpack.lu(A)

            L = L.todense()
            U = U.todense()
            A = A.todense()
            if not do_recip:
                R = 1.0/R
            R = matrix(diag(R))
            P = eye(A.shape[0])[P,:]
            Q = eye(A.shape[1])[:,Q]

            assert_array_almost_equal(P*R*A*Q,L*U) 

def test_real_lu(self):
        """Getting factors of real matrix"""
        umfpack = um.UmfpackContext("di")

        for A in self.real_matrices:
            umfpack.numeric(A)

            (L,U,P,Q,R,do_recip) = umfpack.lu(A)

            L = L.todense()
            U = U.todense()
            A = A.todense()
            if not do_recip:
                R = 1.0/R
            R = matrix(diag(R))
            P = eye(A.shape[0])[P,:]
            Q = eye(A.shape[1])[:,Q]

            assert_array_almost_equal(P*R*A*Q,L*U) 

def test_complex_lu(self):
        # Getting factors of complex matrix
        umfpack = um.UmfpackContext("zi")

        for A in self.complex_matrices:
            umfpack.numeric(A)

            (L,U,P,Q,R,do_recip) = umfpack.lu(A)

            L = L.todense()
            U = U.todense()
            A = A.todense()
            if not do_recip:
                R = 1.0/R
            R = matrix(diag(R))
            P = eye(A.shape[0])[P,:]
            Q = eye(A.shape[1])[:,Q]

            assert_array_almost_equal(P*R*A*Q,L*U) 

def test_complex_int64_lu(self):
        # Getting factors of complex matrix with long indices
        umfpack = um.UmfpackContext("zl")

        for A in self.complex_int64_matrices:
            umfpack.numeric(A)

            (L,U,P,Q,R,do_recip) = umfpack.lu(A)

            L = L.todense()
            U = U.todense()
            A = A.todense()
            if not do_recip:
                R = 1.0/R
            R = matrix(diag(R))
            P = eye(A.shape[0])[P,:]
            Q = eye(A.shape[1])[:,Q]

            assert_array_almost_equal(P*R*A*Q,L*U) 

def test_real_lu(self):
        # Getting factors of real matrix
        umfpack = um.UmfpackContext("di")

        for A in self.real_matrices:
            umfpack.numeric(A)

            (L,U,P,Q,R,do_recip) = umfpack.lu(A)

            L = L.todense()
            U = U.todense()
            A = A.todense()
            if not do_recip:
                R = 1.0/R
            R = matrix(diag(R))
            P = eye(A.shape[0])[P,:]
            Q = eye(A.shape[1])[:,Q]

            assert_array_almost_equal(P*R*A*Q,L*U) 

def test_real_int64_lu(self):
        # Getting factors of real matrix with long indices
        umfpack = um.UmfpackContext("dl")

        for A in self.real_int64_matrices:
            umfpack.numeric(A)

            (L,U,P,Q,R,do_recip) = umfpack.lu(A)

            L = L.todense()
            U = U.todense()
            A = A.todense()
            if not do_recip:
                R = 1.0/R
            R = matrix(diag(R))
            P = eye(A.shape[0])[P,:]
            Q = eye(A.shape[1])[:,Q]

            assert_array_almost_equal(P*R*A*Q,L*U) 

def test_trivial():
    n = 5
    X = ones((n, 1))
    A = eye(n)
    compare_solutions(A, None, n) 

def test_trivial():
    n = 5
    X = ones((n, 1))
    A = eye(n)
    compare_solutions(A, None, n) 

def test_diagonal():
    # This test was moved from '__main__' in lobpcg.py.
    # Coincidentally or not, this is the same eigensystem
    # required to reproduce arpack bug
    # http://forge.scilab.org/index.php/p/arpack-ng/issues/1397/
    # even using the same n=100.

    np.random.seed(1234)

    # The system of interest is of size n x n.
    n = 100

    # We care about only m eigenpairs.
    m = 4

    # Define the generalized eigenvalue problem Av = cBv
    # where (c, v) is a generalized eigenpair,
    # and where we choose A to be the diagonal matrix whose entries are 1..n
    # and where B is chosen to be the identity matrix.
    vals = np.arange(1, n+1, dtype=float)
    A = scipy.sparse.diags([vals], [0], (n, n))
    B = scipy.sparse.eye(n)

    # Let the preconditioner M be the inverse of A.
    M = scipy.sparse.diags([np.reciprocal(vals)], [0], (n, n))

    # Pick random initial vectors.
    X = np.random.rand(n, m)

    # Require that the returned eigenvectors be in the orthogonal complement
    # of the first few standard basis vectors.
    m_excluded = 3
    Y = np.eye(n, m_excluded)

    eigs, vecs = lobpcg(A, X, B, M=M, Y=Y, tol=1e-4, maxiter=40, largest=False)

    assert_allclose(eigs, np.arange(1+m_excluded, 1+m_excluded+m))
    _check_eigen(A, eigs, vecs, rtol=1e-3, atol=1e-3) 

def test_hermitian():
    np.random.seed(1234)

    sizes = [3, 10, 50]
    ks = [1, 3, 10, 50]
    gens = [True, False]

    for size, k, gen in itertools.product(sizes, ks, gens):
        if k > size:
            continue

        H = np.random.rand(size, size) + 1.j * np.random.rand(size, size)
        H = 10 * np.eye(size) + H + H.T.conj()

        X = np.random.rand(size, k)

        if not gen:
            B = np.eye(size)
            w, v = lobpcg(H, X, maxiter=5000)
            w0, v0 = eigh(H)
        else:
            B = np.random.rand(size, size) + 1.j * np.random.rand(size, size)
            B = 10 * np.eye(size) + B.dot(B.T.conj())
            w, v = lobpcg(H, X, B, maxiter=5000)
            w0, v0 = eigh(H, B)

        for wx, vx in zip(w, v.T):
            # Check eigenvector
            assert_allclose(np.linalg.norm(H.dot(vx) - B.dot(vx) * wx) / np.linalg.norm(H.dot(vx)),
                            0, atol=5e-4, rtol=0)

            # Compare eigenvalues
            j = np.argmin(abs(w0 - wx))
            assert_allclose(wx, w0[j], rtol=1e-4) 

def __MR_affinity_matrix(self,img,labels):        
        W,D = self.__MR_W_D_matrix(img,labels)
        aff = pinv(D-self.weight_parameters['alpha']*W)
        aff[sp.eye(sp.amax(labels)+1).astype(bool)] = 0.0 # diagonal elements to 0
        return aff 

def build_from_adjacency_matrix(name,
                                  adj,
                                  features,
                                  train_mask,
                                  val_mask,
                                  test_mask,
                                  labels,
                                  row_normalize=False):
    """Build from adjacency matrix."""
    # Extract train, val, test, unlabeled indices.
    train_indices = np.where(train_mask)[0]
    test_indices = np.where(test_mask)[0]
    val_indices = np.where(val_mask)[0]
    unlabeled_mask = np.logical_not(train_mask | test_mask | val_mask)
    unlabeled_indices = np.where(unlabeled_mask)[0]

    # Extract node features.
    if row_normalize:
      features = GCNDataset.preprocess_features(features)

    features = np.float32(features.todense())

    # Extract labels.
    labels = np.argmax(labels, axis=-1)
    num_classes = max(labels) + 1

    # Extract edges.
    adj = scipy.sparse.coo_matrix(adj)
    edges = [
        GCNDataset.Edge(src, tgt, val)
        for src, tgt, val in zip(adj.row, adj.col, adj.data)
    ]

    # Preprocessing of adjacency matrix for simple GCN model and conversion to
    # tuple representation.
    adj_normalized = GCNDataset.normalize_adj(adj +
                                              scipy.eye(adj.shape[0])).astype(
                                                  np.float32)
    support = GCNDataset.sparse_to_tuple(adj_normalized)

    features_matrix = (
        GCNDataset.row_normalize(features).astype(np.float32)
        if row_normalize else features.astype(np.float32))
    features_sparse = GCNDataset.sparse_to_tuple(features_matrix)
    num_features_nonzero = features_sparse[1].shape
    return GCNDataset(
        name=name,
        features=features_matrix,
        support=support,
        num_features_nonzero=num_features_nonzero,
        features_sparse=features_sparse,
        labels=labels,
        edges=edges,
        indices_train=train_indices,
        indices_test=test_indices,
        indices_val=val_indices,
        indices_unlabeled=unlabeled_indices,
        num_classes=num_classes,
        feature_preproc_fn=lambda x: x) 

def ldpred_inf(beta_hats, h2=0.1, n=1000, inf_shrink_matrices=None, 
               reference_ld_mats=None, genotypes=None, ld_window_size=100, verbose=False):
    """
    Apply the infinitesimal shrink w LD (which requires LD information).
    
    If reference_ld_mats are supplied, it uses those, otherwise it uses the LD in the genotype data.
    
    If genotypes are supplied, then it assumes that beta_hats and the genotypes are synchronized.

    """
    n = float(n)
    if verbose:
        print('Doing LD correction')
    t0 = time.time()
    m = len(beta_hats)
    updated_betas = sp.empty(m)

    for i, wi in enumerate(range(0, m, ld_window_size)):
        start_i = wi
        stop_i = min(m, wi + ld_window_size)
        curr_window_size = stop_i - start_i
        if inf_shrink_matrices!=None:
            A_inv = inf_shrink_matrices[i]
        else:
            if reference_ld_mats != None:
                D = reference_ld_mats[i]
            else:
                if genotypes != None:
                    X = genotypes[start_i: stop_i]
                    num_indivs = X.shape[1]
                    D = sp.dot(X, X.T) / num_indivs
                else:
                    raise NotImplementedError
            A = ((m / h2) * sp.eye(curr_window_size) + (n / (1.0)) * D)
            A_inv = linalg.pinv(A)
        updated_betas[start_i: stop_i] = sp.dot(A_inv * n , beta_hats[start_i: stop_i])  # Adjust the beta_hats

        if verbose:
            sys.stdout.write('\r%0.2f%%' % (100.0 * (min(1, float(wi + ld_window_size) / m))))
            sys.stdout.flush()

    t1 = time.time()
    t = (t1 - t0)
    if verbose:
        print('\nIt took %d minutes and %0.2f seconds to perform the Infinitesimal LD shrink' % (t / 60, t % 60))
    return updated_betas