@@ -201,6 +201,7 @@ def run_single_bootstrap(
     extra_diameter_homology_equivalence: float = DEFAULT_EXTRA_DIAM_EQUIVALENCE,
     filter_column_homology_equivalence: bool = True,
     column_trim_method: ColumnTrimMethod = DEFAULT_COLUMN_TRIM_METHOD,
+    death_scale: float | None = None,
     n_neighbors_column_trim: int = DEFAULT_N_NEIGHBORS_COLUMN_TRIM,
     full_pairwise_distance_matrix: csr_matrix | None = None,
     full_vertex_ids: list[int] | None = None,
@@ -453,6 +454,8 @@ def run_single_bootstrap(
                 max_column_diameter=max_column_diameter,
                 cocycle_edge_mask=cocycle_edge_masks[source_idx],
                 column_trim_method=column_trim_method,
+                embedding=embedding,
+                death_scale=death_scale,
                 column_scores=(
                     source_loop.column_proximity_scores(
                         original_boundary_matrix_d1,
@@ -496,6 +499,10 @@ def run_bootstrap_pipeline(
     **kwargs,
 ) -> list[BootstrapResult]:
     results: list[BootstrapResult] = []
+    _original_deaths = [lc.death for lc in original_loop_classes if lc is not None]
+    global_death_scale = (
+        max(_original_deaths) if _original_deaths else meta.bootstrap.threshold_homology
+    )
 
     full_pairwise_distance_matrix: csr_matrix | None = None
     full_vertex_ids: list[int] | None = None
@@ -512,9 +519,7 @@ def run_bootstrap_pipeline(
     if kwargs.get("column_trim_method", DEFAULT_COLUMN_TRIM_METHOD) == "loop_proximity":
         assert meta.preprocess is not None
         assert meta.preprocess.embedding_method is not None
-        warm_embedding = np.array(
-            adata.obsm[f"X_{meta.preprocess.embedding_method}"]
-        )
+        warm_embedding = np.array(adata.obsm[f"X_{meta.preprocess.embedding_method}"])
         for loop_class in original_loop_classes:
             if loop_class is not None and loop_class.representatives is not None:
                 loop_class.column_proximity_scores(
@@ -537,6 +542,7 @@ def run_bootstrap_pipeline(
                 meta=meta,
                 original_loop_classes=original_loop_classes,
                 original_boundary_matrix_d1=original_boundary_matrix_d1,
+                death_scale=global_death_scale,
                 full_pairwise_distance_matrix=full_pairwise_distance_matrix,
                 full_vertex_ids=full_vertex_ids,
                 reconstruct_on_full_data=reconstruct_on_full_data,
 
@@ -3,8 +3,19 @@
 
 from typing import Mapping, Sequence
 
+import numpy as np
+
 from ..data.types import HomotopyCoherenceMethod
 from ..data.utils import edge_idx_decode
+from ..utils.distance_metrics import compute_loop_frechet
+
+
+def global_h1_death_scale(persistence_diagram: list | None) -> float | None:
+    if not persistence_diagram or len(persistence_diagram) < 2:
+        return None
+    deaths = np.asarray(persistence_diagram[1][1], dtype=float)
+    finite = deaths[np.isfinite(deaths)]
+    return float(finite.max()) if finite.size else None
 
 
 def _flip_triangle(
@@ -33,25 +44,65 @@ def _flip_triangle(
     return cycle ^ triangle
 
 
-def _perimeter(cycle: frozenset[int], edge_lengths: Mapping[int, float]) -> float:
-    return sum(edge_lengths[edge] for edge in cycle)
+def _cycle_vertices(cycle: frozenset[int], num_vertices: int) -> list[int]:
+    adjacency: dict[int, list[int]] = {}
+    for edge in cycle:
+        tail, head = edge_idx_decode(edge, num_vertices)
+        adjacency.setdefault(tail, []).append(head)
+        adjacency.setdefault(head, []).append(tail)
+    start = next(iter(adjacency))
+    order = [start]
+    previous, current = -1, start
+    while len(order) <= len(adjacency):
+        neighbours = adjacency[current]
+        if len(neighbours) != 2:
+            break
+        first, second = neighbours
+        following = first if first != previous else second
+        if following == start:
+            break
+        order.append(following)
+        previous, current = current, following
+    return order
+
+
+def _cycle_coords(
+    cycle: frozenset[int], num_vertices: int, embedding: np.ndarray
+) -> np.ndarray:
+    vertices = _cycle_vertices(cycle, num_vertices)
+    return np.ascontiguousarray(embedding[vertices], dtype=np.float64)
+
+
+def _frechet_to_target(
+    cycle: frozenset[int],
+    target_coords: np.ndarray,
+    num_vertices: int,
+    embedding: np.ndarray,
+) -> float:
+    cycle_coords = _cycle_coords(cycle, num_vertices, embedding)
+    return float(
+        min(
+            compute_loop_frechet(cycle_coords, target_coords),
+            compute_loop_frechet(target_coords, cycle_coords),
+        )
+    )
 
 
 def _greedy_reversal(
     source: frozenset[int],
     target: frozenset[int],
     triangles: tuple[frozenset[int], ...],
-    edge_lengths: Mapping[int, float],
     num_vertices: int,
+    embedding: np.ndarray,
+    death_scale: float | None = None,
 ) -> float | None:
-    """Follow one valid triangle order greedily by target perimeter distance."""
-
-    target_length = _perimeter(target, edge_lengths)
-    source_length = _perimeter(source, edge_lengths)
+    target_coords = _cycle_coords(target, num_vertices, embedding)
+    source_coords = _cycle_coords(source, num_vertices, embedding)
     cycle = source
-    length = source_length
     reversal = 0.0
     total_step_abs = 0.0
+    max_step_abs = 0.0
+    n_steps = 0
     remaining = set(range(len(triangles)))
 
     while remaining:
@@ -63,47 +114,59 @@ def _greedy_reversal(
             if next_cycle is None:
                 continue
 
-            current_remaining_cost = abs(target_length - length)
-            next_length = _perimeter(next_cycle, edge_lengths)
-            next_remaining_cost = abs(target_length - next_length)
-            step_reversal = (
-                max(
-                    next_remaining_cost - current_remaining_cost,
-                    0.0,
-                )
-                / length
+            current_remaining_cost = _frechet_to_target(
+                cycle, target_coords, num_vertices, embedding
+            )
+            next_remaining_cost = _frechet_to_target(
+                next_cycle, target_coords, num_vertices, embedding
+            )
+            current_departed_cost = _frechet_to_target(
+                cycle, source_coords, num_vertices, embedding
+            )
+            next_departed_cost = _frechet_to_target(
+                next_cycle, source_coords, num_vertices, embedding
+            )
+            step_reversal = max(
+                next_remaining_cost - current_remaining_cost,
+                0.0,
+            )
+            step_abs = abs(next_remaining_cost - current_remaining_cost) + abs(
+                next_departed_cost - current_departed_cost
             )
-            step_abs = abs(next_length - length) / length
             key = (next_remaining_cost, triangle_index)
             if best_key is None or key < best_key:
                 best_key = key
                 best_step = (
                     triangle_index,
                     next_cycle,
-                    next_length,
                     step_reversal,
                     step_abs,
                 )
 
         if best_step is None:
             return None
 
-        triangle_index, cycle, length, step_reversal, step_abs = best_step
+        triangle_index, cycle, step_reversal, step_abs = best_step
         remaining.remove(triangle_index)
         reversal += step_reversal
         total_step_abs += step_abs
+        max_step_abs = max(max_step_abs, step_abs)
+        n_steps += 1
 
-    return 1 - reversal / total_step_abs if cycle == target else None
+    if cycle != target:
+        return None
+    # multiply by 2 due to two way frechet
+    return 1 - min(max_step_abs / (2 * death_scale), 1)
 
 
 def path_finding_coherence(
     source_edges: Sequence[int],
     target_edges: Sequence[int],
     triangles: Sequence[Sequence[int]],
-    edge_lengths: Mapping[int, float],
     num_vertices: int,
+    embedding: np.ndarray,
+    death_scale: float | None = None,
 ) -> float | None:
-    """Return coherence from a greedy valid triangle order."""
 
     source = frozenset(source_edges)
     target = frozenset(target_edges)
@@ -112,35 +175,40 @@ def path_finding_coherence(
         source,
         target,
         triangle_sets,
-        edge_lengths,
         num_vertices,
+        embedding,
+        death_scale,
     )
+
     reverse = _greedy_reversal(
         target,
         source,
         triangle_sets,
-        edge_lengths,
         num_vertices,
+        embedding,
+        death_scale,
     )
 
     candidates = [cost for cost in (forward, reverse) if cost is not None]
     if not candidates:
         return None
-    return max(candidates)
+    return min(candidates)
 
 
 def compute_coherence(
     source_edges: Sequence[int],
     target_edges: Sequence[int],
     triangles: Sequence[Sequence[int]],
-    edge_lengths: Mapping[int, float],
     num_vertices: int,
+    embedding: np.ndarray,
+    death_scale: float | None = None,
     method: HomotopyCoherenceMethod = "path_finding",
 ) -> float | None:
     return path_finding_coherence(
         source_edges,
         target_edges,
         triangles,
-        edge_lengths,
         num_vertices,
+        embedding,
+        death_scale,
     )
@@ -157,6 +157,7 @@ def check_homological_equivalence(
     column_trim_method: ColumnTrimMethod = DEFAULT_COLUMN_TRIM_METHOD,
     column_scores: np.ndarray | None = None,
     embedding: np.ndarray | None = None,
+    death_scale: float | None = None,
     n_neighbors_column_trim: int = DEFAULT_N_NEIGHBORS_COLUMN_TRIM,
 ) -> LoopClassEquivalence:
     if len(source_loops) == 0 or len(target_loops) == 0:
@@ -208,7 +209,7 @@ def check_homological_equivalence(
         column_scores=column_scores,
     )
 
-    if not compute_homotopy_coherence:
+    if not compute_homotopy_coherence or embedding is None:
         return result
 
     assert isinstance(loop_edges_a, LoopEdges)
@@ -243,6 +244,8 @@ def check_homological_equivalence(
                 ],
                 edge_lengths=edge_lengths,
                 num_vertices=boundary_matrix_d1.num_vertices,
+                embedding=embedding,
+                death_scale=death_scale,
                 method=homotopy_coherence_method,
             )
         )
@@ -262,16 +265,18 @@ def check_homological_equivalence(
                 source_edges=tuple(source ^ relaxation_edges),
                 target_edges=tuple(target),
                 triangles=triangles,
-                edge_lengths=edge_lengths,
                 num_vertices=boundary_matrix_d1.num_vertices,
+                embedding=embedding,
+                death_scale=death_scale,
                 method=homotopy_coherence_method,
             ),
             compute_coherence(
                 source_edges=tuple(source),
                 target_edges=tuple(target ^ relaxation_edges),
                 triangles=triangles,
-                edge_lengths=edge_lengths,
                 num_vertices=boundary_matrix_d1.num_vertices,
+                embedding=embedding,
+                death_scale=death_scale,
                 method=homotopy_coherence_method,
             ),
         ]
 
@@ -27,6 +27,7 @@
     compute_persistence_diagram_and_cocycles,
 )
 from ..computing.loops import compute_loop_representatives
+from ..computing.coherence import global_h1_death_scale
 from ..computing.matching import (
     check_homological_equivalence,
     cocycle_to_edge_mask,
@@ -815,6 +816,8 @@ def _assess_bootstrap_homology_equivalence(
             max_column_diameter=max_column_diameter,
             cocycle_edge_mask=cocycle_edge_mask,
             column_trim_method=column_trim_method,
+            embedding=embedding,
+            death_scale=global_h1_death_scale(self.persistence_diagram),
             column_scores=(
                 source_loop_class.column_proximity_scores(
                     self.boundary_matrix_d1,
 
@@ -14,6 +14,7 @@
     FeatureSelectionMethod,
     IndexListDownSample,
     Percent_t,
+    PositiveFloat,
     Size_t,
 )
 
@@ -172,6 +173,8 @@ def from_hdf5_group(cls, group: h5py.Group) -> PreprocessMeta:
 class BootstrapMeta(BaseModel):
     indices_resample: list[IndexListDownSample] | None = None
     life_pct: Percent_t | None = None
+    threshold_homology: PositiveFloat | None = None
+    threshold_boundary: PositiveFloat | None = None
 
     def to_hdf5_group(self, group: h5py.Group, compress: bool = True) -> None:
         group.attrs["_type"] = "BootstrapMeta"
 
@@ -153,6 +153,11 @@ def _attribute_to_loop_class(
                     + extra_diameter * max_lifetime,
                     cocycle_edge_mask=cocycle_edge_mask,
                     column_scores=column_scores,
+                    embedding=side.embedding,
+                    death_scale=max(
+                        (lc.death for lc in side.loop_classes if lc is not None),
+                        default=None,
+                    ),
                     **kwargs_equivalence,
                 )
                 if equivalence.is_equivalent(relax=with_relaxation):
 
@@ -147,6 +147,8 @@ def find_loops(
         boundary_thresh = threshold_boundary
         if boundary_thresh is None:
             boundary_thresh = threshold_homology
+        meta.bootstrap.threshold_homology = threshold_homology
+        meta.bootstrap.threshold_boundary = threshold_boundary
         hd._compute_boundary_matrix_d1(
             adata=adata, thresh=boundary_thresh, verbose=verbose
         )
Original file line number	Diff line number	Diff line change
`@@ -147,6 +147,8 @@ def find_loops(`
`147`	`147`	`boundary_thresh = threshold_boundary`
`148`	`148`	`if boundary_thresh is None:`
`149`	`149`	`boundary_thresh = threshold_homology`
	`150`	`+ meta.bootstrap.threshold_homology = threshold_homology`
	`151`	`+ meta.bootstrap.threshold_boundary = threshold_boundary`
`150`	`152`	`hd._compute_boundary_matrix_d1(`
`151`	`153`	`adata=adata, thresh=boundary_thresh, verbose=verbose`
`152`	`154`	`)`