stability-shuffle null fix / preprocess reorder

CHANGELOG.md

@@ -5,6 +5,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 While the major version is `0.x`, breaking API changes may land in any release;
 treat the package as alpha-stage until `1.0`.
 
+## [0.1.2] — 2026-05-06
+
+### Fixed
+
+- Stability shuffle no longer injects `0.0` for degenerate iterations;
+  the p-value denominator now uses the count of valid shuffles
+  (Phipson & Smyth 2010). Previously anti-conservative for stable units
+  on sparse or short sessions.
+
+### Changed
+
+- Arena `preprocess_behavior` step order: `perspective → Hampel → clip
+  → mm` (was Hampel first). Bundle outputs unchanged for existing data.
+- `oasis-deconv` `ImportError` reformatted with install commands on
+  separate lines.
+
 ## [0.1.1] — 2026-05-06
 
 ### Fixed

camap/analysis/spatial_1d.py

@@ -353,7 +353,8 @@ def compute_stability_score_1d(
     (sensitive to session-long drift); ``n_split_blocks>=4`` interleaves
     the halves so within-session drift averages out. Combining both
     (e.g. ``stability_splits=[2, 10]``) tests stability at two different
-    timescales.
+    timescales. Degenerate shuffle iterations are skipped rather than
+    counted as 0 (Phipson & Smyth 2010); see :func:`compute_stability_score`.
 
     Returns
     -------
@@ -492,8 +493,8 @@ def compute_half_rate_map(
         if min_shift_frames >= n_frames // 2:
             min_shift_frames = 0
 
-        shuffled_corrs = np.empty(n_shuffles)
-        for i in range(n_shuffles):
+        valid_corrs: list[float] = []
+        for _ in range(n_shuffles):
             shifted = np.roll(aligned_events, _draw_shift(rng, n_frames, min_shift_frames))
             rm1 = _shuffled_rate_map_1d(
                 traj_pos_first,
@@ -515,16 +516,19 @@ def compute_half_rate_map(
             )
 
             if not np.any(bv):
-                shuffled_corrs[i] = 0.0
                 continue
             v1, v2 = rm1[bv], rm2[bv]
             fm = np.isfinite(v1) & np.isfinite(v2)
             if np.sum(fm) < 3:
-                shuffled_corrs[i] = 0.0
                 continue
-            shuffled_corrs[i] = np.corrcoef(v1[fm], v2[fm])[0, 1]
+            valid_corrs.append(float(np.corrcoef(v1[fm], v2[fm])[0, 1]))
 
-        stability_p_val = float((np.sum(shuffled_corrs >= corr) + 1) / (n_shuffles + 1))
+        shuffled_corrs = np.array(valid_corrs)
+        n_valid = len(shuffled_corrs)
+        if n_valid == 0:
+            stability_p_val = np.nan
+        else:
+            stability_p_val = float((np.sum(shuffled_corrs >= corr) + 1) / (n_valid + 1))
     else:
         stability_p_val = np.nan
         shuffled_corrs = np.array([])

camap/analysis/spatial_2d.py

@@ -475,7 +475,10 @@ def compute_stability_score(
 
     Optionally runs a shuffle significance test: circularly shifts
     events and computes the split-half correlation for each shuffle to
-    build a null distribution.
+    build a null distribution. Degenerate iterations (no overlapping
+    valid bins, or fewer than 3 jointly-finite bins) are skipped rather
+    than counted as 0; the p-value uses the count of valid shuffles in
+    its denominator (Phipson & Smyth 2010).
 
     Parameters
     ----------
@@ -657,8 +660,8 @@ def compute_half_rate_map(
         if min_shift_frames >= n_frames // 2:
             min_shift_frames = 0
 
-        shuffled_corrs = np.empty(n_shuffles)
-        for i in range(n_shuffles):
+        valid_corrs: list[float] = []
+        for _ in range(n_shuffles):
             shifted = np.roll(aligned_events, _draw_shift(rng, n_frames, min_shift_frames))
             rm1 = _shuffled_rate_map(
                 traj_x_first,
@@ -682,16 +685,19 @@ def compute_half_rate_map(
             )
 
             if not np.any(bv):
-                shuffled_corrs[i] = 0.0
                 continue
             v1, v2 = rm1[bv], rm2[bv]
             fm = np.isfinite(v1) & np.isfinite(v2)
             if np.sum(fm) < 3:
-                shuffled_corrs[i] = 0.0
                 continue
-            shuffled_corrs[i] = np.corrcoef(v1[fm], v2[fm])[0, 1]
+            valid_corrs.append(float(np.corrcoef(v1[fm], v2[fm])[0, 1]))
 
-        stability_p_val = float((np.sum(shuffled_corrs >= corr) + 1) / (n_shuffles + 1))
+        shuffled_corrs = np.array(valid_corrs)
+        n_valid = len(shuffled_corrs)
+        if n_valid == 0:
+            stability_p_val = np.nan
+        else:
+            stability_p_val = float((np.sum(shuffled_corrs >= corr) + 1) / (n_valid + 1))
     else:
         stability_p_val = np.nan
         shuffled_corrs = np.array([])

camap/dataset/arena.py

@@ -100,8 +100,12 @@ def preprocess_behavior(self) -> None:
         """Apply geometric corrections to the behavior trajectory.
 
         When ``arena_bounds`` is configured:
-            jump removal → perspective correction → boundary clipping →
-            unit conversion (px → mm).
+            perspective correction → Hampel jump removal → boundary
+            clipping → unit conversion (px → mm).
+
+        Hampel runs after perspective so its threshold reflects the
+        post-correction trajectory; it stays before clipping so genuine
+        out-of-bounds outliers remain visible to the filter.
 
         When ``arena_bounds`` is **not** configured:
             warnings are logged and the trajectory remains in pixels.
@@ -142,21 +146,7 @@ def preprocess_behavior(self) -> None:
             self._preprocess_steps: dict[str, pd.DataFrame] = {}
             self._preprocess_steps["Raw"] = self.trajectory[["x", "y"]].copy()
 
-            # 1. Hampel-filter outlier removal on the raw 2D trajectory.
-            self.trajectory, n_jumps = remove_position_jumps(
-                self.trajectory,
-                window_frames=bcfg.hampel_window_frames,
-                n_sigmas=bcfg.hampel_n_sigmas,
-            )
-            logger.info(
-                "Hampel jump removal: %d frames interpolated (window=%d, n_sigmas=%.1f)",
-                n_jumps,
-                bcfg.hampel_window_frames,
-                bcfg.hampel_n_sigmas,
-            )
-            self._preprocess_steps["Jump removal"] = self.trajectory[["x", "y"]].copy()
-
-            # 2. Perspective correction
+            # 1. Perspective correction
             self.trajectory = correct_perspective(
                 self.trajectory,
                 arena_bounds=dbcfg.arena_bounds,
@@ -172,6 +162,22 @@ def preprocess_behavior(self) -> None:
             )
             self._preprocess_steps["Perspective"] = self.trajectory[["x", "y"]].copy()
 
+            # 2. Hampel jump removal — runs on perspective-corrected positions
+            # so the threshold reflects the post-correction trajectory; before
+            # clipping so genuine out-of-bounds outliers remain visible.
+            self.trajectory, n_jumps = remove_position_jumps(
+                self.trajectory,
+                window_frames=bcfg.hampel_window_frames,
+                n_sigmas=bcfg.hampel_n_sigmas,
+            )
+            logger.info(
+                "Hampel jump removal: %d frames interpolated (window=%d, n_sigmas=%.1f)",
+                n_jumps,
+                bcfg.hampel_window_frames,
+                bcfg.hampel_n_sigmas,
+            )
+            self._preprocess_steps["Jump removal"] = self.trajectory[["x", "y"]].copy()
+
             # 3. Boundary clipping
             self.trajectory = clip_to_arena(self.trajectory, arena_bounds=dbcfg.arena_bounds)
             logger.info("Boundary clipping to arena_bounds=%s", dbcfg.arena_bounds)