mod update corresponding examples

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)
+# points, seed = random_graph(n=5000)
-g, weight = spatial_graph(points)
+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})")
+plot.graph_plot(fig_graph, ax_graph, g, vp, convex_hull, f"mibitof")
-fig_graph.savefig("Pointcloud_graph.svg", format='svg')
+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, None, weight, convex_hull, ax1, closeness, "Closeness")
-apply(g, seed, weight, convex_hull, ax2, pagerank, "PageRank")
+apply(g, None, weight, convex_hull, ax2, pagerank, "PageRank")
-apply(g, seed, weight, convex_hull, ax3, betweenness, "Betweeness")
+apply(g, None, weight, convex_hull, ax3, betweenness, "Betweeness")
-apply(g, seed, weight, convex_hull, ax4, eigenvector, "Eigenvector")
+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, None, weight, convex_hull, ax1, katz, "Katz")
-apply(g, seed, weight, convex_hull, ax2, hits, "Hits")
+apply(g, None, weight, convex_hull, ax2, hits, "Hits")
-apply(g, seed, weight, convex_hull, ax3, leverage, "Leverage")
+apply(g, None, weight, convex_hull, ax3, leverage, "Leverage")
-apply(g, seed, weight, convex_hull, ax4, degree, "Degree")
+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')

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)
+        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
 
     # normalization
-    min_val, max_val = vp.a.min(), vp.a.max()
+    min_val, max_val = ep.a.min(), ep.a.max()
-    vp.a = (vp.a - min_val) / (max_val - min_val)
+    ep.a = (ep.a - min_val) / (max_val - min_val)
 
-    # generate model based on convex hull and associated centrality values
+    quantification = plot.quantification_data_edges(g, ep, convex_hull)
-    quantification = plot.quantification_data(g, vp, 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`
+# draw without any centrality measure `ep`
-vp = g.new_vertex_property("double")
+ep = g.new_edge_property("double")
-plot.graph_plot(fig_graph, ax_graph, g, vp, convex_hull, f"Pointcould (seed: {seed})")
+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, ax1, betweenness, "Betweeness")
-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
+# ax1, ax2, ax3, ax4 = row2
-apply(g, None, weight, convex_hull, ax1, katz, "Katz")
+# apply(g, None, weight, convex_hull, ax1, katz, "Katz")
-apply(g, None, weight, convex_hull, ax2, hits, "Hits")
+# apply(g, None, weight, convex_hull, ax2, hits, "Hits")
-apply(g, None, weight, convex_hull, ax3, leverage, "Leverage")
+# apply(g, None, weight, convex_hull, ax3, leverage, "Leverage")
-apply(g, None, weight, convex_hull, ax4, degree, "Degree")
+# apply(g, None, weight, convex_hull, ax4, degree, "Degree")
 
 fig.savefig(f"Comparison_edge_centralities_artificial_.svg", format='svg')

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":
+    if method_name == "Closeness":
-        vp, ep = method(g, weight=weight)
+        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":
+    if method_name == "Closeness":
-        vp, ep = method(g, weight=weight)
+        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")
+    vp_closeness_corrected = apply_corrected(g_sub, seed, weight_sub, convex_hull, None, closeness, "Closeness")
-    plot.graph_plot(fig_sub, ax2, g_sub, vp_betweenness_corrected, convex_hull, f"{percentage:.0%} of Pointcloud with applied prediction")
+    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_betweenness_5000_{percentage * 100:.0f}_percent.svg", format='svg')
+    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
+        # calculate for closeness
-        value = vp_betweenness_original[key]
+        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')

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

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)

src/plot.py

@@ -60,15 +60,17 @@ 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))
+            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)
-    sc = ax.scatter(x, y, s=1, cmap=plt.cm.plasma, c=c) # map closeness values as color mapping on the verticies
+    else:
-    ax.set_title(name)
+        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)
 
 
@@ -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: