Python scipy.sparse 模块，triu() 实例源码

def test_cholmod(self):
        A, chol_L, _, cv = pickle.load(vc('/unittest/linalg/cholmod.pkl'))

        c_data = np.ones(len(cv))/len(cv)
        c_rows = cv.flatten()
        c_cols = (np.zeros(len(cv))).astype(np.int32)
        c = sp.csc_matrix((c_data, (c_rows, c_cols)), shape=(A.shape[0], 1))
        Ac = sp.hstack([A, c], format='csc')

        AAc = Ac.dot(Ac.T)

        [chol_L_comp, L_nonpsd, chol_S_comp] = lchol.lchol(AAc)

        right = chol_S_comp.T.dot(AAc.dot(chol_S_comp))
        left = chol_L_comp.dot(chol_L_comp.T)

        self.assertTrue(sum((abs(right-left)).data))  # it's a reordered LLt decomposition
        self.assertEqual(sp.triu(chol_L, k=1).nnz, 0) # it's lower triangular'
        self.assertEqual(L_nonpsd, 0)                 # the input is positive definite
        # self.assertTrue(sum((abs(chol_L - chol_L_comp)).data) < 1e-1)
        # self.assertTrue(sum((abs(chol_S - chol_S_comp)).data) < 1e-1)

def sqrtvc(m):
    mup=m
    mdown=mup.transpose()
    mdown.setdiag(0)
    mtogether=mup+mdown
    sums_sq=np.sqrt(mtogether.sum(axis=1)) 
    D_sq = sps.spdiags(1.0/sums_sq.flatten(), [0], mtogether.get_shape()[0], mtogether.get_shape()[1], format='csr')
    return sps.triu(D_sq.dot(mtogether.dot(D_sq)))

def hichip_add_diagonal(m):
    mup=m
    mdown=mup.transpose()
    mdown.setdiag(0)
    mtogether=mup+mdown
    sums=mtogether.sum(axis=1)
    max_sum=np.max(sums)
    to_add=1.0*max_sum-1.0*sums
    to_add_values=[]
    for i in range(m.shape[0]):
        to_add_values.append(to_add[i,0])
    mtogether.setdiag(np.array(to_add_values))
    D = sps.spdiags(1.0/sums.flatten(), [0], mtogether.get_shape()[0], mtogether.get_shape()[1], format='csr')
    return sps.triu(D.dot(mtogether))

def coverage_norm(m):
    mup=m
    mdown=mup.transpose()
    mdown.setdiag(0)
    mtogether=mup+mdown
    sums=mtogether.sum(axis=1)
    D = sps.spdiags(1.0/sums.flatten(), [0], mtogether.get_shape()[0], mtogether.get_shape()[1], format='csr')
    return sps.triu(D.dot(mtogether.dot(D)))

#assumes matrix is upper triangular

def array_2_coverageVector(m):
    assert np.allclose(m, np.triu(m))
    m_sym=m+m.T-m.diagonal()
    return m_sym.sum(axis=0)

def subsample_to_depth_array_upperTri(m,seq_depth):
    m=np.triu(m)
    subsampled_data=np.zeros(m.shape)
    depthm=m.sum()
    assert seq_depth<=depthm
    subsampling_prob=seq_depth/depthm
    for i in range(m.shape[0]):
        for j in range(m.shape[1]):
            if j<=i:
                continue
            n=m[i,j]
            subsampled_data[i,j]=np.random.binomial(n,subsampling_prob,1)[0]
    return subsampled_data

def omega_index(clustering1, clustering2):
    '''
    Omega index

    (cite)
    '''

    A1 = make_overlapping_membership_matrix(clustering1)
    A2 = make_overlapping_membership_matrix(clustering2)

    M = clustering1.n_elements * (clustering1.n_elements - 1) / 2.0

    maxNover = max(max(A1.diagonal()), max(A2.diagonal())) + 1

    Anot = spsparse.triu((A1 != A2), k = 1).sum()

    omega_u = 1.0 - Anot.sum() / M

    t_0_1 = M - spsparse.triu((A1 != 0), k = 1).sum()
    t_0_2 = M - spsparse.triu((A2 != 0), k = 1).sum()

    t_k_1 = [spsparse.triu((A1 == i), k = 1).sum() for i in xrange(1, maxNover)]
    t_k_2 = [spsparse.triu((A2 == i), k = 1).sum() for i in xrange(1, maxNover)]

    omega_e = (t_0_1*t_0_2 + np.dot(t_k_1, t_k_2) ) / M**2

    return (omega_u - omega_e) / (1.0 - omega_e)

def get_upper_triangle(X, k=1, is_sparse=True):
    if is_sparse:
        return sparse.triu(X, k=k).toarray().flatten()
    else:
        return np.triu(X, k=k).flatten()

def _triu(a, sparse):
    if sparse:
        return sp.triu(a, k=1)
    return np.triu(a, k=1)