Python Euclid类代码示例

def get_response_content(fs):
    # check input compatibility
    if fs.nvertices < fs.naxes+1:
        msg_a = 'attempting to plot too many eigenvectors '
        msg_b = 'for the given number of vertices'
        raise ValueError(msg_a + msg_b)
    # define the requested physical size of the images (in pixels)
    physical_size = (640, 480)
    # get the points
    L = create_laplacian_matrix(fs.nvertices)
    D = Euclid.laplacian_to_edm(L)
    HSH = Euclid.edm_to_dccov(D)
    W, VT = np.linalg.eigh(HSH)
    V = VT.T.tolist()
    if fs.eigenvalue_scaling:
        vectors = [np.array(v)*w for w, v in list(reversed(sorted(zip(np.sqrt(W), V))))[:-1]]
    else:
        vectors = [np.array(v) for w, v in list(reversed(sorted(zip(np.sqrt(W), V))))[:-1]]
    X = np.array(zip(*vectors))
    # transform the points to eigenfunctions such that the first point is positive
    F = X.T[:fs.naxes]
    for i in range(fs.naxes):
        if F[i][0] < 0:
            F[i] *= -1
    # draw the image
    try:
        ext = Form.g_imageformat_to_ext[fs.imageformat]
        return create_image_string(ext, physical_size, F, fs.xaxis_length)
    except CairoUtil.CairoUtilError as e:
        raise HandlingError(e)

def process(tree_string):
    """
    @param tree_string: a newick string
    @return: a multi-line string that summarizes the results
    """
    np.set_printoptions(linewidth=200)
    out = StringIO()
    # build the newick tree from the string
    tree = NewickIO.parse(tree_string, FelTree.NewickTree)
    # get ordered names and ids
    ordered_ids, ordered_names = get_ordered_ids_and_names(tree)
    # get the distance matrix with ordered indices including all nodes in the tree
    nvertices = len(list(tree.preorder()))
    nleaves = len(list(tree.gen_tips()))
    id_to_index = dict((myid, i) for i, myid in enumerate(ordered_ids))
    D = np.array(tree.get_partial_distance_matrix(ordered_ids))
    # define mass vectors
    m_uniform_unscaled = [1]*nvertices
    m_degenerate_unscaled = [1]*nleaves + [0]*(nvertices-nleaves)
    m_uniform = np.array(m_uniform_unscaled, dtype=float) / sum(m_uniform_unscaled)
    m_degenerate = np.array(m_degenerate_unscaled, dtype=float) / sum(m_degenerate_unscaled)
    # show some of the distance matrices
    print >> out, 'ordered names:'
    print >> out, ordered_names
    print >> out
    print >> out, 'embedded points with mass uniformly distributed among all vertices:'
    print >> out, Euclid.edm_to_weighted_points(D, m_uniform)
    print >> out
    print >> out, 'embedded points with mass uniformly distributed among the leaves:'
    print >> out, Euclid.edm_to_weighted_points(D, m_degenerate)
    print >> out
    # return the response
    return out.getvalue().strip()

def process():
    """
    @return: a multi-line string that summarizes the results
    """
    np.set_printoptions(linewidth=200)
    out = StringIO()
    # define a degenerate mass vector
    m_degenerate = np.array([0.25, 0.25, 0.25, 0.25, 0, 0])
    # define some distance matrices
    D_leaves = Euclid.g_D_b
    D_all = Euclid.g_D_c
    nvertices = 6
    nleaves = 4
    # get the projection and the weighted multidimensional scaling
    X = Euclid.edm_to_points(D_all)
    Y = Euclid.edm_to_weighted_points(D_all, m_degenerate)
    D_X = np.array([[np.dot(pb-pa, pb-pa) for pa in X] for pb in X])
    D_Y = np.array([[np.dot(pb-pa, pb-pa) for pa in Y] for pb in Y])
    # get the embedding using only the leaves
    print >> out, 'embedding of leaves from the leaf distance matrix:'
    print >> out, Euclid.edm_to_points(D_leaves)
    print >> out, 'projection of all vertices onto the MDS space of the leaves:'
    print >> out, do_projection(D_all, nleaves)
    print >> out, 'embedding of all vertices using uniform weights:'
    print >> out, X
    print >> out, 'corresponding distance matrix:'
    print >> out, D_X
    print >> out, 'embedding of all vertices using degenerate weights:'
    print >> out, Y
    print >> out, 'corresponding distance matrix:'
    print >> out, D_Y
    return out.getvalue().strip()

def get_response_content(fs):
    # build the newick tree from the string
    tree = NewickIO.parse(fs.tree_string, FelTree.NewickTree)
    nvertices = len(list(tree.preorder()))
    nleaves = len(list(tree.gen_tips()))
    ninternal = nvertices - nleaves
    # get ordered ids with the internal nodes first
    ordered_ids = get_ordered_ids(tree)
    leaf_ids = [id(node) for node in tree.gen_tips()]
    # get the distance matrix and the augmented distance matrix
    D_leaf = np.array(tree.get_partial_distance_matrix(leaf_ids))
    D = np.array(tree.get_partial_distance_matrix(ordered_ids))
    D_aug = get_augmented_distance(D, nleaves, fs.ndups)
    # analyze the leaf distance matrix
    X_leaf = Euclid.edm_to_points(D_leaf)
    # get the eigendecomposition of the centered augmented distance matrix
    X_aug = Euclid.edm_to_points(D_aug, nvertices-1)
    # explicitly compute the points for the given number of dups using weights
    m = [1]*ninternal + [1+fs.ndups]*nleaves
    m = np.array(m, dtype=float) / sum(m)
    X_weighted = Euclid.edm_to_weighted_points(D, m)
    # explicitly compute the points for 10x dups
    m = [1]*ninternal + [1+fs.ndups*10]*nleaves
    m = np.array(m, dtype=float) / sum(m)
    X_weighted_10x = Euclid.edm_to_weighted_points(D, m)
    # explicitly compute the limiting points as the number of dups increases
    X = Euclid.edm_to_points(D)
    X -= np.mean(X[-nleaves:], axis=0)
    XL = X[-nleaves:]
    U, s, Vt = np.linalg.svd(XL)
    Z = np.dot(X, Vt.T)
    # report the results
    np.set_printoptions(linewidth=300, threshold=10000)
    out = StringIO()
    print >> out, 'leaf distance matrix:'
    print >> out, D_leaf
    print >> out
    print >> out, 'points derived from the leaf distance matrix'
    print >> out, '(the first column is proportional to the Fiedler vector):'
    print >> out, X_leaf
    print >> out
    if fs.show_aug:
        print >> out, 'augmented distance matrix:'
        print >> out, D_aug
        print >> out
    print >> out, 'points derived from the augmented distance matrix'
    print >> out, '(the first column is proportional to the Fiedler vector):'
    print >> out, get_ugly_matrix(X_aug, ninternal, nleaves)
    print >> out
    print >> out, 'points computed using masses:'
    print >> out, X_weighted
    print >> out
    print >> out, 'points computed using masses with 10x dups:'
    print >> out, X_weighted_10x
    print >> out
    print >> out, 'limiting points:'
    print >> out, Z
    print >> out
    return out.getvalue()

def get_response_content(fs):
    locations = get_locations()
    np_locs = [np.array(p) for p in locations]
    edges = get_edges()
    npoints = len(locations)
    # start writing the response
    np.set_printoptions(linewidth=200)
    out = StringIO()
    # print the layout data
    print >> out, 'POINTS'
    for i, (x, y) in enumerate(locations):
        print >> out, i, x, y
    print >> out, 'EDGES'
    for i, j in edges:
        print >> out, i, j
    print >> out
    # show the unweighted adjacency matrix
    UA = np.zeros((npoints, npoints))
    for i, j in edges:
        UA[i, j] = 1
        UA[j, i] = 1
    print >> out, 'unweighted adjacency matrix:'
    print >> out, UA
    print >> out
    # show the unweighted laplacian matrix
    UL = Euclid.adjacency_to_laplacian(UA)
    print >> out, 'unweighted laplacian matrix:'
    print >> out, UL
    print >> out
    # show the weighted adjacency matrix
    WA = np.zeros((npoints, npoints))
    for i, j in edges:
        d = np.linalg.norm(np_locs[i] - np_locs[j]) / math.sqrt(2.0)
        w = 1.0 / d
        WA[i, j] = w
        WA[j, i] = w
    print >> out, 'weighted adjacency matrix:'
    print >> out, WA
    print >> out
    # show the weighted laplacian matrix
    WL = Euclid.adjacency_to_laplacian(WA)
    print >> out, 'weighted laplacian matrix:'
    print >> out, WL
    print >> out
    # remove the two internal nodes by schur complementation
    ntips = 4
    schur_L = SchurAlgebra.schur_helper(WL, 2)
    X = Euclid.dccov_to_points(np.linalg.pinv(schur_L))
    print >> out, 'schur graph layout:'
    print >> out, 'POINTS'
    for i, v in enumerate(X):
        print >> out, i, v[0], v[1]
    print >> out, 'EDGES'
    for i in range(ntips):
        for j in range(i+1, ntips):
            print >> out, i, j
    # return the response
    return out.getvalue()

def get_response_content(fs):
    # get the tree
    tree = NewickIO.parse(fs.tree_string, FelTree.NewickTree)
    # get information about the tree topology
    internal = [id(node) for node in tree.gen_internal_nodes()]
    tips = [id(node) for node in tree.gen_tips()]
    vertices = internal + tips
    ntips = len(tips)
    ninternal = len(internal)
    nvertices = len(vertices)
    # get the ordered ids with the leaves first
    ordered_ids = vertices
    # get the full weighted adjacency matrix
    A = np.array(tree.get_affinity_matrix(ordered_ids))
    # compute the weighted adjacency matrix of the decorated tree
    p = ninternal
    q = ntips
    N = fs.N
    if fs.weight_n:
        weight = float(N)
    elif fs.weight_sqrt_n:
        weight = math.sqrt(N)
    A_aug = get_A_aug(A, weight, p, q, N)
    # compute the weighted Laplacian matrix of the decorated tree
    L_aug = Euclid.adjacency_to_laplacian(A_aug)
    # compute the eigendecomposition
    w, vt = np.linalg.eigh(L_aug)
    # show the output
    np.set_printoptions(linewidth=1000, threshold=10000)
    out = StringIO()
    if fs.lap:
        print >> out, 'Laplacian of the decorated tree:'
        print >> out, L_aug
        print >> out
    if fs.eigvals:
        print >> out, 'eigenvalues:'
        for x in w:
            print >> out, x
        print >> out
    if fs.eigvecs:
        print >> out, 'eigenvector matrix:'
        print >> out, vt
        print >> out
    if fs.compare:
        # get the distance matrix for only the original tips
        D_tips = np.array(tree.get_partial_distance_matrix(tips))
        X_tips = Euclid.edm_to_points(D_tips)
        # wring the approximate points out of the augmented tree
        X_approx = vt[p:p+q].T[1:1+q-1].T / np.sqrt(w[1:1+q-1])
        # do the comparison
        print >> out, 'points from tip-only MDS:'
        print >> out, X_tips
        print >> out
        print >> out, 'approximate points from decorated tree:'
        print >> out, X_approx
        print >> out
    return out.getvalue()

def process():
    """
    @return: a multi-line string that summarizes the results
    """
    np.set_printoptions(linewidth=200)
    # define the adjacency matrix
    A = g_A
    n = 6
    # define some mass distributions
    m_uniform = np.ones(n) / float(n)
    m_weighted = np.array([102, 102, 102, 102, 1, 1], dtype=float) / 410
    # make the response
    out = StringIO()
    # look at the eigendecomposition of -(1/2)HDH where D is the leaf distance matrix
    HSH = Euclid.edm_to_dccov(Euclid.g_D_b)
    W_HSH, VT_HSH = np.linalg.eigh(HSH)
    print >> out, 'W for -(1/2)HDH of the leaf distance matrix:'
    print >> out, W_HSH
    print >> out, 'VT for -(1/2)HDH of the leaf distance matrix:'
    print >> out, VT_HSH
    # look at the eigendecomposition of S given a degenerate mass distribution on the full tree
    m_degenerate = np.array([.25, .25, .25, .25, 0, 0])
    S = Euclid.edm_to_weighted_cross_product(Euclid.g_D_c, m_degenerate)
    W_S, VT_S = np.linalg.eigh(S)
    print >> out, 'W for -(1/2)(Xi)D(Xi)^T of the full distance matrix with degenerate masses:'
    print >> out, W_S
    print >> out, 'VT for -(1/2)(Xi)D(Xi)^T of the full distance matrix with degenerate masses:'
    print >> out, VT_S
    # look at the effects of various mass distributions on the MDS of the full tree
    for m in (m_uniform, m_weighted):
        # the mass distribution should sum to 1
        if not np.allclose(np.sum(m), 1):
            raise ValueError('masses should sum to 1')
        # to compute the perturbed laplacian matrix first get weighted sums
        v = np.dot(m, A)
        # now divide elementwise by the masses
        v /= m
        # subtract the adjacency matrix from the diagonal formed by elements of this vector
        Lp = np.diag(v) - A
        # now get the eigendecomposition of the pseudoinverse of the perturbed laplacian
        W_Lp_pinv, VT_Lp_pinv = np.linalg.eigh(np.linalg.pinv(Lp))
        # look at the eigendecomposition of the S matrix associated with the distance matrix of this tree
        D = Euclid.g_D_c
        S = Euclid.edm_to_weighted_cross_product(D, m)
        W_S, VT_S = np.linalg.eigh(S)
        print >> out, 'perturbed laplacian:'
        print >> out, Lp
        print >> out, 'm:', m
        print >> out, 'W for the pseudoinverse of the perturbed laplacian:'
        print >> out, W_Lp_pinv
        print >> out, 'VT for the pseudoinverse of the perturbed laplacian:'
        print >> out, VT_Lp_pinv
        print >> out, 'W for the cross product matrix:'
        print >> out, W_S
        print >> out, 'VT for the cross product matrix:'
        print >> out, VT_S
    return out.getvalue().strip()

def produce_CRT():
    (p,q) = produce_p_q()
    n = p*q
    Euler = (p-1)*(q-1) #欧拉函数
    d = Euclid.extended_Euclid(e,Euler)#求出e模Euler的逆元d，e*d=1mod（Euler）
    dP = Euclid.extended_Euclid(e,p-1)
    dQ = Euclid.extended_Euclid(e,q-1)
    qInv = Euclid.extended_Euclid(q,p)
    return (p,q,n,d,dP,dQ,qInv)

def get_response_content(fs):
    # build the newick tree from the string
    tree = NewickIO.parse(fs.tree_string, FelTree.NewickTree)
    nvertices = len(list(tree.preorder()))
    nleaves = len(list(tree.gen_tips()))
    # get ordered ids with the leaves first
    ordered_ids = get_ordered_ids(tree)
    # get the adjacency matrix and the augmented adjacency matrix
    A = np.array(tree.get_affinity_matrix(ordered_ids))
    A_aug = get_augmented_adjacency(A, nleaves, fs.ndups, fs.strength)
    # get the laplacian matrices
    L = Euclid.adjacency_to_laplacian(A)
    L_aug = Euclid.adjacency_to_laplacian(A_aug)
    # get the schur complement
    R = SchurAlgebra.mschur(L, set(range(nleaves, nvertices)))
    R_pinv = np.linalg.pinv(R)
    vals, vecs = EigUtil.eigh(R_pinv)
    # get the scaled Fiedler vector for the Schur complement
    w, v = EigUtil.principal_eigh(R_pinv)
    fiedler = v * math.sqrt(w)
    # get the eigendecomposition of the augmented Laplacian
    L_aug_pinv = np.linalg.pinv(L_aug)
    vals_aug, vecs_aug = EigUtil.eigh(L_aug_pinv)
    # get the scaled Fiedler vector for the augmented Laplacian
    w_aug, v_aug = EigUtil.principal_eigh(L_aug_pinv)
    fiedler_aug = v_aug * math.sqrt(w_aug)
    # report the results
    np.set_printoptions(linewidth=300)
    out = StringIO()
    print >> out, 'Laplacian matrix:'
    print >> out, L
    print >> out
    print >> out, 'Schur complement of Laplacian matrix:'
    print >> out, R
    print >> out
    print >> out, 'scaled Fiedler vector of Schur complement:'
    print >> out, fiedler
    print >> out
    print >> out, 'eigenvalues of pinv of Schur complement:'
    print >> out, vals
    print >> out
    print >> out, 'corresponding eigenvectors of pinv of Schur complement:'
    print >> out, np.array(vecs).T
    print >> out
    print >> out
    print >> out, 'augmented Laplacian matrix:'
    print >> out, L_aug
    print >> out
    print >> out, 'scaled Fiedler vector of augmented Laplacian:'
    print >> out, fiedler_aug
    print >> out
    print >> out, 'eigenvalues of pinv of augmented Laplacian:'
    print >> out, vals_aug
    print >> out
    print >> out, 'rows are eigenvectors of pinv of augmented Laplacian:'
    print >> out, np.array(vecs_aug)
    return out.getvalue()

def get_response_content(fs):
    # build the newick tree from the string
    tree = NewickIO.parse(fs.tree_string, FelTree.NewickTree)
    nvertices = len(list(tree.preorder()))
    nleaves = len(list(tree.gen_tips()))
    # get ordered ids with the leaves first
    ordered_ids = get_ordered_ids(tree)
    # get the distance matrix and the augmented distance matrix
    D = np.array(tree.get_partial_distance_matrix(ordered_ids))
    D_aug = get_augmented_distance(D, nleaves, fs.ndups)
    # get the laplacian matrix
    L = Euclid.edm_to_laplacian(D)
    # get the schur complement
    R = SchurAlgebra.mschur(L, set(range(nleaves, nvertices)))
    R_pinv = np.linalg.pinv(R)
    vals, vecs = EigUtil.eigh(R_pinv)
    # get the scaled Fiedler vector for the Schur complement
    w, v = EigUtil.principal_eigh(R_pinv)
    fiedler = v * math.sqrt(w)
    # get the eigendecomposition of the centered augmented distance matrix
    L_aug_pinv = Euclid.edm_to_dccov(D_aug)
    vals_aug, vecs_aug = EigUtil.eigh(L_aug_pinv)
    # get the scaled Fiedler vector for the augmented Laplacian
    w_aug, v_aug = EigUtil.principal_eigh(L_aug_pinv)
    fiedler_aug = v_aug * math.sqrt(w_aug)
    # report the results
    np.set_printoptions(linewidth=300, threshold=10000)
    out = StringIO()
    print >> out, "Laplacian matrix:"
    print >> out, L
    print >> out
    print >> out, "Schur complement of Laplacian matrix:"
    print >> out, R
    print >> out
    print >> out, "scaled Fiedler vector of Schur complement:"
    print >> out, fiedler
    print >> out
    print >> out, "eigenvalues of pinv of Schur complement:"
    print >> out, vals
    print >> out
    print >> out, "corresponding eigenvectors of pinv of Schur complement:"
    print >> out, np.array(vecs).T
    print >> out
    print >> out
    print >> out, "augmented distance matrix:"
    print >> out, D_aug
    print >> out
    print >> out, "scaled Fiedler vector of augmented Laplacian limit:"
    print >> out, fiedler_aug
    print >> out
    print >> out, "eigenvalues of pinv of augmented Laplacian limit:"
    print >> out, vals_aug
    print >> out
    print >> out, "rows are eigenvectors of pinv of augmented Laplacian limit:"
    print >> out, np.array(vecs_aug)
    return out.getvalue()

def get_response_content(fs):
    # read the lat-lon points from the input
    lines = Util.get_stripped_lines(fs.datalines.splitlines())
    rows = parse_lines(lines)
    latlon_points = []
    city_names = []
    for city, latd, latm, lond, lonm in rows:
        lat = math.radians(GPS.degrees_minutes_to_degrees(latd, latm))
        lon = math.radians(GPS.degrees_minutes_to_degrees(lond, lonm))
        latlon_points.append((lat, lon))
        city_names.append(city)
    npoints = len(latlon_points)
    # start writing the response
    np.set_printoptions(linewidth=200)
    out = StringIO()
    radius = GPS.g_earth_radius_miles
    for dfunc, name in (
            (GPS.get_arc_distance, 'great arc'),
            (GPS.get_euclidean_distance, 'euclidean')):
        # define the edm whose elements are squared euclidean-like distances
        edm = np.zeros((npoints, npoints))
        D = np.zeros((npoints, npoints))
        for i, pointa in enumerate(latlon_points):
            for j, pointb in enumerate(latlon_points):
                D[i, j] = dfunc(pointa, pointb, radius)
                edm[i, j] = D[i, j]**2
        print >> out, name, 'distances:'
        print >> out, D
        print >> out
        print >> out, name, 'EDM:'
        print >> out, edm
        print >> out
        G = Euclid.edm_to_dccov(edm)
        print >> out, name, 'Gower centered matrix:'
        print >> out, G
        print >> out
        spectrum = np.array(list(reversed(sorted(np.linalg.eigvals(G)))))
        print >> out, name, 'spectrum of Gower centered matrix:'
        for x in spectrum:
            print >> out, x
        print >> out
        print >> out, name, 'rounded spectrum:'
        for x in spectrum:
            print >> out, '%.1f' % x
        print >> out
        mds_points = Euclid.edm_to_points(edm)
        print >> out, '2D MDS coordinates:'
        for name, mds_point in zip(city_names, mds_points):
            x = mds_point[0]
            y = mds_point[1]
            print >> out, '\t'.join(str(x) for x in [name, x, y])
        print >> out
        # break between distance methods
        print >> out
    # return the response
    return out.getvalue()

def produce_e_d():
    """产生(e,d)"""
    e = 3
    while True:
        d = Euclid.extended_Euclid(e,m)#e*d=1modm
        if Euclid.gcd(m,e) == 1 and d > 0:
            break
        else:
            e += 2
    return (e,d)

def update_using_laplacian(D, index_set):
    """
    Update the distance matrix by summing rows and columns of the removed indices.
    @param D: the distance matrix
    @param index_set: the set of indices that will be removed from the updated distance matrix
    @return: an updated distance matrix
    """
    L = Euclid.edm_to_laplacian(D)
    L_small = SchurAlgebra.mmerge(L, index_set)
    D_small = Euclid.laplacian_to_edm(L_small)
    return D_small

def produce_e_d():
    """产生(e, d)"""
    # Generate a number e so that gcd(e, m) = 1, start with e = 3
    e = 3
    while True:
        d = Euclid.extended_Euclid(e, m)
        if Euclid.gcd(m, e) == 1 and d > 0:
            break
        else:
            e += 2
    return (e, d)

def get_splits(initial_distance_matrix, split_function, update_function, on_label_split=None):
    """
    This is the most external of the functions in this module.
    Get the set of splits implied by the tree that would be reconstructed.
    @param initial_distance_matrix: a distance matrix
    @param split_function: takes a distance matrix and returns an index split
    @param update_function: takes a distance matrix and an index subset and returns a distance matrix
    @param on_label_split: notifies the caller of the label split induced by an index split
    @return: a set of splits
    """
    n = len(initial_distance_matrix)
    # keep a stack of (label_set_per_vertex, distance_matrix) pairs
    initial_state = ([set([i]) for i in range(n)], initial_distance_matrix)
    stack = [initial_state]
    # process the stack in a depth first manner, building the split set
    label_split_set = set()
    while stack:
        label_sets, D = stack.pop()
        # if the matrix is small then we are done
        if len(D) < 4:
            continue
        # split the indices using the specified function
        try:
            index_split = split_function(D)
            # convert the index split to a label split
            label_split = index_split_to_label_split(index_split, label_sets)
            # notify the caller if a callback is requested
            if on_label_split:
                on_label_split(label_split)
            # add the split to the master set of label splits
            label_split_set.add(label_split)
            # for large matrices create the new label sets and the new conformant distance matrices
            a, b = index_split
            for index_selection, index_complement in ((a, b), (b, a)):
                if len(index_complement) > 2:
                    next_label_sets = SchurAlgebra.vmerge(label_sets, index_selection)
                    next_D = update_function(D, index_selection)
                    next_state = (next_label_sets, next_D)
                    stack.append(next_state)
        except DegenerateSplitException, e:
            # we cannot recover from a degenerate split unless there are more than four indices
            if len(D) <= 4:
                continue
            # with more than four indices we can fall back to partial splits
            index_set = set([e.index])
            # get the next label sets
            next_label_sets = SchurAlgebra.vdelete(label_sets, index_set)
            # get the next conformant distance matrix by schur complementing out the offending index
            L = Euclid.edm_to_laplacian(D)
            L_small = SchurAlgebra.mschur(L, index_set)
            next_D = Euclid.laplacian_to_edm(L_small)
            next_state = (next_label_sets, next_D)
            stack.append(next_state)