From a89c6d4833eafaf28fb96b044b73f964c824009f Mon Sep 17 00:00:00 2001 From: Yves Biener Date: Thu, 9 Apr 2026 08:48:07 +0200 Subject: [PATCH] mod update corresponding examples --- comparison.py | 88 ++++++++++++++++----------------------- comparison_edge_scores.py | 40 +++++++----------- diff_comparison.py | 35 ++++++++-------- distance_types.py | 20 +++++---- src/fitting.py | 2 + src/plot.py | 81 +++++++++++++++++++++++++++++++---- 6 files changed, 155 insertions(+), 111 deletions(-) diff --git a/comparison.py b/comparison.py index b0e8eac..2d58203 100644 --- a/comparison.py +++ b/comparison.py @@ -2,6 +2,7 @@ import math import matplotlib.pyplot as plt import numpy as np +import squidpy as sq from graph_tool.all import * from src import centrality @@ -9,6 +10,23 @@ from src import plot from src import fitting +def merfish(): + """ + Merfish dataset from `squidpy`. + """ + adata = sq.datasets.merfish() + adata = adata[adata.obs.Bregma == -9].copy() + return adata + + +def mibitof(): + """ + Mibitof dataset from `squidpy`. + """ + adata = sq.datasets.mibitof() + return adata + + def degree(g, weight): # VertexPropertyMap vp = g.new_vertex_property("double") @@ -25,6 +43,8 @@ def leverage(g, weight): li = 0.0 neighbours = g.get_all_neighbours(v) ki = len(neighbours) + if ki == 0: + continue # sum for nv in neighbours: other_neighbours = g.get_all_neighbours(nv) @@ -48,46 +68,6 @@ def random_graph(n=5000, seed=None): return rng.random((n, 2)), seed -def random_graph_favor_border(n=3000, seed = None): - if seed is None: - import secrets - seed = secrets.randbits(128) - rng = np.random.default_rng(seed=seed) - vps = np.zeros((n, 2)) - for i in range(0, n): - r_x = rng.random() - if rng.random() > 0.5: - while (r_x > 0.3 and r_x < 0.7): - r_x = rng.random() - r_y = rng.random() - if rng.random() > 0.5: - while (r_y > 0.3 and r_y < 0.7): - r_y = rng.random() - vps[i][0] = r_x - vps[i][1] = r_y - return vps, seed - - -def random_graph_favor_center(n=3000, seed = None): - if seed is None: - import secrets - seed = secrets.randbits(128) - rng = np.random.default_rng(seed=seed) - vps = np.zeros((n, 2)) - for i in range(0, n): - r_x = rng.random() - if rng.random() > 0.7: - while (r_x < 0.4 or r_x > 0.6): - r_x = rng.random() - r_y = rng.random() - if rng.random() > 0.7: - while (r_y < 0.4 or r_y > 0.6): - r_y = rng.random() - vps[i][0] = r_x - vps[i][1] = r_y - return vps, seed - - def spatial_graph(adata): """ Generate the spatial graph using delaunay for the given `adata`. @@ -102,6 +82,7 @@ def spatial_graph(adata): weight[e] = math.sqrt(sum(map(abs, pos[e.source()].a - pos[e.target()].a)))**2 return g, weight + def apply(g, seed, weight, convex_hull, ax, method, method_name): # calculate centrality values vp = None @@ -147,8 +128,9 @@ def apply(g, seed, weight, convex_hull, ax, method, method_name): # - apply centrality measure to the next axis # - Draw the corresponding resulting models into a grid # -points, seed = random_graph(n=5000) -g, weight = spatial_graph(points) +# points, seed = random_graph(n=5000) +adata = mibitof() +g, weight = spatial_graph(adata.obsm['spatial']) g = GraphView(g) # calculate convex hull convex_hull = centrality.convex_hull(g) @@ -157,23 +139,23 @@ convex_hull = centrality.convex_hull(g) fig_graph, ax_graph = plt.subplots(figsize=(15, 12)) # draw without any centrality measure `vp` vp = g.new_vertex_property("double") -plot.graph_plot(fig_graph, ax_graph, g, vp, convex_hull, f"Pointcloud (seed: {seed})") -fig_graph.savefig("Pointcloud_graph.svg", format='svg') +plot.graph_plot(fig_graph, ax_graph, g, vp, convex_hull, f"mibitof") +fig_graph.savefig(f"mibitof_graph.svg", format='svg') fig = plt.figure(figsize=(15, 12)) row1, row2 = fig.subplots(2, 4) ax1, ax2, ax3, ax4 = row1 # TODO select corresponding centrality measure method -apply(g, seed, weight, convex_hull, ax1, closeness, "Closeness") -apply(g, seed, weight, convex_hull, ax2, pagerank, "PageRank") -apply(g, seed, weight, convex_hull, ax3, betweenness, "Betweeness") -apply(g, seed, weight, convex_hull, ax4, eigenvector, "Eigenvector") +apply(g, None, weight, convex_hull, ax1, closeness, "Closeness") +apply(g, None, weight, convex_hull, ax2, pagerank, "PageRank") +apply(g, None, weight, convex_hull, ax3, betweenness, "Betweeness") +apply(g, None, weight, convex_hull, ax4, eigenvector, "Eigenvector") ax1, ax2, ax3, ax4 = row2 -apply(g, seed, weight, convex_hull, ax1, katz, "Katz") -apply(g, seed, weight, convex_hull, ax2, hits, "Hits") -apply(g, seed, weight, convex_hull, ax3, leverage, "Leverage") -apply(g, seed, weight, convex_hull, ax4, degree, "Degree") +apply(g, None, weight, convex_hull, ax1, katz, "Katz") +apply(g, None, weight, convex_hull, ax2, hits, "Hits") +apply(g, None, weight, convex_hull, ax3, leverage, "Leverage") +apply(g, None, weight, convex_hull, ax4, degree, "Degree") -fig.savefig(f"Comparison_Pointcloud.svg", format='svg') +fig.savefig(f"Comparison_node_centralities_mibitof_.svg", format='svg') diff --git a/comparison_edge_scores.py b/comparison_edge_scores.py index a0cbb2d..2d20cf5 100644 --- a/comparison_edge_scores.py +++ b/comparison_edge_scores.py @@ -85,23 +85,18 @@ def spatial_graph(adata): def apply(g, seed, weight, convex_hull, ax, method, method_name): # calculate centrality values - vp = None + ep = None if method_name == "Betweeness": vp, ep = method(g, weight=weight) - elif method_name == "Eigenvector": - ep, vp = method(g, weight=weight) - elif method_name == "Hits": - ep, vp, hub_centrality = method(g, weight=weight) else: - vp = method(g, weight=weight) - vp.a = np.nan_to_num(vp.a) # correct floating point values + ep = method(g, weight=weight) + ep.a = np.nan_to_num(ep.a) # correct floating point values # normalization - min_val, max_val = vp.a.min(), vp.a.max() - vp.a = (vp.a - min_val) / (max_val - min_val) + min_val, max_val = ep.a.min(), ep.a.max() + ep.a = (ep.a - min_val) / (max_val - min_val) - # generate model based on convex hull and associated centrality values - quantification = plot.quantification_data(g, vp, convex_hull) + quantification = plot.quantification_data_edges(g, ep, convex_hull) # optimize model's piece-wise linear function d = quantification[:, 0] @@ -129,16 +124,16 @@ def apply(g, seed, weight, convex_hull, ax, method, method_name): # - Draw the corresponding resulting models into a grid # points, seed = random_graph(n=5000) -g, weight = spatial_graph(adata.obsm['spatial']) +g, weight = spatial_graph(points) g = GraphView(g) # calculate convex hull convex_hull = centrality.convex_hull(g) # plot graph with convex_hull fig_graph, ax_graph = plt.subplots(figsize=(15, 12)) -# draw without any centrality measure `vp` -vp = g.new_vertex_property("double") -plot.graph_plot(fig_graph, ax_graph, g, vp, convex_hull, f"Pointcould (seed: {seed})") +# draw without any centrality measure `ep` +ep = g.new_edge_property("double") +plot.graph_plot(fig_graph, ax_graph, g, ep, convex_hull, f"Pointcould (seed: {seed})", True) # draw edges fig_graph.savefig(f"comparison_edge_scores_artificial_graph.svg", format='svg') fig = plt.figure(figsize=(15, 12)) @@ -148,15 +143,12 @@ row1, row2 = fig.subplots(2, 4) # - some share similarities to the node based counter parts ax1, ax2, ax3, ax4 = row1 -apply(g, None, weight, convex_hull, ax1, closeness, "Closeness") -apply(g, None, weight, convex_hull, ax2, pagerank, "PageRank") -apply(g, None, weight, convex_hull, ax3, betweenness, "Betweeness") -apply(g, None, weight, convex_hull, ax4, eigenvector, "Eigenvector") +apply(g, None, weight, convex_hull, ax1, betweenness, "Betweeness") -ax1, ax2, ax3, ax4 = row2 -apply(g, None, weight, convex_hull, ax1, katz, "Katz") -apply(g, None, weight, convex_hull, ax2, hits, "Hits") -apply(g, None, weight, convex_hull, ax3, leverage, "Leverage") -apply(g, None, weight, convex_hull, ax4, degree, "Degree") +# ax1, ax2, ax3, ax4 = row2 +# apply(g, None, weight, convex_hull, ax1, katz, "Katz") +# apply(g, None, weight, convex_hull, ax2, hits, "Hits") +# apply(g, None, weight, convex_hull, ax3, leverage, "Leverage") +# apply(g, None, weight, convex_hull, ax4, degree, "Degree") fig.savefig(f"Comparison_edge_centralities_artificial_.svg", format='svg') diff --git a/diff_comparison.py b/diff_comparison.py index 2c480e3..f105f3c 100644 --- a/diff_comparison.py +++ b/diff_comparison.py @@ -1,6 +1,7 @@ import math import matplotlib.pyplot as plt +from matplotlib.colors import TwoSlopeNorm import numpy as np from graph_tool.all import * @@ -87,7 +88,7 @@ def plot_graph_diff(G, c, fig, ax, name, cmap=plt.cm.plasma): x.append(ver[0]) y.append(ver[1]) - sc = ax.scatter(x, y, s=1, cmap=cmap, c=c) # map closeness values as color mapping on the verticies + sc = ax.scatter(x, y, s=1, cmap=cmap, norm=TwoSlopeNorm(0), c=c) # map closeness values as color mapping on the verticies ax.set_title(name) fig.colorbar(sc, ax=ax) @@ -95,8 +96,8 @@ def plot_graph_diff(G, c, fig, ax, name, cmap=plt.cm.plasma): def apply(g, seed, weight, convex_hull, ax, method, method_name): # calculate centrality values vp = None - if method_name == "Betweenness": - vp, ep = method(g, weight=weight) + if method_name == "Closeness": + vp = method(g, weight=weight) elif method_name == "Eigenvector": ep, vp = method(g, weight=weight) elif method_name == "Hits": @@ -136,8 +137,8 @@ def apply_corrected(g, seed, weight, convex_hull, ax, method, method_name): # calculate centrality values vp = None ep = None - if method_name == "Betweenness": - vp, ep = method(g, weight=weight) + if method_name == "Closeness": + vp = method(g, weight=weight) elif method_name == "Eigenvector": ep, vp = method(g, weight=weight) elif method_name == "Hits": @@ -183,7 +184,7 @@ def apply_corrected(g, seed, weight, convex_hull, ax, method, method_name): # - apply centrality measure to the next axis # - Draw the corresponding resulting models into a grid # -points, seed = random_graph(n=5000) +points, seed = random_graph(n=5000, seed=303437129487698362622376224319354280305) g, weight = spatial_graph(points) g = GraphView(g) @@ -193,15 +194,15 @@ convex_hull = centrality.convex_hull(g) # plot graph with convex_hull fig_graph, ax_graph = plt.subplots(figsize=(15, 12)) # draw without any centrality measure `vp` -vp, ep = betweenness(g, weight=weight) +vp = closeness(g, weight=weight) plot.graph_plot(fig_graph, ax_graph, g, vp, convex_hull, f"Pointcloud (seed: {seed})") -fig_graph.savefig("model_prediction_graph_original_betweenness_5000.svg", format='svg') +fig_graph.savefig("model_prediction_graph_original_closeness_5000.svg", format='svg') # normalization min_val, max_val = vp.a.min(), vp.a.max() vp.a = (vp.a - min_val) / (max_val - min_val) -vp_betweenness_original = vp +vp_closeness_original = vp for percentage in np.arange(0.1, 1, 0.1, dtype=float): print(f"Percentage: {percentage:.0%}") @@ -211,15 +212,15 @@ for percentage in np.arange(0.1, 1, 0.1, dtype=float): # draw subgraph fig_sub = plt.figure(figsize=(25, 12)) ax1, ax2 = fig_sub.subplots(1, 2) - vp, ep = betweenness(g_sub, weight=weight_sub) + vp = closeness(g_sub, weight=weight_sub) plot.graph_plot(fig_sub, ax1, g_sub, vp, convex_hull, f"{percentage:.0%} of Pointcloud (seed: {seed})") min_val, max_val = vp.a.min(), vp.a.max() vp.a = (vp.a - min_val) / (max_val - min_val) - vp_betweenness_corrected = apply_corrected(g_sub, seed, weight_sub, convex_hull, None, betweenness, "Betweenness") - plot.graph_plot(fig_sub, ax2, g_sub, vp_betweenness_corrected, convex_hull, f"{percentage:.0%} of Pointcloud with applied prediction") - fig_sub.savefig(f"model_prediction_subgraph_betweenness_5000_{percentage * 100:.0f}_percent.svg", format='svg') + vp_closeness_corrected = apply_corrected(g_sub, seed, weight_sub, convex_hull, None, closeness, "Closeness") + plot.graph_plot(fig_sub, ax2, g_sub, vp_closeness_corrected, convex_hull, f"{percentage:.0%} of Pointcloud with applied prediction") + fig_sub.savefig(f"model_prediction_subgraph_closeness_5000_{percentage * 100:.0f}_percent.svg", format='svg') distance_of_center = 0.5 * percentage @@ -243,10 +244,10 @@ for percentage in np.arange(0.1, 1, 0.1, dtype=float): # sub_position = g_sub.vp["pos"][sub_key] # print(f"position: {position} | sub_position: {sub_position}") - # calculate for betweenness - value = vp_betweenness_original[key] + # calculate for closeness + value = vp_closeness_original[key] pre_prediction = vp[sub_key] - sub_value = vp_betweenness_corrected[sub_key] + sub_value = vp_closeness_corrected[sub_key] scores.append(value) raw_sub_scores.append(pre_prediction) @@ -291,4 +292,4 @@ for percentage in np.arange(0.1, 1, 0.1, dtype=float): plot_graph_diff(g, diff_scores, fig, plot_sub_graph_ax, "Differences after correction of sub graph compared to original graph", plt.cm.seismic) plot_graph_diff(g, vp.a, fig, plot_sub_graph_before_ax, "Sub Graph (extracted region of original graph) without correction") - fig.savefig(f"model_prediction_subgraph_betweenness_5000_{percentage * 100:.0f}_percentage_diff.svg", format='svg') + fig.savefig(f"model_prediction_subgraph_closeness_5000_{percentage * 100:.0f}_percentage_diff.svg", format='svg') diff --git a/distance_types.py b/distance_types.py index d8218e2..c4b8081 100644 --- a/distance_types.py +++ b/distance_types.py @@ -6,7 +6,6 @@ from graph_tool.all import * from src import centrality from src import plot -from src import fitting def random_graph(n=5000, seed=None): @@ -40,19 +39,21 @@ def spatial_graph(adata): def apply(g, seed, weight, convex_hull, ax, ax2, method): # calculate centrality values vp, ep = method(g, weight=weight) - vp.a = np.nan_to_num(vp.a) # correct floating point values + ep.a = np.nan_to_num(ep.a) # correct floating point values + min_val, max_val = ep.a.min(), ep.a.max() + ep.a = (ep.a - min_val) / (max_val - min_val) # euklidian distance - quantification = plot.quantification_data(g, vp, convex_hull) + quantification = plot.quantification_data(g, ep, convex_hull) plot.quantification_plot(ax, quantification, None, None, "Euklidian Distance", None) # generate model based on convex hull and associated centrality values # path distance - quantification = plot.quantification_data_path_distance(g, weight, vp, convex_hull) + quantification = plot.quantification_data_path_distance(g, weight, ep, convex_hull) plot.quantification_plot(ax2, quantification, None, None, "Shortest Path Distance", None) -points, seed = random_graph(n=5000) +points, seed = random_graph(n=5000, seed=77011137629244244786159039169016117129) g, weight = spatial_graph(points) g = GraphView(g) # calculate convex hull @@ -62,12 +63,13 @@ fig = plt.figure(figsize=(21, 5)) ax1, ax2, ax3 = fig.subplots(1, 3) # plot graph with convex_hull -# draw without any centrality measure `vp` vp, ep = betweenness(g, weight=weight) -vp.a = np.nan_to_num(vp.a) # correct floating point values +ep.a = np.nan_to_num(ep.a) # correct floating point values +min_val, max_val = ep.a.min(), ep.a.max() +ep.a = (ep.a - min_val) / (max_val - min_val) -plot.graph_plot(fig, ax1, g, vp, convex_hull, f"Pointcloud (seed: {seed})") +plot.graph_plot(fig, ax1, g, ep, convex_hull, f"Pointcloud (seed: {seed})", True) apply(g, seed, weight, convex_hull, ax2, ax3, betweenness) -fig.savefig(f"Distance_5000_betweenness_euklidian.svg", format='svg') +fig.savefig(f"Distance_5000_betweenness_edge_euklidian.svg", format='svg') diff --git a/src/fitting.py b/src/fitting.py index 1cabd10..f99e041 100644 --- a/src/fitting.py +++ b/src/fitting.py @@ -52,6 +52,8 @@ def fit_piece_wise_linear(d, C, M=1000): model.addConstr((1 - z[i]) * M >= d[i] - b) model.optimize() + + # FIXME does not work with real world data? what am I doing wrong? # AIC k = 4 aic = 2. * k + n * math.log(model.ObjVal) diff --git a/src/plot.py b/src/plot.py index 70e8a06..49cab6f 100644 --- a/src/plot.py +++ b/src/plot.py @@ -60,16 +60,18 @@ def graph_plot(fig, ax, G, measures, convex_hull, name, show_edges=False): c = measures.get_array() # convex hull -> Bounding-Box ch = LineCollection([convex_hull], colors=['g'], linewidths=1) + ax.set_title(name) ax.add_collection(ch) if show_edges: - for e in G.edges(): + cmap = plt.cm.plasma.resampled(G.num_edges()) + for idx, e in enumerate(G.edges()): ex = [pos[e.source()][0], pos[e.target()][0]] ey = [pos[e.source()][1], pos[e.target()][1]] - ax.add_collection(LineCollection([np.column_stack([ex, ey])], colors=['k'], linewidths=0.1)) - - sc = ax.scatter(x, y, s=1, cmap=plt.cm.plasma, c=c) # map closeness values as color mapping on the verticies - ax.set_title(name) - fig.colorbar(sc, ax=ax) + ax.add_collection(LineCollection([np.column_stack([ex, ey])], colors=cmap(c[idx]), linewidths=0.5)) + fig.colorbar(plt.cm.ScalarMappable(norm=None, cmap=cmap), ax=ax) + else: + sc = ax.scatter(x, y, s=1, c=c, cmap=plt.cm.plasma) # map closeness values as color mapping on the verticies + fig.colorbar(sc, ax=ax) def graph_plot_effected(fig, ax, G, measures, convex_hull, b, name, show_edges=False): @@ -145,6 +147,10 @@ def quantification_data(G, measures, convex_hull): min_distance = math.inf key = next(keys) for edge in convex_hull: + # TODO isn't there the dot product missing? + # -> such that there might be a shorter path? + # -> for each `point` take its each of its two neighbours (idx - 1 & idx + 1) + # and create another vector on which you project the verticies too? vector = Vector.vec(point, edge) distance = Vector.vec_len(vector) if distance < min_distance: @@ -156,6 +162,50 @@ def quantification_data(G, measures, convex_hull): return np.array(quantification) +def quantification_data_edges(G, measures, convex_hull): + # calculate distance based on the median of the distances of the two verticies an edge connects + quantification = [] + pos = G.vp["pos"] + x = [] + y = [] + for v in G.vertices(): + ver = pos[v] + x.append(ver[0]) + y.append(ver[1]) + + measures = measures.a + keys = iter(measures) + + points = np.stack((np.array(x), np.array(y)), axis=-1) + for e in G.edges(): + min_distance_source = math.inf + min_distance_target = math.inf + key = next(keys) + for point in convex_hull: + # TODO isn't there the dot product missing? + # -> such that there might be a shorter path? + # -> for each `point` take its each of its two neighbours (idx - 1 & idx + 1) + # and create another vector on which you project the verticies too? + vector = Vector.vec(pos[e.source()], point) + distance = Vector.vec_len(vector) + if distance < min_distance_source: + min_distance_source = distance + for point in convex_hull: + # TODO isn't there the dot product missing? + # -> such that there might be a shorter path? + # -> for each `point` take its each of its two neighbours (idx - 1 & idx + 1) + # and create another vector on which you project the verticies too? + vector = Vector.vec(pos[e.target()], point) + distance = Vector.vec_len(vector) + if distance < min_distance_target: + min_distance_target = distance + quantification.append([(min_distance_target + min_distance_source) / 2, key]) + + # sort by distance + quantification.sort(key=lambda entry: entry[0]) + return np.array(quantification) + + def quantification_data_path_distance(G, weights, measures, convex_hull): quantification = [] pos = G.vp["pos"] @@ -172,11 +222,26 @@ def quantification_data_path_distance(G, weights, measures, convex_hull): keys = iter(measures) points = np.stack((np.array(x), np.array(y)), axis=-1) - for v in G.vertices(): + for idx, e in enumerate(G.edges()): + ex = [pos[e.source()][0], pos[e.target()][0]] + ey = [pos[e.source()][1], pos[e.target()][1]] min_distance = math.inf key = next(keys) + # short cut + if e.source() in convex_hull_verticies or e.target() in convex_hull_verticies: + quantification.append([0, key]) + continue + + # for either side of the edge (source) for h in convex_hull_verticies: - vertices, edges = graph_tool.topology.shortest_path(G, v, h, weights=weights) + vertices, edges = graph_tool.topology.shortest_path(G, e.source(), h, weights=weights) + # TODO calculate the total distance + path_length = sum([weights[edge] for edge in edges]) + if path_length < min_distance: + min_distance = path_length + # for either side of the edge (target) + for h in convex_hull_verticies: + vertices, edges = graph_tool.topology.shortest_path(G, e.target(), h, weights=weights) # TODO calculate the total distance path_length = sum([weights[edge] for edge in edges]) if path_length < min_distance: