Removed CV and Auttocorrelation method, replaced with new method of finding equilibrium when upward or downward trending stops

Author: Abou1308

scripts/warmup_study.py

@@ -89,80 +89,69 @@ class WarmupStudyConfig:
 # EQUILIBRATION DETECTION
 # =============================================================================
 
-def estimate_equilibration_cv(
+def estimate_equilibration_trend(
     time_series: np.ndarray,
     sample_interval: int,
     window: int = 10,
     threshold: float = 0.02,
 ) -> int:
     """
-    Estimate equilibration step using coefficient of variation (CV).
-    
-    Returns the step when rolling CV drops below threshold for a sustained period.
+    Estimate equilibration step by detecting when the mean stops trending.
+    
+    This method checks if consecutive rolling means are stable (not changing
+    significantly). Unlike CV-based detection, this is robust to grid size
+    because it measures actual trend, not fluctuation magnitude.
+    
+    Parameters
+    ----------
+    time_series : np.ndarray
+        Population density or count over time.
+    sample_interval : int
+        Number of simulation steps between samples.
+    window : int
+        Size of rolling window for computing means.
+    threshold : float
+        Maximum allowed relative change between consecutive window means.
+        Default 0.02 means the mean can change by at most 2% between windows.
+    
+    Returns
+    -------
+    int
+        Estimated equilibration step.
     """
-    if len(time_series) < window * 2:
+    if len(time_series) < window * 3:
         return len(time_series) * sample_interval
 
     # Skip initial transient (first 10% of data)
-    start_idx = max(window, len(time_series) // 10)
+    start_idx = max(window * 2, len(time_series) // 10)
 
-    for i in range(start_idx, len(time_series)):
-        recent = time_series[i-window:i]
-        mean_val = np.mean(recent)
-        
-        if mean_val > 0:
-            cv = np.std(recent) / mean_val
-            if cv < threshold:
-                # Check if it stays stable for another window
-                if i + window < len(time_series):
-                    next_window = time_series[i:i+window]
-                    next_cv = np.std(next_window) / np.mean(next_window) if np.mean(next_window) > 0 else 1.0
-                    if next_cv < threshold * 1.5:
-                        return i * sample_interval
-                else:
-                    return i * sample_interval
+    # Number of consecutive stable windows required
+    required_stable = 3
+    stable_count = 0
 
-    return len(time_series) * sample_interval
-
-
-def estimate_equilibration_autocorr(
-    time_series: np.ndarray,
-    sample_interval: int,
-    lag_threshold: float = 0.1,
-) -> Tuple[int, float]:
-    """
-    Estimate equilibration via autocorrelation decay.
-    
-    Returns (equilibration_step, integrated_autocorrelation_time).
-    """
-    n = len(time_series)
-    if n < 10:
-        return n * sample_interval, 0.0
-    
-    mean = np.mean(time_series)
-    var = np.var(time_series)
-    
-    if var == 0:
-        return 0, 0.0
-    
-    # Compute autocorrelation function
-    acf = np.correlate(time_series - mean, time_series - mean, mode='full')
-    acf = acf[n-1:] / (var * n)
-    
-    # Integrated autocorrelation time
-    tau_int = 1.0
-    for lag in range(1, min(n // 4, 100)):
-        if acf[lag] > 0:
-            tau_int += 2 * acf[lag]
+    for i in range(start_idx, len(time_series) - window):
+        # Compare means of two consecutive windows
+        window1 = time_series[i - window:i]
+        window2 = time_series[i:i + window]
+        
+        mean1 = np.mean(window1)
+        mean2 = np.mean(window2)
+        
+        if mean1 > 0 and mean2 > 0:
+            # Relative change in mean
+            relative_change = abs(mean2 - mean1) / mean1
+            
+            if relative_change < threshold:
+                stable_count += 1
+                if stable_count >= required_stable:
+                    # Return the step where stability began
+                    return (i - required_stable * window // 2) * sample_interval
+            else:
+                stable_count = 0  # Reset if trend detected
         else:
-            break
+            stable_count = 0  # Reset if population crashed
 
-    # Find where ACF drops below threshold
-    for lag, corr in enumerate(acf):
-        if corr < lag_threshold:
-            return lag * sample_interval, tau_int
-    
-    return len(acf) * sample_interval, tau_int
+    return len(time_series) * sample_interval
 
 
 # =============================================================================
@@ -200,9 +189,7 @@ def set_numba_seed(seed): pass
 
         size_results = {
             'time_per_step': [],
-            'equilibration_steps_cv': [],
-            'equilibration_steps_acf': [],
-            'integrated_autocorr_time': [],
+            'equilibration_steps': [],
             'final_prey_density': [],
             'final_pred_density': [],
         }
@@ -250,43 +237,33 @@ def set_numba_seed(seed): pass
             prey_densities = np.array(prey_densities)
             pred_densities = np.array(pred_densities)
 
-            # Estimate equilibration using both methods
-            eq_steps_cv = estimate_equilibration_cv(
+            # Estimate equilibration (trend-based, robust to grid size)
+            eq_steps = estimate_equilibration_trend(
                 prey_densities,
                 cfg.sample_interval,
                 window=cfg.equilibration_window,
                 threshold=cfg.equilibration_threshold,
             )
 
-            eq_steps_acf, tau_int = estimate_equilibration_autocorr(
-                prey_densities,
-                cfg.sample_interval,
-            )
-            
             size_results['time_per_step'].append(time_per_step)
-            size_results['equilibration_steps_cv'].append(eq_steps_cv)
-            size_results['equilibration_steps_acf'].append(eq_steps_acf)
-            size_results['integrated_autocorr_time'].append(tau_int)
+            size_results['equilibration_steps'].append(eq_steps)
             size_results['final_prey_density'].append(prey_densities[-1])
             size_results['final_pred_density'].append(pred_densities[-1])
 
             if (rep + 1) % max(1, cfg.n_replicates // 5) == 0:
                 logger.info(f"  Replicate {rep+1}/{cfg.n_replicates}: "
-                           f"eq_steps={eq_steps_cv}, time/step={time_per_step*1000:.2f}ms")
+                           f"eq_steps={eq_steps}, time/step={time_per_step*1000:.2f}ms")
 
         # Aggregate results
         results[L] = {
             'grid_size': L,
             'grid_cells': L * L,
             'mean_time_per_step': float(np.mean(size_results['time_per_step'])),
             'std_time_per_step': float(np.std(size_results['time_per_step'])),
-            'mean_eq_steps_cv': float(np.mean(size_results['equilibration_steps_cv'])),
-            'std_eq_steps_cv': float(np.std(size_results['equilibration_steps_cv'])),
-            'mean_eq_steps_acf': float(np.mean(size_results['equilibration_steps_acf'])),
-            'std_eq_steps_acf': float(np.std(size_results['equilibration_steps_acf'])),
-            'mean_tau_int': float(np.mean(size_results['integrated_autocorr_time'])),
+            'mean_eq_steps': float(np.mean(size_results['equilibration_steps'])),
+            'std_eq_steps': float(np.std(size_results['equilibration_steps'])),
             'mean_total_warmup_time': float(
-                np.mean(size_results['equilibration_steps_cv']) * 
+                np.mean(size_results['equilibration_steps']) * 
                 np.mean(size_results['time_per_step'])
             ),
             'mean_final_prey_density': float(np.mean(size_results['final_prey_density'])),
@@ -297,8 +274,8 @@ def set_numba_seed(seed): pass
         logger.info(f"\n  Summary for L={L}:")
         logger.info(f"    Time per step: {results[L]['mean_time_per_step']*1000:.2f} ± "
                    f"{results[L]['std_time_per_step']*1000:.2f} ms")
-        logger.info(f"    Equilibration (CV): {results[L]['mean_eq_steps_cv']:.0f} ± "
-                   f"{results[L]['std_eq_steps_cv']:.0f} steps")
+        logger.info(f"    Equilibration steps: {results[L]['mean_eq_steps']:.0f} ± "
+                   f"{results[L]['std_eq_steps']:.0f}")
         logger.info(f"    Total warmup time: {results[L]['mean_total_warmup_time']:.2f} s")
 
     return results
@@ -342,17 +319,11 @@ def plot_warmup_scaling(
 
     # Panel 2: Equilibration steps vs L (log-log)
     ax = axes[1]
-    eq_steps = [results[L]['mean_eq_steps_cv'] for L in sizes]
-    eq_stds = [results[L]['std_eq_steps_cv'] for L in sizes]
 
+    eq_steps = [results[L]['mean_eq_steps'] for L in sizes]
+    eq_stds = [results[L]['std_eq_steps'] for L in sizes]
     ax.errorbar(sizes, eq_steps, yerr=eq_stds, fmt='o-', capsize=5, 
-                linewidth=2, color='forestgreen', markersize=8, label='CV method')
-    
-    # Also plot ACF method
-    eq_steps_acf = [results[L]['mean_eq_steps_acf'] for L in sizes]
-    eq_stds_acf = [results[L]['std_eq_steps_acf'] for L in sizes]
-    ax.errorbar(sizes, eq_steps_acf, yerr=eq_stds_acf, fmt='s--', capsize=5,
-                linewidth=2, color='darkorange', markersize=6, alpha=0.7, label='ACF method')
+                linewidth=2, color='forestgreen', markersize=8)
 
     ax.set_xscale('log')
     ax.set_yscale('log')
@@ -426,9 +397,10 @@ def plot_scaling_summary(
     # Plot equilibration steps normalized by theoretical scaling
     # Try different z values
     for z, color, style in [(1.0, 'green', '--'), (1.5, 'orange', '-.'), (2.0, 'red', ':')]:
-        eq_normalized = [results[L]['mean_eq_steps_cv'] / (L**z) for L in sizes]
+        eq_normalized = [results[L]['mean_eq_steps'] / (L**z) for L in sizes]
         # Normalize to first point for comparison
-        eq_normalized = [x / eq_normalized[0] for x in eq_normalized]
+        if eq_normalized[0] > 0:
+            eq_normalized = [x / eq_normalized[0] for x in eq_normalized]
         ax.plot(sizes, eq_normalized, style, color=color, linewidth=2, alpha=0.7,
                 label=f'Eq. steps / L^{z:.1f} (normalized)')
 
@@ -554,15 +526,26 @@ def main():
 
     # Compute scaling exponents
     if len(sizes) >= 2:
-        eq_steps = [results[L]['mean_eq_steps_cv'] for L in sizes]
+        eq_steps = [results[L]['mean_eq_steps'] for L in sizes]
         total_times = [results[L]['mean_total_warmup_time'] for L in sizes]
 
-        log_L = np.log(sizes)
-        log_eq = np.log(eq_steps)
-        log_T = np.log(total_times)
+        # Filter out any zero or negative values for log
+        valid_eq = [(L, eq) for L, eq in zip(sizes, eq_steps) if eq > 0]
+        valid_T = [(L, T) for L, T in zip(sizes, total_times) if T > 0]
 
-        z_eq, _, r_eq, _, _ = linregress(log_L, log_eq)
-        z_total, _, r_total, _, _ = linregress(log_L, log_T)
+        if len(valid_eq) >= 2:
+            log_L_eq = np.log([x[0] for x in valid_eq])
+            log_eq = np.log([x[1] for x in valid_eq])
+            z_eq, _, r_eq, _, _ = linregress(log_L_eq, log_eq)
+        else:
+            z_eq, r_eq = 0, 0
+        
+        if len(valid_T) >= 2:
+            log_L_T = np.log([x[0] for x in valid_T])
+            log_T = np.log([x[1] for x in valid_T])
+            z_total, _, r_total, _, _ = linregress(log_L_T, log_T)
+        else:
+            z_total, r_total = 0, 0
 
         logger.info(f"Equilibration steps scaling: t_eq ~ L^{z_eq:.2f} (R² = {r_eq**2:.3f})")
         logger.info(f"Total warmup time scaling: T_warmup ~ L^{z_total:.2f} (R² = {r_total**2:.3f})")