From 3acf54a000854277c64791efdff57f94ee2b3129 Mon Sep 17 00:00:00 2001 From: Yves Biener Date: Thu, 16 Apr 2026 07:20:19 +0200 Subject: [PATCH] WIP --- comparison.py | 119 +++++++++++++++++++++++++++++++------- comparison_edge_scores.py | 20 +++++++ distance_types.py | 8 +-- point_cloud_example.py | 35 +++++++++-- src/plot.py | 23 +++----- 5 files changed, 161 insertions(+), 44 deletions(-) diff --git a/comparison.py b/comparison.py index 130b59d..62b7f8b 100644 --- a/comparison.py +++ b/comparison.py @@ -1,8 +1,10 @@ import math import matplotlib.pyplot as plt +from matplotlib.collections import LineCollection import numpy as np import squidpy as sq +import scipy from graph_tool.all import * from src import centrality @@ -56,6 +58,25 @@ def leverage(g, weight): return vp +def laplacian(g, weight): + vp = g.new_vertex_property("double") + lap_g = graph_tool.spectral.laplacian(g, weight=weight) + elap_g = sum(l**2 for l in scipy.linalg.eigvals(lap_g.toarray())) + for v in g.vertices(): + gv = g.copy() + gv.remove_vertex(v, True) + # pos = gv.vp["pos"] + + # weight_gv = gv.new_edge_property("double") + # for e in gv.edges(): + # weight_gv[e] = math.sqrt(sum(map(abs, pos[e.source()].a - pos[e.target()].a)))**2 + + lap_gv = graph_tool.spectral.laplacian(gv, weight=gv.ep["weight"]) + elap_gv = sum(l**2 for l in scipy.linalg.eigvals(lap_gv.toarray())) + vp[v] = (elap_g - elap_gv) / elap_g + return vp + + def random_graph(n=5000, seed=None): """ Uniformly random point cloud generation. @@ -81,6 +102,7 @@ def spatial_graph(adata): weight = g.new_edge_property("double") for e in g.edges(): weight[e] = math.sqrt(sum(map(abs, pos[e.source()].a - pos[e.target()].a)))**2 + g.ep["weight"] = weight return g, weight @@ -121,6 +143,25 @@ def apply(g, seed, weight, convex_hull, ax, method, method_name): plot.quantification_plot(ax, quantification, d_curve, C_curve, method_name, aic_opt) +def draw_graph(G, ax, name): + pos = G.vp["pos"] + x = [] + y = [] + for v in G.vertices(): + ver = pos[v] + x.append(ver[0]) + y.append(ver[1]) + + # edges + for e in 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)) + + ax.scatter(x, y, s=1) + ax.set_title(name) + + # # - Create a random point cloud and calculate a triangulation on it # - For that graph calculate the convex hull @@ -129,34 +170,68 @@ 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) -adata = mibitof() -g, weight = spatial_graph(adata.obsm['spatial']) +points, seed = random_graph(n=3000) +# adata = merfish() +# g, weight = spatial_graph(adata.obsm['spatial']) +g, weight = spatial_graph(points) g = GraphView(g) + +# NOTE remove duplicated node that has is an isolated node +# only relevant for `mibitof` +# for v in g.vertices(): +# neighbours = g.get_all_neighbours(v) +# if len(neighbours) == 0: +# g.remove_vertex(v) +# break + +# pos = g.vp["pos"] +# weight = g.new_edge_property("double") +# for e in g.edges(): +# weight[e] = math.sqrt(sum(map(abs, pos[e.source()].a - pos[e.target()].a)))**2 + # calculate convex hull convex_hull = centrality.convex_hull(g) -# plot graph with convex_hull +# plot graph 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"mibitof") -fig_graph.savefig(f"mibitof_graph.svg", format='svg') +draw_graph(g, ax_graph, f"Artifical (n=3000)\n(seed = {seed})") +fig_graph.savefig(f"Comparison_node_artificial_3000_graph.svg", format='svg') -fig = plt.figure(figsize=(15, 12)) -row1, row2 = fig.subplots(2, 4) +# | Closeness | PageRank | Eigenvector | Leverage | +# | Betweenness | Katz | Laplacian | Degree | +# | | Hits | | | +fig, ax = plt.subplots(figsize=(15, 12)) +apply(g, None, weight, convex_hull, ax, closeness, "Closeness") +fig.savefig(f"Comparison_node_closeness_artifical_3000.svg", format='svg') -ax1, ax2, ax3, ax4 = row1 -# TODO select corresponding centrality measure method -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") +fig, ax = plt.subplots(figsize=(15, 12)) +apply(g, None, weight, convex_hull, ax, betweenness, "Betweeness") +fig.savefig(f"Comparison_node_betweenness_artifical_3000.svg", format='svg') -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, ax = plt.subplots(figsize=(15, 12)) +apply(g, None, weight, convex_hull, ax, pagerank, "PageRank") +fig.savefig(f"Comparison_node_pagerank_artifical_3000.svg", format='svg') -fig.savefig(f"Comparison_node_centralities_mibitof_.svg", format='svg') +fig, ax = plt.subplots(figsize=(15, 12)) +apply(g, None, weight, convex_hull, ax, eigenvector, "Eigenvector") +fig.savefig(f"Comparison_node_eigenvector_artifical_3000.svg", format='svg') + +fig, ax = plt.subplots(figsize=(15, 12)) +apply(g, None, weight, convex_hull, ax, hits, "Hits") +fig.savefig(f"Comparison_node_hits_artifical_3000.svg", format='svg') + +fig, ax = plt.subplots(figsize=(15, 12)) +apply(g, None, weight, convex_hull, ax, katz, "Katz") +fig.savefig(f"Comparison_node_katz_artifical_3000.svg", format='svg') + +fig, ax = plt.subplots(figsize=(15, 12)) +apply(g, None, weight, convex_hull, ax, degree, "Degree") +fig.savefig(f"Comparison_node_degree_artifical_3000.svg", format='svg') + +fig, ax = plt.subplots(figsize=(15, 12)) +apply(g, None, weight, convex_hull, ax, leverage, "Leverage") +fig.savefig(f"Comparison_node_leverage_artifical_3000.svg", format='svg') + +fig, ax = plt.subplots(figsize=(15, 12)) +apply(g, None, weight, convex_hull, ax, laplacian, "Laplacian") +fig.savefig(f"Comparison_node_laplacian_artifical_3000.svg", format='svg') diff --git a/comparison_edge_scores.py b/comparison_edge_scores.py index 2d20cf5..86f5e0e 100644 --- a/comparison_edge_scores.py +++ b/comparison_edge_scores.py @@ -55,6 +55,26 @@ def leverage(g, weight): return vp +def path(g, weight): + # NOTE is this not just betweenness? + ep = g.new_vertex_property("double") + for v in g.vertices(): + for u in g.vertices(): + if (v == u): + continue + paths = graph_tool.topology.all_shortest_paths(g, v, u, weights=weight, edges=True) + for edges in paths: + for edge in edges: + for idx, g_e in enumerate(g.edges()): + if (g_e == edge): + # NOTE we end up counting twice! + ep[idx] += 0.5; + break + # for e in g.edges(): + # ep[e] /= 2; + return ep + + def random_graph(n=5000, seed=None): """ Uniformly random point cloud generation. diff --git a/distance_types.py b/distance_types.py index 2e46ffb..d36a0bc 100644 --- a/distance_types.py +++ b/distance_types.py @@ -38,7 +38,7 @@ def spatial_graph(adata): def apply(g, seed, weight, convex_hull, ax, ax2, method): # calculate centrality values - vp = method(g, weight=weight) + vp, ep = method(g, weight=weight) vp.a = np.nan_to_num(vp.a) # correct floating point values min_val, max_val = vp.a.min(), vp.a.max() vp.a = (vp.a - min_val) / (max_val - min_val) @@ -63,13 +63,13 @@ fig = plt.figure(figsize=(21, 5)) ax1, ax2, ax3 = fig.subplots(1, 3) # plot graph with convex_hull -vp = closeness(g, weight=weight) +vp, ep = betweenness(g, weight=weight) vp.a = np.nan_to_num(vp.a) # correct floating point values min_val, max_val = vp.a.min(), vp.a.max() vp.a = (vp.a - min_val) / (max_val - min_val) plot.graph_plot(fig, ax1, g, vp, convex_hull, f"Pointcloud (seed: {seed})", False) -apply(g, seed, weight, convex_hull, ax2, ax3, closeness) +apply(g, seed, weight, convex_hull, ax2, ax3, betweenness) -fig.savefig(f"Distance_5000_node_closeness.svg", format='svg') +fig.savefig(f"Distance_5000_node_betweenness.svg", format='svg') diff --git a/point_cloud_example.py b/point_cloud_example.py index aa78c90..8dd1c66 100644 --- a/point_cloud_example.py +++ b/point_cloud_example.py @@ -3,6 +3,7 @@ import math import matplotlib.pyplot as plt from matplotlib.collections import LineCollection import numpy as np +import squidpy as sq from graph_tool.all import * @@ -19,6 +20,14 @@ def random_graph(n=5000, seed=None): return rng.random((n, 2)), seed +def mibitof(): + """ + Mibitof dataset from `squidpy`. + """ + adata = sq.datasets.mibitof() + return adata + + def spatial_graph(adata): """ Generate the spatial graph using delaunay for the given `adata`. @@ -55,6 +64,11 @@ def draw_graph(G, ax, name): ax.add_collection(LineCollection([np.column_stack([ex, ey])], colors=['k'], linewidths=0.1)) ax.scatter(x, y, s=1) # map closeness values as color mapping on the verticies + # for v in G.vertices(): + # neighbours = g.get_all_neighbours(v) + # if len(neighbours) == 0: + # ax.scatter(pos[v][0], pos[v][1], s=1, color=['r']) + ax.set_title(name) @@ -66,11 +80,24 @@ def draw_graph(G, ax, name): # - apply centrality measure to the next axis # - Draw the corresponding resulting models into a grid # -points, seed = random_graph(n=3000) -g, weight = spatial_graph(points) +# points, seed = random_graph(n=3000) +# g, weight = spatial_graph(points) +adata = mibitof() +g, weight = spatial_graph(adata.obsm['spatial']) g = GraphView(g) +for v in g.vertices(): + neighbours = g.get_all_neighbours(v) + if len(neighbours) == 0: + g.remove_vertex(v) + break + +pos = g.vp["pos"] +weight = g.new_edge_property("double") +for e in g.edges(): + weight[e] = math.sqrt(sum(map(abs, pos[e.source()].a - pos[e.target()].a)))**2 + # plot graph with convex_hull fig_graph, ax_graph = plt.subplots(figsize=(15, 12)) -draw_graph(g, ax_graph, f"Pointcould (seed: {seed} | n: 500)") -fig_graph.savefig("point_cloud_example.svg", format='svg') +draw_graph(g, ax_graph, f"mibitof") +fig_graph.savefig("mibitof_graph.svg", format='svg') diff --git a/src/plot.py b/src/plot.py index 092e600..84e9065 100644 --- a/src/plot.py +++ b/src/plot.py @@ -165,8 +165,7 @@ def quantification_data(G, measures, convex_hull): def quantification_data_node_path_distance(G, weights, measures, convex_hull): quantification = [] pos = G.vp["pos"] - x = [] - y = [] + convex_hull_verticies = [] for v in G.vertices(): ver = pos[v] @@ -214,18 +213,16 @@ def quantification_data_edges(G, measures, convex_hull): min_distance_source = math.inf min_distance_target = math.inf key = next(keys) - for idx, point in enumerate(convex_hull): - hull_line = Vector.vec(convex_hull[idx - 1], point) - a = point - b = convex_hull[idx - 1] - distance = abs((a[1] - b[1]) * pos[e.source()][0] - (a[0] - b[0]) * pos[e.source()][1] + a[1]*b[0] - b[1]*a[0])/Vector.vec_len(hull_line) + for point in convex_hull: + v = np.stack((np.array([pos[e.source()][0]]), np.array([pos[e.source()][1]])), axis=-1) + vector = Vector.vec(v[0], edge) + distance = Vector.vec_len(vector) if distance < min_distance_source: min_distance_source = distance for point in convex_hull: - hull_line = Vector.vec(convex_hull[idx - 1], point) - a = point - b = convex_hull[idx - 1] - distance = abs((a[1] - b[1]) * pos[e.target()][0] - (a[0] - b[0]) * pos[e.target()][1] + a[1]*b[0] - b[1]*a[0])/Vector.vec_len(hull_line) + v = np.stack((np.array([pos[e.target()][0]]), np.array([pos[e.target()][1]])), axis=-1) + vector = Vector.vec(v[0], edge) + distance = Vector.vec_len(vector) if distance < min_distance_target: min_distance_target = distance quantification.append([(min_distance_target + min_distance_source) / 2, key]) @@ -238,8 +235,7 @@ def quantification_data_edges(G, measures, convex_hull): def quantification_data_path_distance(G, weights, measures, convex_hull): quantification = [] pos = G.vp["pos"] - x = [] - y = [] + convex_hull_verticies = [] for v in G.vertices(): ver = pos[v] @@ -250,7 +246,6 @@ def quantification_data_path_distance(G, weights, measures, convex_hull): measures = measures.a keys = iter(measures) - points = np.stack((np.array(x), np.array(y)), axis=-1) 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]]