update: bootstrap track store

Author: stanfish06

.github/workflows/test-import.yml

@@ -27,7 +27,7 @@ jobs:
       run: |
         export PATH="$HOME/.cargo/bin:$PATH"
         make build
-        make sync-fresh
+        make fresh-sync
 
     - name: Test import
       run: |

Makefile

@@ -21,7 +21,7 @@ build: build-m4ri
 	CPLUS_INCLUDE_PATH=$(PROJECT_ROOT)/src/scloop/data:$(DM_PREFIX) \
 		uv build
 
-sync-fresh: build-m4ri
+fresh-sync: build-m4ri
 	CPLUS_INCLUDE_PATH=$(PROJECT_ROOT)/src/scloop/data:$(DM_PREFIX)/discrete-frechet-distance uv sync
 
 sync: clean build-m4ri

src/scloop/data/analysis_containers.py

@@ -1,17 +1,45 @@
 # Copyright 2025 Zhiyuan Yu (Heemskerk's lab, University of Michigan)
 from __future__ import annotations
 
+from typing import Optional
+
 import numpy as np
 from pydantic import ConfigDict, Field
 from pydantic.dataclasses import dataclass
 
 
+@dataclass(config=ConfigDict(arbitrary_types_allowed=True))
+class LoopMatch:
+    idx_bootstrap: int
+    birth_bootstrap: float
+    death_bootstrap: float
+    target_class_idx: int
+    geometric_distance: Optional[float] = None
+    neighbor_rank: Optional[int] = None
+    extra: dict = Field(default_factory=dict)
+
+
+@dataclass(config=ConfigDict(arbitrary_types_allowed=True))
+class LoopTrack:
+    source_class_idx: int
+    birth_root: float
+    death_root: float
+    matches: list[LoopMatch] = Field(default_factory=list)
+
+    def presence_prob(self, n_bootstraps: int) -> float:
+        if n_bootstraps == 0:
+            return 0.0
+        hit_boots = {m.idx_bootstrap for m in self.matches}
+        return len(hit_boots) / n_bootstraps
+
+
 @dataclass
 class BootstrapAnalysis:
     num_bootstraps: int = 0
     persistence_diagrams: list[list] = Field(default_factory=list)
     cocycles: list[list] = Field(default_factory=list)
     loop_representatives: list[list[list[list[int]]]] = Field(default_factory=list)
+    loop_tracks: dict[int, LoopTrack] = Field(default_factory=dict)
 
 
 @dataclass(config=ConfigDict(arbitrary_types_allowed=True))

src/scloop/data/containers.py

@@ -14,10 +14,15 @@
     compute_loop_homological_equivalence,
     compute_persistence_diagram_and_cocycles,
 )
-from .analysis_containers import BootstrapAnalysis, HodgeAnalysis
+from .analysis_containers import (
+    BootstrapAnalysis,
+    HodgeAnalysis,
+    LoopMatch,
+    LoopTrack,
+)
 from .loop_reconstruction import reconstruct_n_loop_representatives
 from .metadata import BootstrapMeta, ScloopMeta
-from .types import Diameter_t, Index_t, IndexListDownSample, Size_t
+from .types import Diameter_t, Index_t, IndexListDownSample, LoopDistMethod, Size_t
 from .utils import (
     decode_edges,
     decode_triangles,
@@ -239,7 +244,6 @@ def _compute_loop_representatives(
             return [], []
         top_k = min(top_k, loop_births.size)
 
-        # get top k homology classes
         indices_top_k = np.argpartition(loop_deaths - loop_births, -top_k)[-top_k:]
 
         dm_upper = triu(pairwise_distance_matrix, k=1).tocoo()
@@ -298,9 +302,11 @@ def _assess_bootstrap_geometric_equivalence(
         source_class_idx: int,
         target_class_idx: int,
         idx_bootstrap: int = 0,
+        method: LoopDistMethod = "hausdorff",
     ) -> tuple[int, int, float]:
         assert self.loop_representatives is not None
         assert self.bootstrap_data is not None
+        assert self.meta.preprocess is not None
         assert self.meta.preprocess.embedding_method is not None
 
         if idx_bootstrap >= len(self.bootstrap_data.loop_representatives):
@@ -334,7 +340,7 @@ def _assess_bootstrap_geometric_equivalence(
                     distances.append(np.nan)
 
         mean_distance = np.nanmean(distances) if distances else np.nan
-        return (source_class_idx, target_class_idx, mean_distance)
+        return (source_class_idx, target_class_idx, float(mean_distance))
 
     def _assess_bootstrap_homology_equivalence(
         self,
@@ -372,6 +378,26 @@ def _assess_bootstrap_homology_equivalence(
         )
         return (source_class_idx, target_class_idx, any(r == 0 for r in results))
 
+    def _ensure_loop_tracks(self) -> None:
+        if self.bootstrap_data is None:
+            return
+        if self.bootstrap_data.loop_tracks:
+            return
+        if self.persistence_diagram is None or self.loop_representatives is None:
+            return
+        loop_births = np.array(self.persistence_diagram[1][0], dtype=np.float32)
+        loop_deaths = np.array(self.persistence_diagram[1][1], dtype=np.float32)
+        top_k = len(self.loop_representatives)
+        if top_k == 0:
+            return
+        indices_top_k = np.argpartition(loop_deaths - loop_births, -top_k)[-top_k:]
+        for track_idx, loop_idx in enumerate(indices_top_k):
+            birth = float(loop_births[loop_idx])
+            death = float(loop_deaths[loop_idx])
+            self.bootstrap_data.loop_tracks[track_idx] = LoopTrack(
+                source_class_idx=track_idx, birth_root=birth, death_root=death
+            )
+
     def _bootstrap(
         self,
         adata: AnnData,
@@ -392,7 +418,7 @@ def _bootstrap(
         verbose: bool = True,
         **nei_kwargs,
     ) -> None:
-        self.bootstrap_data = BootstrapAnalysis(num_bootstraps=n_bootstrap)
+        self.bootstrap_data = BootstrapAnalysis()
         if self.meta.bootstrap is None:
             self.meta.bootstrap = BootstrapMeta(indices_resample=[])
         else:
@@ -446,6 +472,7 @@ def _bootstrap(
                     for j in range(n_bootstrap_loop_classes):
                         task = executor.submit(
                             self._assess_bootstrap_geometric_equivalence,
+                            adata,
                             i,
                             j,
                             idx_bootstrap,
@@ -459,7 +486,6 @@ def _bootstrap(
             neighbor_indices, neighbor_distances = nearest_neighbor_per_row(
                 pairwise_result_matrix, k_neighbors_check_equivalence
             )
-
             """
             ========= homological matching =========
             - gf2 regression
@@ -478,10 +504,32 @@ def _bootstrap(
                                 idx_bootstrap,
                                 n_pairs_check_equivalence,
                             )
-                            tasks[task] = (si, tj, neighbor_distances[si, k])
+                            tasks[task] = (si, tj, neighbor_distances[si, k], k)
 
                 for task in as_completed(tasks):
-                    si, tj, geom_dist = tasks[task]
-                    _, _, is_equivalent = task.result()
-                    if is_equivalent:
-                        pass
+                    si, tj, geo_dist, neighbor_rank = tasks[task]
+                    _, _, is_homologically_equivalent = task.result()
+                    if self.bootstrap_data is not None and is_homologically_equivalent:
+                        self._ensure_loop_tracks()
+                        track = self.bootstrap_data.loop_tracks.get(si)
+                        birth_boot = float(
+                            self.bootstrap_data.persistence_diagrams[idx_bootstrap][1][
+                                0
+                            ][tj]
+                        )
+                        death_boot = float(
+                            self.bootstrap_data.persistence_diagrams[idx_bootstrap][1][
+                                1
+                            ][tj]
+                        )
+                        track.matches.append(
+                            LoopMatch(
+                                idx_bootstrap=idx_bootstrap,
+                                birth_bootstrap=birth_boot,
+                                death_bootstrap=death_boot,
+                                target_class_idx=tj,
+                                geometric_distance=float(geo_dist),
+                                neighbor_rank=neighbor_rank,
+                            )
+                        )
+            self.bootstrap_data.num_bootstraps += 1

src/scloop/data/types.py

@@ -6,6 +6,7 @@
 FeatureSelectionMethod = Literal["hvg", "delve", "none"]
 EmbeddingMethod = Literal["pca", "diffmap", "scvi"]
 EmbeddingNeighbors = Literal["pca", "scvi"]
+LoopDistMethod = Literal["hausdorff", "frechet"]
 
 Index_t = Annotated[int, Field(ge=0)]
 Size_t = Annotated[int, Field(ge=0)]