Aarondard5/constrained_multiOT

1from libpysal import weights, examples2from libpysal.cg import voronoi_frames3import geopandas4import networkx as nx5import numpy as np6from itertools import combinations7from math import sqrt8import random9 10 11def planar_graph(12    nnode,13    G0=None,14    subN=None,15    L_max=None,16    L_min=None,17    domain=(0, 0, 1, 1),18    metric=None,19    seed=10,20    connected_component=True,21):22    """Returns a planar random graph modified so that 23    only egdes within a distance are possible"""24 25    prng = np.random.RandomState(seed)26    if G0 is None:27        G = nx.Graph()28        G.add_nodes_from(np.arange(nnode))29        (xmin, ymin, xmax, ymax) = domain30        # Each node gets a uniformly random position in the given rectangle.31        pos = {v: (prng.uniform(xmin, xmax), prng.uniform(ymin, ymax)) for v in G}32        nx.set_node_attributes(G, pos, "pos")33    else:34        print("using G")35        G = nx.Graph()36        G.add_nodes_from(G0.nodes(data=True))37        nnode = G.number_of_nodes()38        pos = nx.get_node_attributes(G, "pos")39 40    if subN is None:41        subN = nnode42 43    if subN < nnode:44        random.seed(seed)45        node2remove = random.sample(list(G.nodes()), k=nnode - subN)46        G.remove_nodes_from(node2remove)47 48    ### Extracts planar graph49    nodes = list(G.nodes())50    coordinates = np.array([(pos[n][0], pos[n][1]) for n in list(G.nodes())])51    cells, generators = voronoi_frames(coordinates, clip="convex hull")52    delaunay = weights.Rook.from_dataframe(cells)53    G = delaunay.to_networkx()54    G = nx.relabel_nodes(G, {i: nodes[i] for i in range(G.number_of_nodes())})55    positions = {n: coordinates[i] for i, n in enumerate(list(G.nodes()))}56    nx.set_node_attributes(G, positions, "pos")57 58    # If no distance metric is provided, use Euclidean distance.59    if metric is None:60        metric = euclidean61 62    if L_max is None:63        L_max = max(metric(x, y) for x, y in combinations(positions.values(), 2))64    if L_min is None:65        L_min = 066 67    def dist(u, v):68        return metric(positions[u], positions[v])69 70    edges2remove = [71        e72        for e in list(G.edges())73        if np.logical_or(dist(*e) > L_max, dist(*e) < L_min) == True74    ]75    G.remove_edges_from(edges2remove)76 77    if connected_component == True:78        Gc = max(nx.connected_components(G), key=len)79        nodes_to_remove = set(G.nodes()).difference(Gc)80        G.remove_nodes_from(list(nodes_to_remove))81 82    return G83 84 85def euclidean(x, y):86    """Returns the Euclidean distance between the vectors 87    ``x`` and ``y``.88 89    Each of ``x`` and ``y`` can be any iterable of numbers. The90    iterables must be of the same length."""91 92    return np.sqrt(sum((a - b) ** 2 for a, b in zip(x, y)))93 94 95