Add min-dwell stability to detect_regimes and regenerate Stage 5 figures

Agent-Logs-Url: https://github.com/dfeen87/Cell-DNA-Dynamics/sessions/73a94cea-919c-44c1-9c6f-e89a944478de

Co-authored-by: dfeen87 <158860247+dfeen87@users.noreply.github.com>

Author: Copilot

figures/stage5_regime_transitions/regime_segmentation.png

@@ -48,11 +48,77 @@
 from src.core.models import RegimeLabels
 
 
+def _enforce_min_dwell(labels: np.ndarray, min_dwell: int) -> np.ndarray:
+    """Merge regime intervals shorter than *min_dwell* samples into neighbours.
+
+    Iteratively finds the first contiguous run whose length is less than
+    *min_dwell* and reassigns it to match the longer of its two adjacent runs
+    (left neighbour wins a tie).  The loop terminates once all runs meet the
+    minimum length requirement.
+
+    Parameters
+    ----------
+    labels : np.ndarray, shape (T,), dtype object
+        Per-time-point regime label array to stabilise in-place copy.
+    min_dwell : int
+        Minimum number of consecutive samples required to keep a regime run.
+        Values ≤ 1 disable the filter.
+
+    Returns
+    -------
+    np.ndarray
+        A new label array with all runs at least *min_dwell* samples long.
+    """
+    if min_dwell <= 1:
+        return labels.copy()
+
+    result = labels.copy()
+
+    while True:
+        # Build run-length encoding: list of [regime, start, stop].
+        runs: list[list] = []
+        cur = result[0]
+        start = 0
+        for i in range(1, len(result)):
+            if result[i] != cur:
+                runs.append([cur, start, i])
+                cur = result[i]
+                start = i
+        runs.append([cur, start, len(result)])
+
+        # Find the first run that is too short.
+        short_idx: int | None = None
+        for k, (_, s, e) in enumerate(runs):
+            if (e - s) < min_dwell:
+                short_idx = k
+                break
+
+        if short_idx is None:
+            break  # All runs satisfy the minimum dwell time.
+
+        k = short_idx
+        left_len = (runs[k - 1][2] - runs[k - 1][1]) if k > 0 else -1
+        right_len = (runs[k + 1][2] - runs[k + 1][1]) if k < len(runs) - 1 else -1
+
+        if left_len >= right_len and k > 0:
+            # Absorb into the left neighbour.
+            result[runs[k][1] : runs[k][2]] = runs[k - 1][0]
+        elif k < len(runs) - 1:
+            # Absorb into the right neighbour.
+            result[runs[k][1] : runs[k][2]] = runs[k + 1][0]
+        else:
+            # Single-run array; nothing to merge.
+            break
+
+    return result
+
+
 def detect_regimes(
     delta_phi: np.ndarray,
     q_low: float = 0.75,
     q_high: float = 0.90,
     detect_recovery: bool = True,
+    min_dwell: int = 10,
 ) -> RegimeLabels:
     """Segment ΔΦ(t) into regime labels following the manuscript's rules.
 
@@ -77,6 +143,14 @@ def detect_regimes(
     the system's return toward the reference stability region (Section 6:
     "cellular recovery trajectories").
 
+    Minimum dwell-time stabilisation (*min_dwell* > 1):
+    After the initial per-sample labelling and optional recovery detection,
+    any contiguous run of identical labels shorter than *min_dwell* samples
+    is merged into its longer adjacent run.  This prevents the classifier
+    from switching regimes on single-sample noise and ensures the
+    segmentation output contains broad, stable regime intervals rather than
+    rapid oscillations around threshold boundaries.
+
     Parameters
     ----------
     delta_phi : array_like, shape (T,)
@@ -93,6 +167,10 @@ def detect_regimes(
         Whether to apply the optional recovery-labelling post-processing
         step.  Set to ``False`` to obtain the three-regime version
         (stable / pre-instability / instability only).
+    min_dwell : int, default 10
+        Minimum number of consecutive samples required to retain a regime
+        run.  Any run shorter than this value is absorbed into its longer
+        neighbour.  Set to 1 to disable minimum-dwell-time stabilisation.
 
     Returns
     -------
@@ -168,6 +246,15 @@ def detect_regimes(
             recovery_mask = labels[post_instability] == "stable"
             labels[post_instability[recovery_mask]] = "recovery"
 
+    # ------------------------------------------------------------------
+    # 4. Minimum dwell-time stabilisation.
+    #
+    #    Merge short regime runs to prevent the classifier from toggling
+    #    on single-sample threshold crossings due to measurement noise.
+    # ------------------------------------------------------------------
+    if min_dwell > 1:
+        labels = _enforce_min_dwell(labels, min_dwell)
+
     return RegimeLabels(
         labels=labels,
         theta_low=theta_low,

src/stage5_regime_transitions/detect_regimes.py

@@ -9,9 +9,9 @@
    - Horizontal dashed lines for θ_low and θ_high thresholds
 
 2. figures/stage5_regime_transitions/regime_segmentation.png
-   - ΔΦ(t) curve over time
-   - Shaded regions for each regime (stable, pre-instability, instability, recovery)
-   - Vertical markers at each regime transition
+   - Broad colored regime bands (stable / pre-instability / instability / recovery)
+   - Clean vertical boundary lines at regime transitions
+   - ΔΦ(t) overlaid as a single curve
 
 Reproducibility settings:
     seed = 0, DPI = 150, figure sizes as specified per figure below.
@@ -47,13 +47,14 @@
 DPI = 150
 np.random.seed(SEED)
 
-# ── Consistent regime colour palette ─────────────────────────────────────────
+# ── Consistent regime colour palette (Stage 2 / Stage 3 compatible) ──────────
 REGIME_COLORS = {
     "stable":          "#2196F3",   # blue
     "pre-instability": "#FF9800",   # orange
     "instability":     "#F44336",   # red
     "recovery":        "#4CAF50",   # green
 }
+REGIME_ALPHA = 0.18
 
 # ── Paths ─────────────────────────────────────────────────────────────────────
 _DATA_PATH = os.path.join(_REPO_ROOT, "data", "synthetic", "eic_timeseries.csv")
@@ -125,46 +126,45 @@ def plot_regime_segmentation(
     regime_labels,
     out_path: str,
 ) -> None:
-    """Plot ΔΦ(t) with shaded regime regions and transition markers and save."""
+    """Plot broad regime bands with ΔΦ(t) overlaid and transition markers."""
     intervals = regime_shading_intervals(regime_labels)
     transitions = extract_transition_events(regime_labels)
 
     fig, ax = plt.subplots(figsize=(14, 5))
 
-    # ── Shaded regime spans ───────────────────────────────────────────────────
+    # ── Broad colored regime bands ────────────────────────────────────────────
+    dt_half = (t[1] - t[0]) / 2
     for iv in intervals:
         t_start = t[iv.start]
-        # Use half-step extension at the right edge to avoid gaps
-        t_end = t[iv.stop - 1] + (t[1] - t[0]) / 2 if iv.stop < len(t) else t[-1]
+        t_end = t[iv.stop - 1] + dt_half if iv.stop < len(t) else t[-1] + dt_half
         ax.axvspan(
             t_start,
             t_end,
-            alpha=0.18,
+            alpha=REGIME_ALPHA,
             color=REGIME_COLORS.get(iv.regime, "#CCCCCC"),
             linewidth=0,
         )
 
-    # ── ΔΦ(t) curve ──────────────────────────────────────────────────────────
+    # ── Clean boundary lines at transitions ──────────────────────────────────
+    for ev in transitions:
+        ax.axvline(
+            t[ev.index],
+            color="#616161",
+            linewidth=0.8,
+            linestyle="--",
+            alpha=0.6,
+            zorder=3,
+        )
+
+    # ── ΔΦ(t) overlaid as a single curve ─────────────────────────────────────
     ax.plot(
         t, delta_phi,
         color="#212121",
-        linewidth=0.9,
-        zorder=3,
+        linewidth=1.0,
+        zorder=4,
         label=r"$\Delta\Phi(t)$",
     )
 
-    # ── Transition markers ────────────────────────────────────────────────────
-    for ev in transitions:
-        x_trans = t[ev.index]
-        ax.axvline(
-            x_trans,
-            color="#9C27B0",
-            linewidth=0.9,
-            linestyle=":",
-            alpha=0.8,
-            zorder=4,
-        )
-
     ax.set_xlabel("Time (s)", fontsize=11)
     ax.set_ylabel(r"$\Delta\Phi(t)$", fontsize=11)
     ax.set_title(
@@ -173,33 +173,35 @@ def plot_regime_segmentation(
     )
 
     # ── Legend ────────────────────────────────────────────────────────────────
+    present_regimes = set(regime_labels)
     regime_patches = [
         Patch(
             facecolor=REGIME_COLORS[r],
-            alpha=0.5,
+            alpha=0.6,
             label=r.capitalize(),
         )
         for r in ["stable", "pre-instability", "instability", "recovery"]
-        if r in set(regime_labels)
+        if r in present_regimes
     ]
     transition_handle = Line2D(
         [0], [0],
-        color="#9C27B0",
-        linewidth=0.9,
-        linestyle=":",
-        alpha=0.8,
+        color="#616161",
+        linewidth=0.8,
+        linestyle="--",
+        alpha=0.6,
         label="Transition",
     )
     delta_phi_handle = Line2D(
         [0], [0],
         color="#212121",
-        linewidth=0.9,
+        linewidth=1.0,
         label=r"$\Delta\Phi(t)$",
     )
     ax.legend(
         handles=[delta_phi_handle] + regime_patches + [transition_handle],
         loc="upper right",
         fontsize=8,
+        framealpha=0.9,
     )
     ax.grid(True, linestyle="--", linewidth=0.4, alpha=0.4)
 
@@ -220,8 +222,8 @@ def main() -> None:
     # ── 2. Compute ΔΦ(t) ─────────────────────────────────────────────────────
     delta_phi = compute_delta_phi(series)
 
-    # ── 3. Detect regimes ─────────────────────────────────────────────────────
-    regime_result = detect_regimes(delta_phi)
+    # ── 3. Detect regimes with minimum dwell-time stability (10 samples) ──────
+    regime_result = detect_regimes(delta_phi, min_dwell=10)
 
     # ── 4. Plot ΔΦ(t) with threshold lines ───────────────────────────────────
     plot_delta_phi_thresholds(