Aarondard5/constrained_multiOT
0
1import networkx as nx2import random3import numpy as np4import matplotlib.pyplot as plt5import dynamics_multiot as dyn6 7from itertools import islice8 9print('using tools.py def newo...')10 11def find_central_nodes(G, pos=None, pos_center_of_mass=None):12 13 if pos is None:14 pos = nx.get_node_attributes(G, "pos")15 if pos_center_of_mass is None:16 positions = np.array(list(pos.values()))17 pos_center_of_mass = np.mean(positions, axis=0)18 print(f"pos_center_of_mass:{pos_center_of_mass}")19 nodes = list(G.nodes)20 distances = [euclidean(pos[n], pos_center_of_mass) for n in nodes]21 nid = np.argmin(distances)22 n = nodes[nid]23 return n24 25def extract_central_node_l2(graph1, node2type=None):26 """27 Only considers nodes in layer 228 """29 if node2type is None:30 node2type = {}31 for n, d in graph1.nodes(data=True):32 if d["ntype"] == "intra":33 nlayer = np.unique(34 [d["etype"] for u, v, d in graph1.edges(n, data=True)]35 )36 assert len(nlayer) == 137 node2type[n] = 038 elif d["ntype"] == "inter":39 node2type[n] = 140 else:41 node2type[n] = 242 nodes = list(graph1.nodes)43 pos = nx.get_node_attributes(graph1, "pos")44 45 positions = np.array([pos[n] for n in pos if node2type[n] == 2])46 pos_center_of_mass = np.mean(positions, axis=0)47 distances = [euclidean(pos[n], pos_center_of_mass) for n in nodes]48 nid = np.argsort(distances)49 50 nid_l2 = [i for i in nid if node2type[nodes[i]] == 2]51 pos_center_of_mass = pos[nodes[nid_l2[0]]]52 53 return pos_center_of_mass, nodes[nid_l2[0]]54 55def forcing_generate_real(G, p=0.0, weigth=10, pos=None, pos_center_of_mass=(0.5, 0.5), seed=10):56 57 prng = np.random.RandomState(seed)58 59 n_center = find_central_nodes(G, pos=pos, pos_center_of_mass=pos_center_of_mass)60 61 nodes = sorted(G.nodes())62 nodes = [n for n in nodes if n != n_center]63 forcing = []64 for source in nodes:65 r = prng.rand()66 if r < p:67 # deterministic candidate list68 candidates = [n for n in nodes if n != source]69 sink = prng.choice(candidates)70 else:71 sink = n_center72 73 forcing.append([source, sink, weigth])74 return forcing75 76def forcing_generate(G, p=0.0, weigth=10, pos=None, pos_center_of_mass=(0.5, 0.5), seed=1077):78 """G is the graph in the first layer"""79 80 prng = np.random.RandomState(seed)81 random.seed(seed) 82 n_center = find_central_nodes(G, pos=pos, pos_center_of_mass=pos_center_of_mass)83 84 nodes = set(G.nodes()) - {n_center} # all nodes except central one85 forcing = []86 for source in nodes:87 if source != n_center:88 r = prng.rand()89 if r < p: # rewire: extract a random node90 sink = random.choice(list(nodes - {source}))91 forcing.append([source, sink, weigth])92 else:93 forcing.append([source, n_center, weigth])94 95 return forcing96 97def euclidean(x, y):98 """Returns the Euclidean distance between the vectors ``x`` and ``y``.99 100 Each of ``x`` and ``y`` can be any iterable of numbers. The101 iterables must be of the same length.102 103 """104 return np.sqrt(sum((a - b) ** 2 for a, b in zip(x, y)))105 106 107def plot_results(108 G,109 graph,110 opttdens,111 optpot,112 optflux,113 length,114 flag_trim=False,115 ns=20,116 outfigure=None,117 tau=None,118 dpi=300,119 draw_source_sink=None,120 wl=False,121 w0=1.0,122 edge_width="linear",123 figsize=(15, 5),124 colors_map={125 0: "b",126 1: "r",127 2: "g",128 3: "magenta",129 "inter": "black",130 "inter-super": "black",131 },132):133 134 if flag_trim == True:135 tau = tau136 (137 graph_final,138 opttdens_final,139 flux_norm_final,140 length_final,141 flux_final,142 ) = dyn.abs_trimming_dyn(graph, opttdens, optpot, length, tau=tau)143 else:144 graph_final = nx.Graph(graph)145 opttdens_final = opttdens.copy()146 flux_norm_final = optflux.copy()147 148 # print(flux_norm_final)149 """150 Assign colors151 """152 etypes = nx.get_edge_attributes(graph_final, "etype")153 nlayer = len(set(etypes.values()).difference(set(["inter", "inter-super"])))154 for e in list(graph_final.edges):155 graph_final[e[0]][e[1]]["color"] = colors_map[etypes[e]]156 colors = nx.get_edge_attributes(graph_final, "color")157 158 print(type(flux_norm_final))159 """160 Assign edge widths161 """162 if edge_width == "linear":163 widths = [164 0 + w0 * flux_norm_final[idx]165 for idx, e in enumerate(list(graph_final.edges()))166 ]167 else:168 widths = [169 0 + w0 * np.log(1 + flux_norm_final[idx])170 for idx, e in enumerate(list(graph_final.edges()))171 ]172 173 """for idx, e in enumerate(list(graph_final.edges(data=True))):174 if e[2]['etype'] == 'inter':175 widths[idx] = 0.1176 if e[2]['etype'] == 'inter-super':177 widths[idx] = 0.1"""178 179 """180 Build graph with the solution181 """182 G_plot = nx.Graph(graph_final)183 # G_plot = nx.relabel_nodes(G_plot, nodeId2Name)184 # colors = nx.get_edge_attributes(G_plot,'color')185 G_plot.number_of_edges(), G_plot.number_of_nodes()186 weights = {187 e: flux_norm_final[idx] for idx, e in enumerate(list(graph_final.edges()))188 }189 inv_weights = {190 e: 1.0 / (flux_norm_final[idx] + 1e-12)191 for idx, e in enumerate(list(graph_final.edges()))192 }193 widths_dict = {e: widths[idx] for idx, e in enumerate(list(graph_final.edges()))}194 nx.set_edge_attributes(G_plot, weights, "flux")195 nx.set_edge_attributes(G_plot, inv_weights, "inv_flux")196 nx.set_edge_attributes(G_plot, widths_dict, "widths")197 198 """199 Plot200 """201 202 figsize = figsize203 204 fig, axes = plt.subplots(nrows=1, ncols=2, figsize=figsize)205 ax = axes.flatten()206 207 edges = G.edges()208 colors = "b" # nx.get_edge_attributes(G,'color')209 210 pos = nx.get_node_attributes(G, "pos")211 nx.draw(212 G,213 pos,214 with_labels=wl,215 node_size=ns,216 edge_color=colors,217 font_weight="bold",218 ax=ax[0],219 )220 221 """edges = H.edges()222 colors = 'r' #[H[u][v]['color'] for u,v in edges]223 pos = nx.get_node_attributes(H,'pos')"""224 225 # nx.draw(H, pos, with_labels=wl, node_size = ns, edge_color=colors, font_weight="bold", node_color='grey', ax=ax[1])226 # ax[1].set_axis_off()227 ax[1].set_xlim([1.0 * x for x in ax[0].get_xlim()])228 ax[1].set_ylim([1.0 * y for y in ax[0].get_ylim()])229 230 edges = G_plot.edges()231 pos_plot = nx.get_node_attributes(G_plot, "pos")232 colors = list(nx.get_edge_attributes(G_plot, "color").values())233 # for n in pos_plot.keys():234 # if G_plot.nodes[n]['ntype'] == 'inter': pos_plot[n] = G_plot.nodes[n]['pos'] * (1 + 0.2 * np.random.rand(2))235 # if G_plot.nodes[n]['ntype'] == 'super': pos_plot[n] = G_plot.nodes[n]['pos'] * (1 + 0.2 * np.random.rand(2))236 237 # nx.draw(G_plot,pos_plot, node_size = ns, with_labels = wl, edge_color='grey', font_weight="bold", alpha=0.15, ax=ax[2])238 nx.draw(239 G_plot,240 pos_plot,241 node_size=ns,242 with_labels=wl,243 edge_color=colors,244 font_weight="bold",245 width=widths,246 ax=ax[1],247 )248 if draw_source_sink is not None:249 nx.draw_networkx_nodes(250 G_plot,251 pos_plot,252 nodelist=draw_source_sink,253 node_size=ns * 100,254 label=True,255 node_color=["green", "magenta"],256 ax=ax[1],257 )258 # print('edgelist',widths )259 260 if outfigure is not None:261 plt.savefig(262 outfigure + ".png",263 dpi=dpi,264 format="png",265 bbox_inches="tight",266 pad_inches=0.1,267 )268 269 plt.show()270 271 return G_plot272 273 274def plot_results_all_algos(275 graph,276 flux_norm1,277 w1=0.01,278 figsize=(20, 7),279 ns=10,280 outfigure=None,281 colors_map={282 0: "b",283 1: "r",284 2: "g",285 3: "magenta",286 "inter": "black",287 "inter-super": "black",288 },289 algo_labels={"MultiOT_UConstr": "UnConstrained", "MultiOT_capacity": "Constraind"},290 fs_title=20,291 dpi=300,292):293 """294 Build G_plot295 """296 G_plot = nx.Graph(graph)297 etypes = nx.get_edge_attributes(graph, "etype")298 nlayer = len(set(etypes.values()).difference(set(["inter", "inter-super"])))299 for e in list(graph.edges):300 G_plot[e[0]][e[1]]["color"] = colors_map[etypes[e]]301 302 algos = list(flux_norm1.keys())303 fig, axes = plt.subplots(nrows=1, ncols=len(algos), figsize=figsize)304 ax = axes.flatten()305 pos_plot = nx.get_node_attributes(G_plot, "pos")306 colors = list(nx.get_edge_attributes(G_plot, "color").values())307 308 for a_id, a in enumerate(algos):309 widths = [w1 * flux_norm1[a][idx] for idx, e in enumerate(list(graph.edges()))]310 nx.draw_networkx(311 G_plot,312 pos_plot,313 node_size=ns,314 with_labels=False,315 edge_color=colors,316 font_weight="bold",317 width=widths,318 ax=ax[a_id],319 )320 ax[a_id].set_title(algo_labels[a], fontsize=fs_title)321 322 if outfigure is not None:323 plt.savefig(324 outfigure + ".png",325 dpi=dpi,326 format="png",327 bbox_inches="tight",328 pad_inches=0.1,329 )330 331 plt.show()332 333 return G_plot, plt.gcf()334 335 336def plot_individual_G_final(337 G_plot,338 graph_final,339 ns=0,340 outfigure=None,341 colors=None,342 wl=False,343 figsize=(5, 5),344 widths=None,345 dpi=300,346 alpha=None,347 beta0=None,348 beta1=None,349 colors_map={350 0: "b",351 1: "r",352 2: "g",353 3: "magenta",354 "inter": "black",355 "inter-super": "black",356 },357):358 359 fig = plt.figure(figsize=figsize)360 pos_plot = nx.get_node_attributes(G_plot, "pos")361 362 """add text"""363 x_min, y_min = np.array(list(pos_plot.values())).min(axis=0)364 x_max, y_max = np.array(list(pos_plot.values())).max(axis=0)365 delta_x = x_max - x_min366 delta_y = y_max - y_min367 368 if colors is None:369 colors = list(nx.get_edge_attributes(G_plot, "color").values())370 if widths is None:371 widths = list(nx.get_edge_attributes(G_plot, "widths").values())372 pos = nx.get_node_attributes(graph_final, "pos")373 nx.draw(374 graph_final,375 pos,376 node_size=0,377 node_color="grey",378 edge_color="grey",379 alpha=0.1,380 with_labels=False,381 )382 # print(graph_final.nodes())383 nx.draw(384 G_plot,385 pos_plot,386 node_size=ns,387 with_labels=False,388 edge_color=colors,389 font_weight="bold",390 width=widths,391 )392 393 props = dict(boxstyle="round", facecolor="white", alpha=1)394 # text = r'$\beta_1=$'+str(beta0)+r', $\beta_2=$'+str(beta1) +r', $w_2=$'+str(alpha)395 # plt.gcf().text(0.25 * delta_x , 1.05 * delta_y , text, fontsize=15,bbox=props)396 plt.xlim([x_min, x_max])397 plt.ylim([y_min, y_max])398 399 if outfigure is not None:400 plt.savefig(401 outfigure + ".png",402 dpi=dpi,403 format="png",404 bbox_inches="tight",405 pad_inches=0.1,406 )407 408 plt.show()409 return plt.gcf()410 411 412def gini_coefficient(x):413 """Compute Gini coefficient of array of values"""414 if isinstance(x, (dict)):415 x = np.array(list(x.values()))416 # diffsum = 0417 # for i, xi in enumerate(x):418 # diffsum += np.sum(np.abs(np.subtract.outer(xi, x)))419 diffsum = np.sum(np.abs(np.subtract.outer(x, x)))420 gini = 0.5 * (diffsum / (len(x) ** 2 * np.mean(x)))421 422 return gini423 424 425def flt(x, d=1):426 return round(x, d)427 428'''429def k_simple_paths(G, source, target, k, weight=None):430 return list(islice(nx.all_simple_paths(G, source, target), k))431 432 433def plot_results2(434 G,435 H,436 graph,437 opttdens,438 optpot,439 optflux,440 length,441 M,442 flag_trim=False,443 ns=20,444 outfigure=None,445 tau=None,446 dpi=300,447 draw_source_sink=None,448 wl=False,449 w0=1.0,450 edge_width="linear",451 figsize=(15, 5),452 colors_map={453 0: "b",454 1: "r",455 2: "g",456 3: "magenta",457 "inter": "black",458 "inter-super": "black",459 },460):461 462 if flag_trim == True:463 tau = tau464 (465 graph_final,466 opttdens_final,467 flux_norm_final,468 length_final,469 flux_final,470 ) = dyn.abs_trimming_dyn(graph, opttdens, optpot, length, tau=tau)471 else:472 graph_final = nx.Graph(graph)473 opttdens_final = opttdens.copy()474 flux_norm_final = optflux.copy()475 476 # print(flux_norm_final)477 """478 Assign colors479 """480 etypes = nx.get_edge_attributes(graph_final, "etype")481 nlayer = len(set(etypes.values()).difference(set(["inter", "inter-super"])))482 for e in list(graph_final.edges):483 graph_final[e[0]][e[1]]["color"] = colors_map[etypes[e]]484 colors = nx.get_edge_attributes(graph_final, "color")485 486 print(type(flux_norm_final))487 """488 Assign edge widths489 """490 if edge_width == "linear":491 widths = [492 0 + w0 * flux_norm_final[idx]493 for idx, e in enumerate(list(graph_final.edges()))494 ]495 else:496 widths = [497 0 + w0 * np.log(1 + flux_norm_final[idx])498 for idx, e in enumerate(list(graph_final.edges()))499 ]500 501 for idx, e in enumerate(list(graph_final.edges(data=True))):502 if e[2]["etype"] == "inter":503 widths[idx] = 0.1504 if e[2]["etype"] == "inter-super":505 widths[idx] = 0.1506 507 """508 Build graph with the solution509 """510 G_plot = nx.Graph(graph_final)511 # G_plot = nx.relabel_nodes(G_plot, nodeId2Name)512 # colors = nx.get_edge_attributes(G_plot,'color')513 G_plot.number_of_edges(), G_plot.number_of_nodes()514 weights = {515 e: flux_norm_final[idx] for idx, e in enumerate(list(graph_final.edges()))516 }517 inv_weights = {518 e: 1.0 / (flux_norm_final[idx] + 1e-12)519 for idx, e in enumerate(list(graph_final.edges()))520 }521 widths_dict = {e: widths[idx] for idx, e in enumerate(list(graph_final.edges()))}522 nx.set_edge_attributes(G_plot, weights, "flux")523 nx.set_edge_attributes(G_plot, inv_weights, "inv_flux")524 nx.set_edge_attributes(G_plot, widths_dict, "widths")525 526 """527 Plot528 """529 530 figsize = figsize531 532 fig, axes = plt.subplots(nrows=1, ncols=3, figsize=figsize)533 ax = axes.flatten()534 535 edges = G.edges()536 colors = "b" # nx.get_edge_attributes(G,'color')537 538 pos = nx.get_node_attributes(G, "pos")539 nx.draw(540 G,541 pos,542 with_labels=wl,543 node_size=ns,544 edge_color=colors,545 font_weight="bold",546 ax=ax[0],547 )548 549 edges = H.edges()550 colors = "r" # [H[u][v]['color'] for u,v in edges]551 pos = nx.get_node_attributes(H, "pos")552 553 nx.draw(554 H,555 pos,556 with_labels=wl,557 node_size=ns,558 edge_color=colors,559 font_weight="bold",560 node_color="grey",561 ax=ax[1],562 )563 # ax[1].set_axis_off()564 ax[1].set_xlim([1.0 * x for x in ax[0].get_xlim()])565 ax[1].set_ylim([1.0 * y for y in ax[0].get_ylim()])566 567 edges = G_plot.edges()568 pos_plot = nx.get_node_attributes(G_plot, "pos")569 colors = list(nx.get_edge_attributes(G_plot, "color").values())570 # for n in pos_plot.keys():571 # if G_plot.nodes[n]['ntype'] == 'inter': pos_plot[n] = G_plot.nodes[n]['pos'] * (1 + 0.2 * np.random.rand(2))572 # if G_plot.nodes[n]['ntype'] == 'super': pos_plot[n] = G_plot.nodes[n]['pos'] * (1 + 0.2 * np.random.rand(2))573 574 # nx.draw(G_plot,pos_plot, node_size = ns, with_labels = wl, edge_color='grey', font_weight="bold", alpha=0.15, ax=ax[2])575 nx.draw(576 G_plot,577 pos_plot,578 node_size=ns,579 with_labels=wl,580 edge_color=colors,581 font_weight="bold",582 width=widths,583 ax=ax[2],584 )585 if draw_source_sink is not None:586 nx.draw_networkx_nodes(587 G_plot,588 pos_plot,589 nodelist=draw_source_sink,590 node_size=ns * 100,591 label=True,592 node_color=["green", "magenta"],593 ax=ax[2],594 )595 # print('edgelist',widths )596 ax[2].set_xlim([1.0 * x for x in ax[0].get_xlim()])597 ax[2].set_ylim([1.0 * y for y in ax[0].get_ylim()])598 599 if outfigure is not None:600 plt.savefig(601 outfigure + ".png",602 dpi=dpi,603 format="png",604 bbox_inches="tight",605 pad_inches=0.1,606 )607 608 plt.show()609 610 return G_plot611 612def plot_gmu(613 g_func,614 not_in_C_id,615 model,616 capacity,617 indices=None,618 k_i=5,619 figsize=(20, 7),620 int_ticks=False,621 xlab="Edges",622 outfig=None,623):624 625 """626 plot function g(mu) = c_e - mu_e627 """628 629 plt.figure(figsize=figsize)630 fs = 23631 632 indices = None633 plt.subplot(1, 2, 1)634 if indices is None:635 plt.plot(g_func[k_i:])636 else:637 plt.plot(indices[k_i:], g_func[k_i:])638 plt.xlabel(xlab)639 plt.ylabel(r"$g = c_e - \mu_e$", fontsize=fs)640 if int_ticks:641 plt.xaxis.set_major_locator(MaxNLocator(integer=True))642 plt.grid()643 plt.title(f"{model}:" + r": [$g(\mu_e) = c_e - \mu_e$]", fontsize=fs)644 plt.axline((0, capacity), slope=0.0, color="b", linestyle=(0, (5, 5)))645 plt.text(-28, capacity, r"$c_e$", {"color": "b", "fontsize": 25})646 647 plt.subplot(1, 2, 2)648 if indices is None:649 plt.plot(g_func[not_in_C_id][k_i:])650 else:651 plt.plot(indices[k_i:], g_func[not_in_C_id][k_i:])652 plt.xlabel(xlab)653 plt.ylabel(r"$g(\mu_e) = c_e - \mu_e < 0$", fontsize=fs)654 if int_ticks:655 plt.xaxis.set_major_locator(MaxNLocator(integer=True))656 plt.grid()657 plt.title(f"{model}:" + r"[$g(\mu_e) = c_e - \mu_e < 0$]", fontsize=fs)658 659 plt.axline((0, capacity), slope=0.0, color="b", linestyle=(0, (5, 5)))660 661 # plt.text(-6, capacity, r'$c_e$', {'color': 'b', 'fontsize': 25})662 663 plt.tight_layout()664 665 if outfig is not None:666 plt.savefig(outfig, dpi=300)667 668 return plt.show()669 670'''671def plot_real_all_algos(672 graph,673 flux_norm1,674 w1=0.01,675 w2=400,676 figsize=(20, 7),677 ns=10,678 outfigure=None,679 colors_map={680 0: "b",681 1: "r",682 2: "g",683 3: "magenta",684 "inter": "black",685 "inter-super": "black",686 },687 algo_labels={"MultiOT_UConstr": "UnConstrained", "MultiOT_capacity": "Constraind"},688 fs_title=20,689 dpi=300,690):691 692 """693 Build G_plot694 """695 G_plot = nx.Graph(graph)696 etypes = nx.get_edge_attributes(graph, "etype")697 nlayer = len(set(etypes.values()).difference(set(["inter", "inter-super"])))698 for e in list(graph.edges):699 G_plot[e[0]][e[1]]["color"] = colors_map[etypes[e]]700 701 plt.figure(figsize=(7, 7))702 pos_plot = nx.get_node_attributes(G_plot, "pos")703 x_min, y_min = np.array(list(pos_plot.values())).min(axis=0)704 x_max, y_max = np.array(list(pos_plot.values())).max(axis=0)705 delta_x = x_max - x_min706 delta_y = y_max - y_min707 708 epsilon = 0.3 * w1709 algos = list(flux_norm1.keys())710 fig, axes = plt.subplots(nrows=1, ncols=len(algos), figsize=figsize)711 ax = axes.flatten()712 pos_plot = nx.get_node_attributes(G_plot, "pos")713 colors = list(nx.get_edge_attributes(G_plot, "color").values())714 ncolors = list(nx.get_node_attributes(G_plot, "ncolor").values())715 716 for a_id, a in enumerate(algos):717 widths_nodes = [718 w2 * np.mean(traffic_u[a][idx]) for idx, e in enumerate(list(graph.nodes()))719 ]720 widths = [721 w1 * np.mean(traffic_e[a][idx]) for idx, e in enumerate(list(graph.edges()))722 ]723 nx.draw_networkx(724 G_plot,725 pos_plot,726 node_size=widths_nodes,727 with_labels=False,728 node_color=ncolors,729 edge_color=colors,730 font_weight="bold",731 width=widths,732 ax=ax[a_id],733 )734 nx.draw(735 G_plot,736 pos_plot,737 node_size=0,738 with_labels=False,739 edge_color=colors,740 font_weight="bold",741 width=widths,742 ax=ax[a_id],743 )744 props = dict(boxstyle="round", facecolor="white", alpha=1)745 roads_bordeaux.plot(ax=ax[a_id], color="#d9d9d9", alpha=0.5, zorder=0)746 747 ax[a_id].set_xlim([x_min, x_max])748 ax[a_id].set_ylim([y_min, y_max])749 ax[a_id].set_title(algo_labels[a], fontsize=fs)750 751 if outfigure is not None:752 plt.savefig(753 outfigure + ".png",754 dpi=dpi,755 format="png",756 bbox_inches="tight",757 pad_inches=0.1,758 )759 760 plt.show()761 762 return G_plot, plt.gcf()763 764