From 0ac876a6fc999a80338c1b2baf66297d5de36a90 Mon Sep 17 00:00:00 2001 From: Stan Date: Tue, 23 Jun 2026 15:23:32 -0400 Subject: [PATCH 1/4] hyperplane update --- src/scloop/plotting/_homology.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/scloop/plotting/_homology.py b/src/scloop/plotting/_homology.py index 7f4a345..b3f445e 100644 --- a/src/scloop/plotting/_homology.py +++ b/src/scloop/plotting/_homology.py @@ -59,6 +59,10 @@ def _get_track_loop( return tracked_pairs +def _compute_loop_hyperplane(data: HomologyData, track_id: int): + pass + + @validate_call(config=ConfigDict(arbitrary_types_allowed=True)) def hist_lifetimes( adata: AnnData, From ad5e0b4040e449400ec4396592a8e4a000183df3 Mon Sep 17 00:00:00 2001 From: Stan Date: Wed, 24 Jun 2026 13:50:10 -0400 Subject: [PATCH 2/4] update --- src/scloop/plotting/_homology.py | 34 ++++++++++++++++++++++++++++++-- 1 file changed, 32 insertions(+), 2 deletions(-) diff --git a/src/scloop/plotting/_homology.py b/src/scloop/plotting/_homology.py index b3f445e..1da725f 100644 --- a/src/scloop/plotting/_homology.py +++ b/src/scloop/plotting/_homology.py @@ -8,6 +8,7 @@ from anndata import AnnData from matplotlib.axes import Axes from pydantic import ConfigDict, validate_call +from scipy.linalg import svd from ..data.analysis_containers import BootstrapAnalysis from ..data.constants import DEFAULT_DPI, DEFAULT_FIGSIZE, SCLOOP_UNS_KEY @@ -59,8 +60,37 @@ def _get_track_loop( return tracked_pairs -def _compute_loop_hyperplane(data: HomologyData, track_id: int): - pass +def _compute_loop_decomposition( + adata: AnnData, + basis: str, + track_ids: int | list[int], + key_homology: str = SCLOOP_UNS_KEY, +): + hdata: HomologyData = adata.uns[key_homology] + emb = adata.obsm[basis] + assert emb is np.ndarray + Us = [] + match track_ids: + case int(): + X = np.concatenate( + hdata._get_loop_embedding( + selector=track_ids, embedding_alt=emb, include_bootstrap=True + ), + axis=0, + ) + U, _, _ = svd(X.T, full_matrices=False) + Us.append(U) + case list(): + for tid in track_ids: + X = np.concatenate( + hdata._get_loop_embedding( + selector=tid, embedding_alt=emb, include_bootstrap=True + ), + axis=0, + ) + U, _, _ = svd(X.T, full_matrices=False) + Us.append(U) + return Us @validate_call(config=ConfigDict(arbitrary_types_allowed=True)) From 522a6604c2db77b80a280e2d0ce5a5fc13b4b9a9 Mon Sep 17 00:00:00 2001 From: Stan Date: Fri, 26 Jun 2026 17:15:21 -0400 Subject: [PATCH 3/4] fix: palette issue --- src/scloop/plotting/_homology.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/scloop/plotting/_homology.py b/src/scloop/plotting/_homology.py index 1da725f..cad3eb6 100644 --- a/src/scloop/plotting/_homology.py +++ b/src/scloop/plotting/_homology.py @@ -23,6 +23,8 @@ "loops", ] +DEFAULT_GLASBEY_BLOCK_SIZE = 5 + # ugly function, fix that def _get_track_loop( @@ -242,7 +244,7 @@ def bar_lifetimes( track_ids = track_ids or [] n_tracks = len(track_ids) if n_tracks > 0: - block_size = 5 + block_size = DEFAULT_GLASBEY_BLOCK_SIZE cmap = glasbey.create_block_palette(block_sizes=[block_size] * n_tracks) cmap = [cmap[i : i + block_size] for i in range(0, len(cmap), block_size)] for i, src_tid in enumerate(track_ids): @@ -355,7 +357,7 @@ def persistence_diagram( track_ids = track_ids or [] n_tracks = len(track_ids) if n_tracks > 0: - block_size = 5 + block_size = DEFAULT_GLASBEY_BLOCK_SIZE cmap = glasbey.create_block_palette(block_sizes=[block_size] * n_tracks) cmap = [cmap[i : i + block_size] for i in range(0, len(cmap), block_size)] for i, src_tid in enumerate(track_ids): @@ -465,7 +467,7 @@ def loops( n_selectors = len(selectors) if n_selectors > 0: - block_size = 5 + block_size = DEFAULT_GLASBEY_BLOCK_SIZE cmap = glasbey.create_block_palette(block_sizes=[block_size] * n_selectors) cmap = [cmap[i : i + block_size] for i in range(0, len(cmap), block_size)] else: @@ -556,7 +558,7 @@ def _loops_for_selector( ax.plot( loop[:, components[0]], loop[:, components[1]], - color=cmap[i][j % block_size], + color=cmap[i][(block_size - 1) - j % block_size], **(kwargs_scatter or {}), ) From e1adebb56e9f00b913e09346f093f9769273b381 Mon Sep 17 00:00:00 2001 From: Stan Date: Fri, 26 Jun 2026 19:44:36 -0400 Subject: [PATCH 4/4] update: consensus diffusion plane --- src/scloop/plotting/__init__.py | 2 + src/scloop/plotting/_homology.py | 66 +++++++++++++++++++------------- 2 files changed, 42 insertions(+), 26 deletions(-) diff --git a/src/scloop/plotting/__init__.py b/src/scloop/plotting/__init__.py index e4c7ed9..6f7c565 100644 --- a/src/scloop/plotting/__init__.py +++ b/src/scloop/plotting/__init__.py @@ -4,6 +4,7 @@ from ._homology import ( bar_lifetimes, hist_lifetimes, + loop_embedding, loops, persistence_diagram, ) @@ -12,6 +13,7 @@ __all__ = [ "bar_lifetimes", "hist_lifetimes", + "loop_embedding", "loops", "loop_edge_embedding", "loop_edge_overlay", diff --git a/src/scloop/plotting/_homology.py b/src/scloop/plotting/_homology.py index cad3eb6..0ab35b0 100644 --- a/src/scloop/plotting/_homology.py +++ b/src/scloop/plotting/_homology.py @@ -20,6 +20,7 @@ "hist_lifetimes", "bar_lifetimes", "persistence_diagram", + "loop_embedding", "loops", ] @@ -62,37 +63,50 @@ def _get_track_loop( return tracked_pairs -def _compute_loop_decomposition( +def loop_embedding( adata: AnnData, basis: str, track_ids: int | list[int], + ndims: int = 2, + key_added: str = "loops", key_homology: str = SCLOOP_UNS_KEY, -): +) -> None: hdata: HomologyData = adata.uns[key_homology] - emb = adata.obsm[basis] - assert emb is np.ndarray - Us = [] - match track_ids: - case int(): - X = np.concatenate( - hdata._get_loop_embedding( - selector=track_ids, embedding_alt=emb, include_bootstrap=True - ), - axis=0, - ) - U, _, _ = svd(X.T, full_matrices=False) - Us.append(U) - case list(): - for tid in track_ids: - X = np.concatenate( - hdata._get_loop_embedding( - selector=tid, embedding_alt=emb, include_bootstrap=True - ), - axis=0, - ) - U, _, _ = svd(X.T, full_matrices=False) - Us.append(U) - return Us + emb = np.asarray(adata.obsm[basis]) + track_list = [track_ids] if isinstance(track_ids, int) else list(track_ids) + + planes = [] + for tid in track_list: + X = np.concatenate( + hdata._get_loop_embedding( + selector=tid, embedding_alt=emb, include_bootstrap=True + ), + axis=0, + ) + X = X - X.mean(axis=0, keepdims=True) + U, _, _ = svd(X.T, full_matrices=False) + planes.append(U[:, :ndims]) + + if len(planes) == 1: + plane = planes[0] + score = 1.0 + else: + """ + Steps to figure out a consensus plane: + 1. Given per-loop diffusion planes Us, the projection matricies (to the target diffusion plane) are U^t. We want to figure out a consensus plane P. + 2. The goal is max sum_k(tr(P^t Uk Uk^t P)) such that the projection of P to all diffusion planes are maximized. + - max sum_k(tr(P^t Uk Uk^t P)) + - max tr(P^t sum_k(Uk Uk^t) P) + 3. Let P = [ v1 v2 ], where v1 and v2 are top eigenvecs of sum_k(Uk Uk^t) + 4. tr(P^t sum_k(Uk Uk^t) P) = tr([ [ l1 0 ] [ 0 l2 ] ]) = l1 + l2 + """ + M = sum(U @ U.T for U in planes) + evals, evecs = np.linalg.eigh(M) # ascending + plane = evecs[:, ::-1][:, :ndims] # top ndims, descending + score = float(evals[::-1][:ndims].sum() / (ndims * len(planes))) + + adata.obsm[f"X_{key_added}"] = (emb - emb.mean(axis=0, keepdims=True)) @ plane + adata.uns[f"{key_added}_loop_embedding"] = {"plane": plane, "score": score} @validate_call(config=ConfigDict(arbitrary_types_allowed=True))