update scripts

Author: bdpedigo

experiments/evaluate_clustering/evaluate_clustering.py

@@ -60,7 +60,7 @@ def stashcsv(df, name, **kws):
     savecsv(df, name, pathname=save_path / "outs", **kws)
 
 
-def sort_mg(mg, level_names, class_order=CLASS_ORDER):
+def sort_mg(mg, level_names, class_order=CLASS_ORDER, ascending=True):
     """Required sorting prior to plotting the dendrograms
 
     Parameters
@@ -85,7 +85,7 @@ def sort_mg(mg, level_names, class_order=CLASS_ORDER):
             meta[f"{sc}_{co}_order"] = meta[sc].map(class_value)
             total_sort_by.append(f"{sc}_{co}_order")
         total_sort_by.append(sc)
-    mg = mg.sort_values(total_sort_by, ascending=False)
+    mg = mg.sort_values(total_sort_by, ascending=ascending)  # TODO used to be False!
     return mg
 
 
@@ -96,13 +96,14 @@ def plot_adjacencies(full_mg, axs, lowest_level=7):
     for level in np.arange(lowest_level + 1):
         ax = axs[0, level]
         adj = binarize(full_mg.adj)
+        # [f"lvl{level}_labels", f"merge_class_sf_order", "merge_class"]
         _, _, top, _ = adjplot(
             adj,
             ax=ax,
             plot_type="scattermap",
             sizes=(0.5, 0.5),
             sort_class=level_names[: level + 1],
-            item_order=[f"{CLASS_KEY}_{CLASS_ORDER}_order", CLASS_KEY, CLASS_ORDER],
+            # item_order=[f"{CLASS_KEY}_{CLASS_ORDER}_order", CLASS_KEY, CLASS_ORDER],
             class_order=CLASS_ORDER,
             meta=full_mg.meta,
             palette=CLASS_COLOR_DICT,
@@ -495,37 +496,106 @@ def plot_clustering_results(
 n_init = 256
 max_hops = 16
 allow_loops = False
+include_reverse = False
 walk_spec = f"gt={graph_type}-n_init={n_init}-hops={max_hops}-loops={allow_loops}"
 walk_meta = pd.read_csv(
-    f"maggot_models/experiments/walk_sort/outs/meta_w_order-{walk_spec}.csv",
+    f"maggot_models/experiments/walk_sort/outs/meta_w_order-{walk_spec}-include_reverse={include_reverse}.csv",
     index_col=0,
 )
 meta["median_node_visits"] = walk_meta["median_node_visits"]
 
 # %%
 # plot results
 lowest_level = 7  # last level to show for dendrograms, adjacencies
-plot_clustering_results(
-    adj,
-    meta,
-    basename,
-    lowest_level=lowest_level,
-    show_adjs=True,
-    show_singles=False,
-    make_flippable=False,
-)
+# plot_clustering_results(
+#     adj,
+#     meta,
+#     basename,
+#     lowest_level=lowest_level,
+#     show_adjs=True,
+#     show_singles=False,
+#     make_flippable=False,
+# )
 
 #%%
-# lowest_level = 7
-# mg = MetaGraph(adj, meta)
-# level_names = [f"lvl{i}_labels" for i in range(lowest_level + 1)]
-# mg = sort_mg(mg, level_names)
-# fig, axs = plt.subplots(
-#     2, lowest_level + 1, figsize=10 * np.array([lowest_level + 1, 2])
-# )
+lowest_level = 7
+mg = MetaGraph(adj, meta)
+level_names = [f"lvl{i}_labels" for i in range(lowest_level + 1)]
+mg = sort_mg(mg, level_names)
+fig, axs = plt.subplots(
+    2, lowest_level + 1, figsize=10 * np.array([lowest_level + 1, 2])
+)
 # for level in np.arange(lowest_level + 1):
-#     plot_adjacencies(mg, axs, lowest_level=lowest_level)
-# stashfig(f"adjplots-lowest={lowest_level}" + basename, fmt="png")
+plot_adjacencies(mg, axs, lowest_level=lowest_level)
+stashfig(f"adjplots-lowest={lowest_level}" + basename, fmt="png")
+#%%
+from matplotlib.colors import ListedColormap
+
+sort_meta = mg.meta.copy()
+fig, axs = plt.subplots(
+    1, 2 * (lowest_level + 1), figsize=(10, 10), gridspec_kw=dict(wspace=0)
+)
+
+# meta = mg.meta
+# sort_class = level_names + ["merge_class"]
+# class_order = [class_order]
+# total_sort_by = []
+# for sc in sort_class:
+#     for co in class_order:
+#         class_value = meta.groupby(sc)[co].mean()
+#         meta[f"{sc}_{co}_order"] = meta[sc].map(class_value)
+#         total_sort_by.append(f"{sc}_{co}_order")
+#     total_sort_by.append(sc)
+# mg = mg.sort_values(total_sort_by, ascending=False)
+
+
+for level in np.arange(lowest_level + 1)[::-1]:
+    # sort_meta = sort_meta.sort_values(
+    #     [
+    #         f"lvl{level}_labels_{CLASS_ORDER}_order",
+    #         f"lvl{level}_labels",
+    #         f"{CLASS_KEY}_{CLASS_ORDER}_order",
+    #         CLASS_KEY,
+    #     ],
+    #     ascending=True,
+    # )
+    sort_meta["inds"] = range(len(sort_meta))
+    firsts = sort_meta.groupby(f"lvl{level}_labels", sort=False)["inds"].first()
+
+    # mean_visits = sort_meta.groupby(
+    #     [
+    #         f"lvl{level}_labels",
+    #         f"{CLASS_KEY}_{CLASS_ORDER}_order",
+    #     ]
+    # )["median_node_visit"].mean()
+    # meta.groupby([leaf_key, "merge_class"], sort=False).size()
+
+    sort_meta[CLASS_KEY].values
+    color_dict = CLASS_COLOR_DICT
+    classes = sort_meta["merge_class"].values
+    uni_classes = np.unique(sort_meta["merge_class"])
+    class_map = dict(zip(uni_classes, range(len(uni_classes))))
+    color_sorted = np.vectorize(color_dict.get)(uni_classes)
+    lc = ListedColormap(color_sorted)
+    class_indicator = np.vectorize(class_map.get)(classes)
+    class_indicator = class_indicator.reshape(len(classes), 1)
+    ax = axs[2 * level + 1]
+    sns.heatmap(
+        class_indicator,
+        cmap=lc,
+        cbar=False,
+        yticklabels=False,
+        # xticklabels=False,
+        square=False,
+        ax=ax,
+    )
+    ax.set(xlabel=level, xticks=[])
+
+    ax = axs[2 * level]
+    ax.axis("off")
+    ax.set(ylim=axs[2 * level + 1].get_ylim())
+    for first_ind in firsts:
+        ax.axhline(first_ind, color="grey", linestyle="--", alpha=1, linewidth=1)
 
 # %% [markdown]
 # # ##

experiments/nblast/nblast.py

@@ -0,0 +1,122 @@
+#%%
+import logging
+import time
+
+import numpy as np
+import pandas as pd
+import pymaid
+from sklearn.preprocessing import QuantileTransformer
+from pathlib import Path
+
+from graspologic.utils import symmetrize
+from navis import NeuronList, TreeNeuron, nblast_allbyall
+from src.data import load_metagraph
+from src.pymaid import start_instance
+
+# REF: https://stackoverflow.com/questions/35326814/change-level-logged-to-ipython-jupyter-notebook
+# logger = logging.getLogger()
+# # assert len(logger.handlers) == 1
+# handler = logger.handlers[0]
+# handler.setLevel(logging.ERROR)
+
+t0 = time.time()
+
+# for pymaid to pull neurons
+start_instance()
+
+out_dir = Path("maggot_models/experiments/nblast/outs")
+
+#%% load connectivity data
+mg = load_metagraph("G")
+meta = mg.meta
+
+#%% define some functions
+
+
+def pairwise_nblast(neuron_ids, point_thresh=5):
+    neuron_ids = [int(n) for n in neuron_ids]
+    neurons = pymaid.get_neuron(neuron_ids)  # load in with pymaid
+
+    # HACK: I am guessing there is a better way to do the below?
+    # TODO: I was also getting some errors about neurons with more that one soma, so I
+    # threw them out for now.
+    treenode_tables = []
+    for neuron_id, neuron in zip(neuron_ids, neurons):
+        treenode_table = pymaid.get_treenode_table(neuron, include_details=False)
+        treenode_tables.append(treenode_table)
+
+    success_neurons = []
+    tree_neurons = []
+    for neuron_id, treenode_table in zip(neuron_ids, treenode_tables):
+        treenode_table.rename(columns={"parent_node_id": "parent_id"}, inplace=True)
+
+        tree_neuron = TreeNeuron(treenode_table)
+        if (tree_neuron.soma is not None) and (len(tree_neuron.soma) > 1):
+            print(f"Neuron {neuron_id} has more than one soma, removing")
+        elif len(treenode_table) < point_thresh:
+            print(f"Neuron {neuron_id} has fewer than {point_thresh} points, removing")
+        else:
+            tree_neurons.append(tree_neuron)
+            success_neurons.append(neuron_id)
+
+    tree_neurons = NeuronList(tree_neurons)
+    print(f"{len(tree_neurons)} neurons ready for NBLAST")
+
+    currtime = time.time()
+    # NOTE: I've had too modify original code to allow smat=None
+    # NOTE: this only works when normalized=False also
+    scores = nblast_allbyall(tree_neurons, smat=None, normalized=False, progress=True)
+    print(f"{time.time() - currtime:.3f} elapsed to run NBLAST.")
+
+    scores = pd.DataFrame(
+        data=scores.values, index=success_neurons, columns=success_neurons
+    )
+
+    return scores
+
+
+def postprocess_nblast(scores):
+    distance = scores.values  # the raw nblast scores are dissimilarities/distances
+    sym_distance = symmetrize(distance)  # the raw scores are not symmetric
+    # make the distances between 0 and 1
+    sym_distance /= sym_distance.max()
+    sym_distance -= sym_distance.min()
+    # and then convert to similarity
+    morph_sim = 1 - sym_distance
+
+    # rank transform the similarities
+    # NOTE this is very different from what native NBLAST does and could likely be
+    # improved upon a lot. I did this becuase it seemed like a quick way of accounting
+    # for difference in scale for different neurons as well as the fact that the raw
+    # distribution of similaritys was skewed low (very few small values)
+    quant = QuantileTransformer()
+    indices = np.triu_indices_from(morph_sim, k=1)
+    transformed_vals = quant.fit_transform(morph_sim[indices].reshape(-1, 1))
+    transformed_vals = np.squeeze(transformed_vals)
+    # this is a discrete version of PTR basically
+    ptr_morph_sim = np.ones_like(morph_sim)
+    ptr_morph_sim[indices] = transformed_vals
+    ptr_morph_sim[indices[::-1]] = transformed_vals
+
+    ptr_morph_sim = pd.DataFrame(
+        data=ptr_morph_sim, index=scores.index, columns=scores.columns
+    )
+
+    return ptr_morph_sim
+
+
+#%% run nblast
+for side in ["left", "right"]:
+    print(f"Processing side: {side}")
+    side_meta = meta[meta[side]]
+
+    scores = pairwise_nblast(side_meta.index.values)
+    scores.to_csv(out_dir / f"{side}-nblast-scores.csv")
+
+    similarity = postprocess_nblast(scores)
+    similarity.to_csv(out_dir / f"{side}-nblast-similarities.csv")
+    print()
+
+#%%
+print("\n\n")
+print(f"{time.time() - t0:.3f} elapsed for whole script.")

experiments/walk_sort/generate_walks.py

@@ -102,7 +102,7 @@ def main(
         ("dVNC",),
         ("dSEZ",),
         ("RGN",),
-        ("m`otor-PaN", "motor-MN", "motor-AN", "motor-VAN"),
+        ("motor-PaN", "motor-MN", "motor-AN", "motor-VAN"),
     ]
     stop_names = ["dVNC", "dSEZ", "RGN", "motor"]
 

experiments/walk_sort/walk_sort.py

@@ -54,6 +54,7 @@ def stashcsv(df, name, **kws):
 n_init = 256
 max_hops = 16
 allow_loops = False
+include_reverse = False
 walk_path = "maggot_models/experiments/walk_sort/outs/walks-"
 walk_spec = f"gt={graph_type}-n_init={n_init}-hops={max_hops}-loops={allow_loops}"
 forward_walk_path = walk_path + walk_spec + "-reverse=False" + ".txt"
@@ -71,15 +72,15 @@ def process_paths(walk_path):
     paths.remove("")
     print(f"# of paths after removing duplicates: {len(paths)}")
 
-    n_subsample = len(paths)  # 2 ** 14
-    choice_inds = np.random.choice(len(paths), n_subsample, replace=False)
-    new_paths = []
-    for i in range(len(paths)):
-        if i in choice_inds:
-            new_paths.append(paths[i])
-    paths = new_paths
+    # n_subsample = len(paths)  # 2 ** 14
+    # choice_inds = np.random.choice(len(paths), n_subsample, replace=False)
+    # new_paths = []
+    # for i in range(len(paths)):
+    #     if i in choice_inds:
+    #         new_paths.append(paths[i])
+    # paths = new_paths
 
-    print(f"# of paths after subsampling: {len(paths)}")
+    # print(f"# of paths after subsampling: {len(paths)}")
     paths = [path.split(" ") for path in paths]
     paths = [[int(node) for node in path] for path in paths]
     # all_nodes = set()
@@ -95,23 +96,23 @@ def process_paths(walk_path):
 
 
 # %%
-
+all_nodes = set()
 node_visits = {}
 for path in forward_paths:
     for i, node in enumerate(path):
         if node not in node_visits:
             node_visits[node] = []
         node_visits[node].append(i / (len(path) - 1))
+[[all_nodes.add(node) for node in path] for path in forward_paths]
 
-for path in backward_paths:
-    for i, node in enumerate(path):
-        if node not in node_visits:
-            node_visits[node] = []
-        node_visits[node].append(1 - (i / (len(path) - 1)))
+if include_reverse:
+    for path in backward_paths:
+        for i, node in enumerate(path):
+            if node not in node_visits:
+                node_visits[node] = []
+            node_visits[node].append(1 - (i / (len(path) - 1)))
+    [[all_nodes.add(node) for node in path] for path in backward_paths]
 
-all_nodes = set()
-[[all_nodes.add(node) for node in path] for path in forward_paths]
-[[all_nodes.add(node) for node in path] for path in backward_paths]
 uni_nodes = np.unique(list(all_nodes))
 
 median_node_visits = {}
@@ -128,7 +129,9 @@ def process_paths(walk_path):
     median_class_visits[node_class] = np.median(all_visits_flat)
 meta["median_class_visits"] = meta["merge_class"].map(median_class_visits)
 
-meta.to_csv(f"maggot_models/experiments/walk_sort/outs/meta_w_order-{walk_spec}.csv")
+meta.to_csv(
+    f"maggot_models/experiments/walk_sort/outs/meta_w_order-{walk_spec}-include_reverse={include_reverse}.csv"
+)
 
 print(f"# of nodes: {len(meta)}")
 unvisit_meta = meta[meta["median_node_visits"].isna()]