update: auto-choose threshold. Fixes #70

Author: stanfish06

src/scloop/computing/homology.py

@@ -435,9 +435,7 @@ def compute_loop_homological_equivalence(
         # replace last columns (already sorted by diameters) with identity columns
         n_extra_edges = min(len(n_hubs_edges), max_n_edges_relaxation)
         if n_extra_edges > 0:
-            one_ridx_A[: 3 * n_extra_edges :] = np.repeat(
-                n_hubs_edges[:n_extra_edges], 3
-            )
+            one_ridx_A[: 3 * n_extra_edges] = np.repeat(n_hubs_edges[:n_extra_edges], 3)
             results_relax, solutions_relax = solve_multiple_gf2_m4ri(
                 one_ridx_A=one_ridx_A.tolist(),
                 one_cidx_A=one_cidx_A.tolist(),

src/scloop/data/constants.py

@@ -29,6 +29,8 @@
 DEFAULT_MAX_COLUMNS_BOUNDARY_MATRIX: int = 10000
 DEFAULT_N_BOOTSTRAP: int = 10
 
+DEFAULT_AUTO_THRESHOLD_FACTOR: float = 1.25
+
 DEFAULT_NOISE_SCALE: float = 1e-6
 DEFAULT_N_REPS_PER_LOOP: int = 4
 DEFAULT_N_COCYCLES_USED: int = 3
@@ -39,7 +41,7 @@
 DEFAULT_K_NEIGHBORS_CHECK_EQUIVALENCE: int = 1
 DEFAULT_EXTRA_DIAM_EQUIVALENCE: float = 1.0
 DEFAULT_WITH_RELAXATION_EQUIVALENCE: bool = True
-DEFAULT_N_HUBS_RELAXATION_EQUIVALENCE: int = 2
+DEFAULT_N_HUBS_RELAXATION_EQUIVALENCE: int = 3
 DEFAULT_MAX_N_EDGES_RELAXATION_EQUIVALENCE: int = 500
 
 DEFAULT_WEIGHT_HODGE: float = 0.5

src/scloop/tools/_loops.py

@@ -17,8 +17,10 @@
     TimeElapsedColumn,
     TimeRemainingColumn,
 )
+from scipy.spatial.distance import pdist
 
 from ..data.constants import (
+    DEFAULT_AUTO_THRESHOLD_FACTOR,
     DEFAULT_K_NEIGHBORS_CHECK_EQUIVALENCE,
     DEFAULT_MAX_COLUMNS_BOUNDARY_MATRIX,
     DEFAULT_MAXITER_EIGENDECOMPOSITION,
@@ -103,6 +105,39 @@ def find_loops(
             meta.bootstrap = BootstrapMeta()
         meta.bootstrap.life_pct = tightness_loops
         hd: HomologyData = HomologyData(meta=meta)
+        """
+        ============ Auto-choose PH threshold ============
+        - upper bound would be the max pw dist
+        - choose a value such that all 1-loop dies
+        ==================================================
+        """
+        if threshold_homology is None:
+            assert meta.preprocess is not None
+            assert meta.preprocess.embedding_method is not None
+            emb = adata.obsm[f"X_{meta.preprocess.embedding_method}"]
+            selected_indices = (
+                meta.preprocess.indices_downsample
+                if meta.preprocess.indices_downsample is not None
+                else list(range(emb.shape[0]))
+            )
+            max_pw_dist = float(pdist(emb[selected_indices]).max())
+            if verbose:
+                logger.info(
+                    f"Auto-threshold: max pairwise distance = {max_pw_dist:.4f}"
+                )
+            hd._compute_homology(adata=adata, thresh=max_pw_dist)
+            # need some room for bootstrap, loops will die later for fewer points
+            auto_factor = (kwargs_bootstrap or {}).get(
+                "auto_threshold_factor", DEFAULT_AUTO_THRESHOLD_FACTOR
+            )
+            max_h1_death = float(np.max(hd.persistence_diagram[1][1]))
+            threshold_homology = max_h1_death * auto_factor
+            if verbose:
+                logger.info(
+                    f"Auto-threshold: max H1 death = {max_h1_death:.4f}, "
+                    f"using threshold = {threshold_homology:.4f} (factor={auto_factor})"
+                )
+
         sparse_dist_mat = hd._compute_homology(adata=adata, thresh=threshold_homology)
         boundary_thresh = threshold_boundary
         if boundary_thresh is None: