cleaning up analysis file

Author: codegithubka

docs/kimon_updates.md

@@ -217,22 +217,17 @@ As a resultl, we can probably use a 1000x1000 grid for our HPC simulation!
 
 ## Testing and HPC Run Update (23/1)
 
-HPC Run Estimate (we are using 32 cores)
+HPC Run Estimate (we are using 32 cores).
 
-bash
-```
-2026-01-22 21:09:55,625 [INFO] ============================================================
-2026-01-22 21:09:55,625 [INFO] PP Evolutionary Analysis - OPTIMIZED VERSION
-2026-01-22 21:09:55,625 [INFO] ============================================================
-2026-01-22 21:09:55,625 [INFO] Mode: full
-2026-01-22 21:09:55,626 [INFO] Output: results
-2026-01-22 21:09:55,626 [INFO] Cores: -1
-2026-01-22 21:09:55,626 [INFO] Numba: ENABLED
-2026-01-22 21:09:55,626 [INFO] Directed hunting: DISABLED
-2026-01-22 21:09:55,626 [INFO] Estimated: 23,000 sims, ~0.1h on 96 cores (~5 core-hours)
-2026-01-22 21:09:55,626 [INFO] Dry run - exiting
-(snellius_venv) [kanagnostopoul@int5 ~]$ 
-```
+1000 x 1000 grid -> 1 million cells
+
+At each step: 500 million operations per simulation
+
+This is multiplied by the number of replicates. 50 reps will result in 22,500 simulations.
+
+By the benchmark, we have 1,182 steps per second (throuput) for a 100x100 grid. If we use a 1000x1000 grid, that implies 11.8 steps/second. So 1000x1000 grid with 50 reps
+
+8.26 hours (not ideal!)
 
 ### Tests
 
@@ -281,15 +276,27 @@ We have 48 tests cases validating the folloiwing:
 - Edge cases with extreme parameter values
 
 
-### Issues to be resolved
+## Issues to be resolved
+
+1. Grid size for HPC run
+2. Number of replicates for statistal power
+3. Directed and/or undirected runs 
+4. Evolving and non-evolving runs?
+5. Mean field baseline or non evolving basiline
+6. Warmup period and measurement steps (i.e how many steps do we need to avoid init bias?)
+7. Measurement frequency for statistical accuracy
+8. Default parameters (Need Storm's input on this one).
+
+Options:
+
+1. Asymmetric repliates for non-evolving runs
+2. Coarse initial parameter sweep grid
+3. Discard non-evo runs and use mean field baseline instead or the opposite
 
-1. Grid size for Snellius run
+NOTE: Without the optimization kernels for a 1000x1000 grid the simulation (using 50 reps for statistical power) would run for 548 hours (approximately 23 days)
 
-2. Number of replicates for statistical power
 
-3. Directed and/or undirected runs
 
-The above decisions will complicate our analysis in terms of noise reduction, finite size scaling, and capturing spatial correlations (hryda effect).
 
 
 

mock_results/analysis.log

@@ -0,0 +1,22 @@
+2026-01-23 11:56:47,627 [INFO] ============================================================
+2026-01-23 11:56:47,627 [INFO] PP Evolutionary Analysis - OPTIMIZED VERSION
+2026-01-23 11:56:47,627 [INFO] ============================================================
+2026-01-23 11:56:47,627 [INFO] Mode: plot
+2026-01-23 11:56:47,627 [INFO] Output: mock_results
+2026-01-23 11:56:47,627 [INFO] Cores: -1
+2026-01-23 11:56:47,627 [INFO] Numba: ENABLED
+2026-01-23 11:56:47,627 [INFO] Directed hunting: DISABLED
+2026-01-23 11:56:47,627 [INFO] Estimated: 7,250 sims, ~0.1h on 8 cores (~1 core-hours)
+2026-01-23 11:56:48,037 [INFO] Loaded 400 results
+2026-01-23 11:56:48,779 [INFO] Saved phase_diagrams.png
+2026-01-23 11:56:49,001 [INFO] Saved hydra_analysis.png
+2026-01-23 11:56:49,226 [INFO] Saved pcf_analysis.png
+2026-01-23 11:56:49,226 [INFO] ============================================================
+2026-01-23 11:56:49,226 [INFO] ANALYSIS SUMMARY
+2026-01-23 11:56:49,226 [INFO] ============================================================
+2026-01-23 11:56:49,226 [INFO] Hydra region: 0 combinations
+2026-01-23 11:56:49,226 [INFO] Max Hydra strength: nan
+2026-01-23 11:56:49,226 [INFO] Mean segregation index: 1.000
+2026-01-23 11:56:49,226 [INFO] Critical point: pb=0.100, pd=0.100
+2026-01-23 11:56:49,227 [INFO] Total runtime: 0.00 hours
+2026-01-23 11:56:49,227 [INFO] Done!

mock_results/config.json

@@ -0,0 +1,47 @@
+{
+  "default_grid": 100,
+  "densities": [
+    0.3,
+    0.15
+  ],
+  "n_prey_birth": 15,
+  "n_prey_death": 15,
+  "prey_birth_min": 0.1,
+  "prey_birth_max": 0.35,
+  "prey_death_min": 0.001,
+  "prey_death_max": 0.1,
+  "predator_death": 0.1,
+  "predator_birth": 0.2,
+  "n_replicates": 15,
+  "warmup_steps": 200.0,
+  "measurement_steps": 300,
+  "cluster_samples": 1,
+  "collect_pcf": true,
+  "pcf_sample_rate": 0.2,
+  "pcf_max_distance": 20.0,
+  "pcf_n_bins": 20,
+  "evolve_sd": 0.1,
+  "evolve_min": 0.001,
+  "evolve_max": 0.1,
+  "fss_grid_sizes": [
+    50,
+    75,
+    100,
+    150
+  ],
+  "fss_replicates": 100,
+  "sensitivity_sd_values": [
+    0.02,
+    0.05,
+    0.1,
+    0.15,
+    0.2
+  ],
+  "sensitivity_replicates": 20,
+  "synchronous": false,
+  "directed_hunting": false,
+  "save_diagnostic_plots": false,
+  "diagnostic_param_sets": 5,
+  "min_analysis_density": 0.002,
+  "n_jobs": -1
+}

mock_results/hydra_analysis.png

@@ -0,0 +1,11 @@
+{
+  "coexistence_no_evo": 0,
+  "hydra_region_size": 0,
+  "max_hydra_strength": NaN,
+  "hydra_region_size_evo": 0,
+  "mean_segregation_index": 1.0,
+  "mean_prey_clustering": 1.5,
+  "critical_prey_birth": 0.1,
+  "critical_prey_death": 0.1,
+  "critical_tau_prey": 2.009731666223138
+}

mock_results/pcf_analysis.png

@@ -0,0 +1,61 @@
+import numpy as np
+from pathlib import Path
+import json
+from dataclasses import asdict
+import sys
+import os
+
+sys.path.append(os.getcwd())
+from scripts.pp_analysis import Config, save_sweep_binary
+
+def generate_mock_data(output_dir: str):
+    cfg = Config(default_grid=100, n_prey_birth=10, n_prey_death=10, n_replicates=2)
+    out_path = Path(output_dir)
+    out_path.mkdir(exist_ok=True)
+
+    prey_births = np.linspace(0.1, 0.4, 10)
+    prey_deaths = np.linspace(0.01, 0.1, 10)
+    results = []
+
+    for pb in prey_births:
+        for pd in prey_deaths:
+            for evo in [False, True]:
+                for rep in range(cfg.n_replicates):
+                    hydra_factor = 5000 * np.exp(-(pd - 0.04)**2 / 0.001) if not evo else 6000
+                    base_pop = (pb * 10000) - (pd * 20000) + hydra_factor
+                    
+                    # Ensure pop doesn't go negative
+                    prey_mean = max(50, base_pop + np.random.normal(0, 100))
+                    
+                    # Mock PCF data
+                    dist = np.linspace(0.5, 20, 20)
+                    # C_cr < 1 indicates segregation
+                    seg_idx = 0.8 + (pd * 2) 
+                    
+                    res = {
+                        "prey_birth": float(pb),
+                        "prey_death": float(pd),
+                        "with_evolution": evo,
+                        "prey_mean": float(prey_mean),
+                        "pred_mean": float(prey_mean * 0.4),
+                        "prey_survived": bool(prey_mean > 100),
+                        "prey_tau": 2.05 + np.random.normal(0, 0.05), 
+                        "evolved_prey_death_mean": float(pd * 1.2) if evo else np.nan,
+                        "segregation_index": float(seg_idx),
+                        "prey_clustering_index": 1.5,
+                        "pcf_distances": dist.tolist(),
+                        "pcf_prey_prey_mean": (1.0 + np.exp(-dist/2)).tolist()
+                    }
+                    results.append(res)
+
+    # Save as .npz to mimic the real output
+    save_sweep_binary(results, out_path / "sweep_results.npz")
+    
+    # Save a mock config
+    with open(out_path / "config.json", "w") as f:
+        json.dump(asdict(cfg), f, indent=2)
+        
+    print(f"Mock data generated in {output_dir}")
+
+if __name__ == "__main__":
+    generate_mock_data("mock_results")

mock_results/phase_diagrams.png

@@ -32,6 +32,7 @@
 from typing import Dict, List, Tuple, Optional
 import warnings
 from tqdm import tqdm
+import hashlib
 
 project_root = str(Path(__file__).parents[1])
 if project_root not in sys.path:
@@ -68,31 +69,33 @@ class Config:
     """Central configuration for analysis."""
 
     # Grid settings
-    default_grid: int = 1000
-    densities: Tuple[float, float] = (0.30, 0.15)
+    default_grid: int = 100 #FIXME: Decide default configuration
+    densities: Tuple[float, float] = (0.30, 0.15) #FIXME: Default densities
 
     # 2D sweep resolution
-    n_prey_birth: int = 15
+    n_prey_birth: int = 15 # FIXME: Decide number of grid points along prey axes
     n_prey_death: int = 15
-    prey_birth_min: float = 0.10
+    prey_birth_min: float = 0.10 # FIXME: Range of prey death to sweep
     prey_birth_max: float = 0.35
     prey_death_min: float = 0.001
     prey_death_max: float = 0.10
 
     # Fixed predator parameters
-    predator_death: float = 0.1
-    predator_birth: float = 0.2
+    predator_death: float = 0.1 # FIXME: Default predator death rate
+    predator_birth: float = 0.2 # FIXME: Default predator birth
 
     # Replicates
-    n_replicates: int = 15
+    n_replicates: int = 15 # FIXME: Decide number of indep. runs per parameter config
 
     # Simulation timing
-    warmup_steps: int = 200
-    measurement_steps: int = 300
+    warmup_steps: int = 200 * (default_grid / 100) # FIXME: Steps to run before measuring
+    measurement_steps: int = 300 # FIXME: Decide measurement steps
 
     # Cluster/PCF sampling
     cluster_samples: int = 1  # Reduced from 3 - PCF is expensive
-    cluster_interval: int = 299  # Sample near end of measurement
+    @property
+    def cluster_interval(self) -> int:
+        return self.measurement_steps - 1 # Sample near end of measurement
 
     # PCF settings
     collect_pcf: bool = True
@@ -101,26 +104,29 @@ class Config:
     pcf_n_bins: int = 20
 
     # Evolution parameters
-    evolve_sd: float = 0.10
+    evolve_sd: float = 0.10 # FIXME: Tune evolution parameters
     evolve_min: float = 0.001
     evolve_max: float = 0.10
 
     # Finite size scaling
-    fss_grid_sizes: Tuple[int, ...] = (50, 75, 100, 150)
+    fss_grid_sizes: Tuple[int, ...] = (50, 75, 100, 150) # FIXME: Grid sizes for FSS
     fss_replicates: int = 100
 
     # Evolution sensitivity analysis
-    sensitivity_sd_values: Tuple[float, ...] = (0.02, 0.05, 0.10, 0.15, 0.20)
+    sensitivity_sd_values: Tuple[float, ...] = (0.02, 0.05, 0.10, 0.15, 0.20) # FIXME: SD values to test
     sensitivity_replicates: int = 20
 
     # Update mode
-    synchronous: bool = False
-    directed_hunting: bool = False 
+    synchronous: bool = False # NOTE: This should always be False for PP model
+    directed_hunting: bool = True # FIXME: With or without directed hunting functionality
 
     # Diagnostic snapshots
     save_diagnostic_plots: bool = False
     diagnostic_param_sets: int = 5
 
+    # Min density required for PCF/Clsuter Analysis
+    min_analysis_density: float = 0.002 # FIXME: Minimum prey density (fraction of grid) to analyze clusters/PCF
+    
     # Parallelization
     n_jobs: int = -1
 
@@ -173,6 +179,15 @@ def estimate_runtime(self, n_cores: int = 32) -> str:
 # HELPER FUNCTIONS
 # =============================================================================
 
+
+def generate_unique_seed(pb: float, pd:float, rep:int) -> int:
+    """Creates a unique, deterministic seed from parameters."""
+    
+    identifier = f"{pb:.6f}_{pd:.6f}_{rep}".encode()
+    hash_hex = hashlib.sha256(identifier).hexdigest()[:8]
+
+    return int(hash_hex, 16)
+
 def count_populations(grid: np.ndarray) -> Tuple[int, int, int]:
     """Count empty, prey, predator cells."""
     return int(np.sum(grid == 0)), int(np.sum(grid == 1)), int(np.sum(grid == 2))
@@ -363,6 +378,9 @@ def run_single_simulation(
 
     sample_counter = 0
 
+    # Calculate threshold based on area
+    min_count = int(cfg.min_analysis_density * (grid_size**2))
+    
     for step in range(cfg.measurement_steps):
         model.update()
         _, prey, pred = count_populations(model.grid)
@@ -376,7 +394,7 @@ def run_single_simulation(
 
         # Cluster and PCF sampling
         if step >= cfg.cluster_interval and sample_counter < cfg.cluster_samples:
-            if prey > 10:
+            if prey >= min_count and pred >= (min_count // 4):
                 prey_clusters.extend(measure_cluster_sizes_fast(model.grid, 1))
                 pred_clusters.extend(measure_cluster_sizes_fast(model.grid, 2))
 
@@ -572,10 +590,10 @@ def run_2d_sweep(cfg: Config, output_dir: Path, logger: logging.Logger) -> List[
     for pb in prey_births:
         for pd in prey_deaths:
             for rep in range(cfg.n_replicates):
-                seed_base = int(pb * 1000) + int(pd * 10000) + rep
+                seed = generate_unique_seed(pb, pd, rep)
                 # Both with and without evolution
-                jobs.append((pb, pd, cfg.default_grid, seed_base, False))
-                jobs.append((pb, pd, cfg.default_grid, seed_base, True))
+                jobs.append((pb, pd, cfg.default_grid, seed, False)) #FIXME: Consider cutting non-evo runs
+                jobs.append((pb, pd, cfg.default_grid, seed, True))
 
     logger.info(f"2D Sweep: {len(jobs):,} simulations")
     logger.info(f"  Grid: {len(prey_births)}×{len(prey_deaths)} parameters")
@@ -939,13 +957,18 @@ def generate_plots(cfg: Config, output_dir: Path, logger: logging.Logger):
 
     ax = axes[2]
     mid_pb_idx = n_pb // 2
-    ax.plot(prey_deaths, grids["prey_pop_no_evo"][:, mid_pb_idx], 'b-o',
-            label=f'No Evo (pb={prey_births[mid_pb_idx]:.2f})', markersize=4)
-    ax.plot(prey_deaths, grids["prey_pop_evo"][:, mid_pb_idx], 'g-s',
-            label=f'With Evo (pb={prey_births[mid_pb_idx]:.2f})', markersize=4)
+    target_pb = prey_births[mid_pb_idx]
+    no_evo_slice = grids["prey_pop_no_evo"][:, mid_pb_idx]
+    evo_slice = grids["prey_pop_evo"][:, mid_pb_idx]
+    
+    ax.plot(prey_deaths, no_evo_slice, 'b-o',
+            label=f'No Evo (pb={target_pb:.2f})', markersize=4)
+    ax.plot(prey_deaths, evo_slice, 'g-s',
+            label=f'With Evo (pb={target_pb:.2f})', markersize=4)
+    
     ax.set_xlabel("Prey Death Rate")
     ax.set_ylabel("Prey Population")
-    ax.set_title("Prey Pop vs Death Rate Slice")
+    ax.set_title(f"Prey Pop vs Death Rate (pb={target_pb:.2f})") # Dynamic title
     ax.legend()
     ax.grid(True, alpha=0.3)