Python MatrixUtil.double_centered方法代码示例

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def get_response_content(fs):
    D = fs.matrix
    L = Euclid.edm_to_laplacian(D)
    S = get_sigma_matrix(D)
    P = get_precision_matrix(S)
    # begin the response
    out = StringIO()
    print >> out, 'the Laplacian matrix:'
    print >> out, MatrixUtil.m_to_string(L)
    print >> out
    print >> out, 'the sigma matrix corresponding to the Q matrix:'
    print >> out, MatrixUtil.m_to_string(S)
    print >> out
    print >> out, 'the precision matrix corresponding to the Q matrix:'
    print >> out, MatrixUtil.m_to_string(P)
    print >> out
    print >> out, 'the precision matrix minus the laplacian matrix:'
    print >> out, MatrixUtil.m_to_string(P-L)
    print >> out
    print >> out, 'the double centered precision matrix minus the laplacian matrix:'
    print >> out, MatrixUtil.m_to_string(MatrixUtil.double_centered(P)-L)
    print >> out
    print >> out, 'the pseudo-inverse of the double centered sigma matrix minus the laplacian matrix:'
    print >> out, MatrixUtil.m_to_string(np.linalg.pinv(MatrixUtil.double_centered(S))-L)
    # write the response
    return out.getvalue()

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def get_response_content(fs):
    # read the distance matrix
    D = fs.matrix
    # if the distances are plane-like then square each element of the distance matrix
    if fs.planelike:
        D = D**2
    # get the A matrix
    A = -0.5 * D
    # get the doubly centered A matrix
    HAH = MatrixUtil.double_centered(A)
    # do the eigendecomposition
    eigenvalues, eigenvector_transposes = np.linalg.eigh(HAH)
    eigenvectors = eigenvector_transposes.T
    eigensystem = [(abs(w), w, v.tolist()) for w, v in zip(eigenvalues, eigenvectors)]
    sorted_eigensystem = list(reversed(sorted(eigensystem)))
    sorted_abs_eigenvalues, sorted_eigenvalues, sorted_eigenvectors = zip(*sorted_eigensystem)
    # get the points that approximate the distance matrix
    top_eigenvalues = sorted_eigenvalues[:2]
    top_eigenvectors = sorted_eigenvectors[:2]
    for eigenvalue in top_eigenvalues:
        if eigenvalue <= 0:
            msg = 'one of the top two eigenvalues was non-positive'
            raise HandlingError(msg)
    axes = []
    for eigenvalue, eigenvector in zip(top_eigenvalues, top_eigenvectors):
        axis = [v * math.sqrt(eigenvalue) for v in eigenvector]
        axes.append(axis)
    points = zip(*axes)
    # begin the response
    out = StringIO()
    for point in points:
        print >> out, '\t'.join(str(v) for v in point)
    # return the response
    return out.getvalue()

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def get_eigendecomposition_report(D):
    """
    @param D: a distance matrix
    @return: a multi-line string
    """
    out = StringIO()
    # get some intermediate matrices and vectors
    L = Euclid.edm_to_laplacian(D)
    laplacian_fiedler = BuildTreeTopology.laplacian_to_fiedler(L)
    distance_fiedler = BuildTreeTopology.edm_to_fiedler(D)
    eigensplit = BuildTreeTopology.eigenvector_to_split(laplacian_fiedler)
    # report the two eigenvalue lists that should be the same
    HDH = MatrixUtil.double_centered(D)
    HSH = -0.5 * HDH
    w_distance, vt_distance = np.linalg.eigh(HSH)
    print >> out, 'the laplacian-derived and distance-derived eigenvalues:'
    w_laplacian, vt_laplacian = np.linalg.eigh(L)
    for a, b in zip(sorted(w_laplacian), sorted(w_distance)):
        print >> out, a, '\t', b
    print >> out
    # report the two fiedler vectors that should be the same
    print >> out, 'the laplacian-derived and distance-derived fiedler vectors:'
    for a, b in zip(laplacian_fiedler, distance_fiedler):
        print >> out, a, '\t', b
    return out.getvalue().strip()

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def get_augmented_gower_selection(D):
    """
    Do a spectral sign split with neighbor joining fallback.
    The first choice is to return indices corresponding to
    positive elements of the dominant eigenvector of the gower matrix.
    If this defines a degenerate bipartition,
    then neighbor joining is used as a fallback.
    @param D: a distance matrix
    @return: the set of selected indices
    """
    n = len(D)
    if n < 4:
        raise ValueError('expected a distance matrix with at least four rows')
    # get the gower matrix
    G = MatrixUtil.double_centered(numpy.array(D))
    # get the dominant eigenvector
    eigenvalues, eigenvector_transposes = linalg.eigh(G)
    eigenvectors = eigenvector_transposes.T
    dominant_value, dominant_vector = max((abs(w), v) for w, v in zip(eigenvalues, eigenvectors))
    # get the bipartition defined by the dominant eigenvector
    selection = set(i for i, x in enumerate(dominant_vector) if x > 0)
    complement = set(range(n)) - selection
    # if the bipartition is degenerate then resort to neighbor joining
    if min(len(selection), len(complement)) < 2:
        selection = set(NeighborJoining.get_neighbors(D))
    return selection

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def get_stability(D):
    """
    The stability is defined as a bound on norms of perturbation matrices.
    If D is perturbed by a matrix whose Frobenius norm
    is less than the stability, then the spectral split remains unchanged.
    @param D: a distance matrix
    @return: the stability of the distance matrix
    """
    HDH = MatrixUtil.double_centered(D)
    # get the eigendecomposition
    w_v_pairs = get_sorted_eigensystem(-HDH)
    # compute the eigengap
    w = [w for w, v in w_v_pairs]
    lambda_1 = w[-1]
    lambda_2 = w[-2]
    eigengap = lambda_1 - lambda_2
    delta = eigengap
    # compute an eigenvector stability statistic
    v = [v for w, v in w_v_pairs]
    dominant_eigenvector = v[-1]
    alpha = min(abs(x) for x in dominant_eigenvector)
    # compute the stability as a function of alpha and delta
    eigenvalue_control = delta / (2*math.sqrt(2))
    eigenvector_control = alpha * delta / (4 + math.sqrt(2)*alpha) 
    stability = min(eigenvalue_control, eigenvector_control)
    return stability

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def edm_to_dccov(D):
    """
    @param D: a Euclidean distance matrix
    @return: a double centered covariance matrix
    """
    MatrixUtil.assert_square(D)
    return -(0.5)*MatrixUtil.double_centered(D)

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def get_response_content(fs):
    np.set_printoptions(linewidth=200)
    n = len(fs.D)
    # create the distance matrix with the extra node added
    D_dup = get_duplicate_edm(fs.D)
    # get the principal axis projection from the distance dup matrix
    HSH = -(0.5) * MatrixUtil.double_centered(D_dup)
    X_w, X_v = EigUtil.principal_eigh(HSH)
    D_dup_x = X_v * math.sqrt(X_w)
    # get masses summing to one
    m = np.array([1]*(n-1) + [2], dtype=float) / (n+1)
    # get the principal axis projection using the weight formula
    M = np.diag(np.sqrt(m))
    I = np.eye(n, dtype=float)
    E = I - np.outer(np.ones(n, dtype=float), m)
    ME = np.dot(M, E)
    Q = -(0.5) * np.dot(ME, np.dot(fs.D, ME.T))
    Q_w, Q_v = EigUtil.principal_eigh(Q)
    Q_x = Q_v * math.sqrt(Q_w) / np.sqrt(m)
    # make the response
    out = StringIO()
    print >> out, 'distance matrix with exact duplicate node:'
    print >> out, D_dup
    print >> out
    print >> out, 'principal axis projection:'
    print >> out, D_dup_x
    print >> out
    print >> out, 'principal axis projection using the weight formula:'
    print >> out, Q_x
    return out.getvalue()

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def get_principal_coordinate(D):
    """
    Return the principal eigenvector of the pseudoinverse of the laplacian.
    @param D: a distance matrix
    @return: the principal coordinate
    """
    L_pinv = -0.5 * MatrixUtil.double_centered(D)
    return get_principal_eigenvector(L_pinv)

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
 def test_distance_to_schur(self):
     leaves = T_to_leaves(g_example_T)
     # Compute the Schur complement Laplacian and the leaf distance matrix.
     L_schur = TB_to_L_schur(g_example_T, g_example_B, leaves)
     Dpp_direct = TB_to_D(g_example_T, g_example_B, leaves)
     # Compute one from the other.
     HDppH = MatrixUtil.double_centered(Dpp_direct)
     L_schur_estimate = np.linalg.pinv(-0.5*HDppH)
     self.assertTrue(np.allclose(L_schur_estimate, L_schur))

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
 def test_degenerate_lfdn(self):
     """
     Make sure that replication with a factor of 1 does not do anything.
     """
     lfdo = tree_string_to_LFDO(g_tree_string)
     N = 1
     lfdn = LFDO_to_LFDN(lfdo, N)
     HDH = MatrixUtil.double_centered(lfdo.M)
     self.assertTrue(np.allclose(lfdn.M, HDH))

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def get_aug_b(D, p, q, N):
    """
    Get the -(1/2)HMDM'H augmented matrix.
    @param D: the full distance matrix
    @param p: the number of internal nodes
    @param q: the number of tips
    @param N: the tip duplication factor
    @return: an augmented matrix for comparison
    """
    M = get_M(p, q, N)
    D_aug = np.dot(M, np.dot(D, M.T))
    return -0.5 * MatrixUtil.double_centered(D_aug)

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
 def test_lfdi_submatrix(self):
     """
     Test a principal submatrix of the LFDI matrix.
     The principal submatrix of the LFDI which corresponds to leaf vertices
     should equal the centered leaf vertex distance matrix.
     """
     lfdo = tree_string_to_LFDO(g_tree_string)
     q = lfdo.q
     DQ = lfdo.M[:q, :q]
     HDQH = MatrixUtil.double_centered(DQ)
     lfdi = LFDO_to_LFDI(lfdo)
     self.assertTrue(np.allclose(lfdi.M[:q, :q], HDQH))

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def TB_to_G(T, B, vertices):
    """
    Get the Gower matrix.
    Note that this should be the pseudoinverse of L_schur.
    @param T: topology
    @param B: branch lengths
    @param vertices: ordered vertices
    @return: the Gower matrix
    """
    D = TB_to_D(T, B, vertices)
    HDH = MatrixUtil.double_centered(D)
    return -0.5 * HDH

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def get_split(D):
    """
    @param D: an exact or perturbed distance matrix
    @return: a set of frozensets of indices
    """
    HDH = MatrixUtil.double_centered(D)
    # get the dominant eigenvector
    w_v_pairs = get_sorted_eigensystem(-HDH)
    v = [v for w, v in w_v_pairs]
    ev_dom = v[-1]
    neg_set = frozenset(i for i, x in enumerate(ev_dom) if x < 0)
    nonneg_set = frozenset(i for i, x in enumerate(ev_dom) if x >= 0)
    return set([neg_set, nonneg_set])

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import double_centered [as 别名]
def get_augmented_spectrum(D_in, ntips, ndup_tip, ndup_internal):
    """
    The tips are the first indices of the original distance matrix.
    @param D_in: the original distance matrix
    @param ntips: the number of tips in the tree
    @param ndup_tip: the total number of repeats per tip node
    @param ndup_internal: the total number of repeats per internal node
    @return: eigenvalues of Gower's centered augmented distance matrix
    """
    vdup = [ndup_tip]*ntips + [ndup_internal]*(len(D_in) - ntips)
    D = get_dup_distance_matrix(D_in, vdup)
    G = -0.5 * MatrixUtil.double_centered(D)
    w = scipy.linalg.eigh(G, eigvals_only=True)
    return (w * ndup_internal) / ndup_tip