update: loop construction and code reorganization

Author: stanfish06

src/scloop/computing/__init__.py

@@ -0,0 +1,8 @@
+"""Lightweight compute helpers used by the public API layer."""
+# Copyright 2025 Zhiyuan Yu (Heemskerk's lab, University of Michigan)
+
+from .homology import (
+    compute_boundary_matrix_data,
+    compute_persistence_diagram,
+    compute_sparse_pairwise_distance,
+)

src/scloop/computing/homology.py

@@ -0,0 +1,67 @@
+# Copyright 2025 Zhiyuan Yu (Heemskerk's lab, University of Michigan)
+from __future__ import annotations
+
+import numpy as np
+from anndata import AnnData
+from scipy.sparse import csr_matrix
+from sklearn.neighbors import radius_neighbors_graph
+
+from ..data.metadata import ScloopMeta
+from ..data.ripser_lib import get_boundary_matrix, ripser
+from ..data.types import Diameter_t, IndexListDistMatrix
+from ..data.utils import encode_triangles_and_edges
+
+
+def compute_sparse_pairwise_distance(
+    adata: AnnData,
+    meta: ScloopMeta,
+    bootstrap: bool = False,
+    thresh: Diameter_t | None = None,
+    **nei_kwargs,
+) -> tuple[csr_matrix, IndexListDistMatrix | None]:
+    assert meta.preprocess is not None
+    assert meta.preprocess.embedding_method is not None
+    boot_idx = None
+    if bootstrap:
+        boot_idx = np.random.choice(
+            adata.shape[0], size=adata.shape[0], replace=True
+        ).tolist()
+    return (
+        radius_neighbors_graph(
+            X=adata.obsm[f"X_{meta.preprocess.embedding_method}"],
+            radius=thresh,
+            **nei_kwargs,
+        ),
+        boot_idx,
+    )
+
+
+def compute_persistence_diagram(
+    adata: AnnData, meta: ScloopMeta, thresh: Diameter_t | None = None, **nei_kwargs
+) -> tuple[list[np.ndarray], IndexListDistMatrix | None, csr_matrix]:
+    sparse_pairwise_distance_matrix, boot_idx = compute_sparse_pairwise_distance(
+        adata=adata, meta=meta, bootstrap=False, thresh=thresh, **nei_kwargs
+    )
+    result = ripser(
+        distance_matrix=sparse_pairwise_distance_matrix,
+        modulus=2,
+        dim_max=1,
+        threshold=thresh,
+        do_cocyles=True,
+    )
+    return result.births_and_deaths_by_dim, boot_idx, sparse_pairwise_distance_matrix
+
+
+def compute_boundary_matrix_data(
+    adata: AnnData, meta: ScloopMeta, thresh: Diameter_t | None = None, **nei_kwargs
+):
+    assert meta.preprocess is not None
+    assert meta.preprocess.num_vertices is not None
+    sparse_pairwise_distance_matrix, boot_idx = compute_sparse_pairwise_distance(
+        adata=adata, meta=meta, bootstrap=False, thresh=thresh, **nei_kwargs
+    )
+    result = get_boundary_matrix(sparse_pairwise_distance_matrix, thresh)
+    edge_ids, trig_ids = encode_triangles_and_edges(
+        np.array(result.triangle_vertices), meta.preprocess.num_vertices
+    )
+    return result, edge_ids, trig_ids, boot_idx, sparse_pairwise_distance_matrix

src/scloop/data/containers.py

@@ -1,19 +1,20 @@
 # Copyright 2025 Zhiyuan Yu (Heemskerk's lab, University of Michigan)
-from pydantic.dataclasses import dataclass
 from abc import abstractmethod
+from pydantic.dataclasses import dataclass
 from pydantic import ConfigDict, BaseModel, Field, field_validator, ValidationInfo
-from scipy.spatial import distance_matrix
 from .metadata import ScloopMeta
 from .analysis_containers import BootstrapAnalysis, HodgeAnalysis
 from anndata import AnnData
-from sklearn.neighbors import radius_neighbors_graph
 from scipy.sparse import csr_matrix
-from .ripser_lib import ripser, get_boundary_matrix
 from .loop_reconstruction import reconstruct_n_loop_representatives
 from .types import IndexListDistMatrix, Diameter_t, Size_t, Index_t
-from .utils import decode_edges, decode_triangles, encode_triangles_and_edges
-from numba import jit
+from .utils import decode_edges, decode_triangles
 import numpy as np
+from ..computing.homology import (
+    compute_sparse_pairwise_distance,
+    compute_persistence_diagram,
+    compute_boundary_matrix_data,
+)
 
 
 class BoundaryMatrix(BaseModel):
@@ -88,53 +89,35 @@ def _compute_sparse_pairwise_distance(
         thresh: Diameter_t | None = None,
         **nei_kwargs,
     ) -> tuple[csr_matrix, IndexListDistMatrix | None]:
-        assert self.meta.preprocess is not None
-        assert self.meta.preprocess.embedding_method is not None
-        boot_idx = None
-        if bootstrap:
-            boot_idx = np.random.choice(
-                adata.shape[0], size=adata.shape[0], replace=True
-            ).tolist()
-        return radius_neighbors_graph(
-            X=adata.obsm[f"X_{self.meta.preprocess.embedding_method}"],
-            radius=thresh,
+        return compute_sparse_pairwise_distance(
+            adata=adata,
+            meta=self.meta,
+            bootstrap=bootstrap,
+            thresh=thresh,
             **nei_kwargs,
-        ), boot_idx
+        )
 
     def _compute_homology(
         self, adata: AnnData, thresh: Diameter_t | None = None, **nei_kwargs
     ) -> None:
-        sparse_pairwise_distance_matrix, _ = self._compute_sparse_pairwise_distance(
-            adata=adata, bootstrap=False, thresh=thresh, **nei_kwargs
+        persistence_diagram, _, _ = compute_persistence_diagram(
+            adata=adata, meta=self.meta, thresh=thresh, **nei_kwargs
         )
-        result = ripser(
-            distance_matrix=sparse_pairwise_distance_matrix,
-            modulus=2,
-            dim_max=1,
-            threshold=thresh,
-            do_cocyles=True,
-        )
-        self.persistence_diagram = result.births_and_deaths_by_dim
+        self.persistence_diagram = persistence_diagram
 
     def _compute_boundary_matrix(
         self, adata: AnnData, thresh: Diameter_t | None = None, **nei_kwargs
     ) -> None:
-        assert self.meta.preprocess
-        assert self.meta.preprocess.num_vertices
-        sparse_pairwise_distance_matrix, _ = self._compute_sparse_pairwise_distance(
-            adata=adata, bootstrap=False, thresh=thresh, **nei_kwargs
-        )
-        result = get_boundary_matrix(sparse_pairwise_distance_matrix, thresh)
-        edge_ids, trig_ids = encode_triangles_and_edges(
-            np.array(result.triangle_vertices), self.meta.preprocess.num_vertices
+        result, *_ = compute_boundary_matrix_data(
+            adata=adata, meta=self.meta, thresh=thresh, **nei_kwargs
         )
         self.boundary_matrix_d1 = BoundaryMatrixD1(
             num_vertices=self.meta.preprocess.num_vertices,
             data=([], []),
             shape=(0, 0),
             row_simplex_ids=[],
             col_simplex_ids=[],
-            col_simplex_diams=result.traingle_diameters,
+            col_simplex_diams=result.triangle_diameters,
         )
 
     def _compute_loop_representatives(self):

src/scloop/data/loop_reconstruction.py

@@ -1,7 +1,11 @@
 # Copyright 2025 Zhiyuan Yu (Heemskerk's lab, University of Michigan)
-import numpy as np
+from __future__ import annotations
+
+import math
+from typing import Iterable, List, Sequence, Tuple
+
 import igraph as ig
-from typing import List, Tuple
+import numpy as np
 
 
 def reconstruct_n_loop_representatives(
@@ -18,4 +22,129 @@ def reconstruct_n_loop_representatives(
     loop_upper_pct: float = 95,
     n_max_cocycles: int = 10,
 ) -> Tuple[List[List[int]], List[float]]:
-    return ([], [])
+    """
+    Reconstruct diverse loop representatives using Yen-style deviation rounds.
+    """
+    if n <= 0 or len(edges) == 0:
+        return [], []
+
+    filt_t = loop_birth + (loop_death - loop_birth) * life_pct
+
+    # Parse cocycle edges (each entry is [[i, j], coeff])
+    cocycle_edges: list[tuple[int, int]] = []
+    for simplex in cocycles_dim1:
+        try:
+            verts, coeff = simplex
+        except ValueError:
+            continue
+        if coeff == 0 or len(verts) != 2:
+            continue
+        cocycle_edges.append((int(verts[0]), int(verts[1])))
+        if len(cocycle_edges) == n_max_cocycles:
+            break
+
+    edge_births = np.asarray(edge_births)
+    mask = edge_births <= filt_t
+    if not np.any(mask):
+        return [], []
+    edges_filt = [e for e, keep in zip(edges, mask) if keep]
+    weights_filt = edge_births[mask].tolist()
+
+    sources = [e[0] for e in edges_filt] + [e[0] for e in cocycle_edges]
+    destinations = [e[1] for e in edges_filt] + [e[1] for e in cocycle_edges]
+    weights = weights_filt + [math.inf] * len(cocycle_edges)
+    if len(sources) == 0:
+        return [], []
+
+    n_vertices = max(max(sources), max(destinations)) + 1
+    g = ig.Graph(n=n_vertices, edges=list(zip(sources, destinations)), directed=False)
+    g.es["weight"] = weights
+
+    cycles_pool: list[list[int]] = []
+    cycles_dist: list[float] = []
+
+    for _ in range(n_force_deviate):
+        paths_this_round: list[list[int]] = []
+        for i, j in cocycle_edges:
+            paths = _k_shortest_paths(g, i, j, n_reps_per_loop)
+            if not paths:
+                continue
+            for path in paths:
+                dist = _path_weight(g, path)
+                cycles_pool.append(path)
+                paths_this_round.append(path)
+                cycles_dist.append(dist)
+
+        # Force deviation
+        for path in paths_this_round:
+            for u, v in zip(path[:-1], path[1:]):
+                try:
+                    eid = g.get_eid(u, v, directed=False)
+                except ig._igraph.InternalError:
+                    continue
+                g.es[eid]["weight"] = math.inf
+
+    return _select_diverse_loops(
+        cycles=cycles_pool,
+        distances=cycles_dist,
+        n=n,
+        lower_pct=loop_lower_pct,
+        upper_pct=loop_upper_pct,
+    )
+
+
+def _k_shortest_paths(g: ig.Graph, source: int, target: int, k: int) -> list[list[int]]:
+    if source == target:
+        return []
+    try:
+        return g.get_k_shortest_paths(
+            source, target, k=k, weights=g.es["weight"], mode="ALL"
+        )
+    except ig._igraph.InternalError:
+        return []
+
+
+def _path_weight(g: ig.Graph, path: Sequence[int]) -> float:
+    if len(path) < 2:
+        return math.inf
+    weight = 0.0
+    for u, v in zip(path[:-1], path[1:]):
+        try:
+            eid = g.get_eid(u, v, directed=False)
+        except ig._igraph.InternalError:
+            return math.inf
+        w = g.es[eid]["weight"]
+        weight += float(w)
+    return weight
+
+
+def _select_diverse_loops(
+    cycles: Iterable[Sequence[int]],
+    distances: Iterable[float],
+    n: int,
+    lower_pct: float,
+    upper_pct: float,
+) -> Tuple[List[List[int]], List[float]]:
+    pairs = sorted(
+        [(float(d), list(c)) for d, c in zip(distances, cycles) if math.isfinite(d)],
+        key=lambda x: x[0],
+    )
+    if not pairs:
+        return [], []
+
+    n_total = len(pairs)
+    n_return = min(n_total, n)
+    if n_return == 1:
+        idxs = [n_total // 2]
+    else:
+        step = (upper_pct - lower_pct) / (n_return - 1)
+        idxs = []
+        for i in range(n_return):
+            pct = (lower_pct + step * i) / 100
+            idx = min(int(math.floor(n_total * pct)), n_total - 1)
+            idxs.append(idx)
+
+    selected = [pairs[i] for i in idxs]
+    dists = [p[0] for p in selected]
+    loops = [p[1] for p in selected]
+    return loops, dists