diff --git a/CHANGELOG.md b/CHANGELOG.md index b9b342a..74039a9 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/camap/analysis/spatial_1d.py b/camap/analysis/spatial_1d.py index f0872b5..3e97948 100644 --- a/camap/analysis/spatial_1d.py +++ b/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([]) diff --git a/camap/analysis/spatial_2d.py b/camap/analysis/spatial_2d.py index 92ffe7d..618609f 100644 --- a/camap/analysis/spatial_2d.py +++ b/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([]) diff --git a/camap/dataset/arena.py b/camap/dataset/arena.py index e74f3e1..b2cdcfc 100644 --- a/camap/dataset/arena.py +++ b/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)