WIP

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')

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.

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')

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')

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]]