Kimon

.gitignore

@@ -12,4 +12,12 @@ results/
 data/
 .pytest_cache/
 
-.DS_Store
+.DS_Store
+
+hpc_data/phase4_18735304
+
+hpc_data/phase6_18780164
+
+
+hpc_data/**/*.jsonl
+hpc_data/**/*.jsonl

README.md

@@ -0,0 +1,296 @@
+## Predator-Prey Cellular Automaton: Model Documentation
+
+### Overview
+
+This project impelments a spatial predator-prey cellular automaton (CA) to investigate the Hydra effect, a counterintuitive ecological phenomenon where increased prey mortality paradoxically leads to higher prey population densities. The model explores the rise of emergent population dynamics through spatial structure and local interactions that differe fundemntally from well-mixed (mean-field) predictions.
+
+The codebase uses Numba JIT compilation for computationally intensive kernels and is tuned for high performance computing environments. It supporrs parameter sweeps, finite size scaling analysis, and self-organized criticality exploration.
+
+---
+
+### Background
+
+#### The Hydra Effect
+
+In classical Lotka-Volterra dynamics, increasing prey mortality always reduces the prey population. However, theoretical work has identified conditions where the opposite effect can be observed. The Hydra effect emerges in spatially structured systems where
+
+1. Local interactions prevent predators from immediately exploiting all available prey.
+2. Spatial refugia from naturally through predator-prey dynamics
+3. Increased prey mortality can disrupt predator populations disproportionally.
+
+This study uses a cellular automaton framework toi study how spatial strcuture generates the Hydra effecr and whether the system exhibits signatures of self-organized criticality at the transition point.
+
+#### Self-Organized Criticality (SOC)
+
+SOC refers to systems that naturally evolve toward a critical state without external tuning. At criticality, such systems exhibit:
+
+- Power-law cluster size distributions: $P(s) \sim s^{-\tau}$
+- Scale-free correlations: No characteristic length scale dominates
+- Critical slowing down: Perturbations decay slowly near the critical point
+- Fractal spatial patterns: Self-similar structure across scales
+
+In the predator-prey context, SOC would manifest as the system self-tuning toward a critical prey mortality rate where population dynamics become scale-free and the Hydra effect is maximized.
+
+---
+
+### Model Description
+
+The model uses a 2D lattice with periodic boundary conditions. Each cell occupies one of three states:
+
+| State | Value | Description |
+|-------|-------|-------------|
+| Empty | 0 | Unoccupied cell |
+| Prey | 1 | Prey organism |
+| Predator | 2 | Predator organism |
+
+The model uses asynchronous updates: cells are processed in random order each timestep, with state changes taking effect immediately. This prevents artificial synchronization artifacts common in parallel update schemes.
+
+For each prey cell, in order:
+
+- Death: With probability `prey_death`, the prey dies and the cell becomes empty
+- Reproduction: If alive and a randomly selected neighbor is empty, with probability `prey_birth`, a new prey is placed in that neighbor cell
+
+For each predator cell, in order:
+
+- Death: With probability `predator_death`, the predator dies (starvation)
+- Hunting/Reproduction: If alive and a randomly selected neighbor contains prey, with probability `predator_birth`, the predator consumes the prey and reproduces into that cell
+
+The model supports both neighborhood types:
+
+- Moore neighborhood: (default): 8 adjacent cells (including diagonals)
+- von Neumann neighborhood: 4 adjacent cells (cardinal directions only)
+
+Moore neighborhoods are used throughout the experiments as they provide more realistic movement and interaction patterns.
+
+---
+
+### Hunting Modes
+
+The model implements two distinct neighbor selection strategies that qualitatively affect dynamics.
+
+#### Random Neighbor Selection (Default)
+
+In the standard mode, each organism selects a single random neighbor for interaction:
+
+```python
+# Prey: pick random neighbor, reproduce if empty
+neighbor_idx = np.random.randint(0, n_neighbors)
+if grid[neighbor] == EMPTY and random() < prey_birth:
+    grid[neighbor] = PREY
+
+# Predator: pick random neighbor, hunt if prey
+if grid[neighbor] == PREY and random() < predator_birth:
+    grid[neighbor] = PREDATOR
+```
+This creates a blind interaction model where organisms are not aware of their surroundings.
+
+#### Directed Hunting Mode
+
+The directed mode implements "intelligent" neighbor selection:
+
+- Prey: Scan all neighbors for empty cells, then randomly select one empty cell for reproduction
+- Predators: Scan all neighbors for prey, then randomly select one prey cell to hunt
+
+```python
+# Directed prey reproduction
+empty_neighbors = [n for n in neighbors if grid[n] == EMPTY]
+if empty_neighbors and random() < prey_birth:
+    target = random.choice(empty_neighbors)
+    grid[target] = PREY
+
+# Directed predator hunting
+prey_neighbors = [n for n in neighbors if grid[n] == PREY]
+if prey_neighbors and random() < predator_birth:
+    target = random.choice(prey_neighbors)
+    grid[target] = PREDATOR
+```
+
+This increases the effective reproduction and predation rates without changing the nominal probability parameters.
+
+---
+
+### Implementation Architecture
+
+#### Class Hierarchy
+
+```
+CA (base class)
+ └── PP (predator-prey specialization)
+      └── PPKernel (Numba-optimized update kernel)
+```
+
+The `CA` base class provides:
+- Grid initialization with specified densities
+- Neighborhood definitions (Moore/von Neumann)
+- Per-cell parameter evolution framework
+- Validation and run loop infrastructure
+
+The `PP` class adds:
+- Predator-prey specific parameters and defaults
+- Synchronous/asynchronous update dispatch
+- Integration with Numba-optimized kernels
+
+#### Basic Usage
+
+```python
+from models.CA import PP
+
+# Initialize model
+model = PP(
+    rows=100,
+    cols=100,
+    densities=(0.30, 0.15),  # (prey_fraction, predator_fraction)
+    params={
+        "prey_birth": 0.2,
+        "prey_death": 0.05,
+        "predator_birth": 0.8,
+        "predator_death": 0.05,
+    },
+    seed=42,
+    synchronous=False,  # Always False for this study
+    directed_hunting=False,
+)
+
+# Run simulation
+for step in range(1000):
+    model.update()
+
+# Access grid state
+prey_count = np.sum(model.grid == 1)
+pred_count = np.sum(model.grid == 2)
+```
+
+#### Evolution Mode
+
+The model supports per-cell parameter evolution, where offspring inherit (with mutation) their parent's parameter values:
+
+```python
+# Enable evolution of prey_death parameter
+model.evolve(
+    "prey_death",
+    sd=0.01,      # Mutation standard deviation
+    min_val=0.0,  # Minimum allowed value
+    max_val=0.20, # Maximum allowed value
+)
+```
+
+When evolution is enabled, each prey cell maintains its own `prey_death` value. Upon reproduction, offspring inherit the parent's value plus Gaussian noise. This allows investigation of whether the population self-organizes toward a critical mortality rate.
+
+---
+
+### Numba Optimization
+
+The computational bottleneck is the update kernel, which must process every occupied cell each timestep. For a 1000×1000 grid with 50% occupancy, this means ~500,000 cell updates per step.
+
+Key optimizations:
+
+- **JIT Compilation**: Core kernels use `@njit(cache=True)` for ahead-of-time compilation
+2. **Pre-allocated Buffers**: The `PPKernel` class maintains reusable arrays to avoid allocation overhead
+3. **Efficient Shuffling**: Fisher-Yates shuffle implemented in Numba for random cell ordering
+4. **Cell Lists for PCF**: Pair correlation functions use spatial hashing for O(N) instead of O(N²) complexity
+
+#### PPKernel Class
+
+```python
+class PPKernel:
+    """Wrapper for predator-prey kernel with pre-allocated buffers."""
+
+    def __init__(self, rows, cols, neighborhood="moore", directed_hunting=False):
+        self.rows = rows
+        self.cols = cols
+        self.directed_hunting = directed_hunting
+        self._occupied_buffer = np.empty((rows * cols, 2), dtype=np.int32)
+
+        # Neighbor offset arrays
+        if neighborhood == "moore":
+            self._dr = np.array([-1, -1, -1, 0, 0, 1, 1, 1], dtype=np.int32)
+            self._dc = np.array([-1, 0, 1, -1, 1, -1, 0, 1], dtype=np.int32)
+        else:
+            self._dr = np.array([-1, 1, 0, 0], dtype=np.int32)
+            self._dc = np.array([0, 0, -1, 1], dtype=np.int32)
+```
+---
+
+### Cluster Detection
+
+Clusters are contiguous groups of same-species cells (using Moore connectivity). The implementation uses stack-based flood fill with periodic boundary conditions.
+
+```python
+from models.numba_optimized import get_cluster_stats_fast
+
+stats = get_cluster_stats_fast(grid, species=1)  # Prey clusters
+print(f"Number of clusters: {stats['n_clusters']}")
+print(f"Largest cluster: {stats['largest']} cells")
+print(f"Largest fraction: {stats['largest_fraction']:.3f}")
+```
+
+Metrics collected:
+- Cluster size distribution: Power-law indicates criticality
+- Largest cluster fraction: Order parameter for percolation transition
+---
+
+### Configuration System
+
+The `Config` dataclass centralizes all experimental parameters:
+
+```python
+from config import PHASE1_CONFIG, Config
+
+# Use predefined phase config
+cfg = PHASE1_CONFIG
+
+# Or create custom config
+cfg = Config(
+    grid_size=500,
+    n_replicates=20,
+    warmup_steps=1000,
+    measurement_steps=2000,
+    collect_pcf=True,
+    pcf_sample_rate=0.5,
+)
+
+# Access parameters
+print(f"Grid: {cfg.grid_size}x{cfg.grid_size}")
+print(f"Estimate: {cfg.estimate_runtime(n_cores=32)}")
+```
+
+Each phase has a predefined configuration (`PHASE1_CONFIG` through `PHASE5_CONFIG`) with appropriate defaults for that analysis.
+
+---
+
+### Output Format
+
+Results are saved in JSONL format (one JSON object per line) for efficient streaming and parallel processing:
+
+```json
+{"prey_birth": 0.2, "prey_death": 0.05, "seed": 12345, "final_prey": 28500, ...}
+{"prey_birth": 0.2, "prey_death": 0.06, "seed": 12346, "final_prey": 27200, ...}
+```
+
+Each result dictionary contains:
+- Input parameters
+- Final population counts
+- Cluster statistics
+- Evolution statistics (if enabled)
+- Time series (if enabled)
+
+Metadata files (`phase{N}_metadata.json`) accompany each results file with configuration details and runtime information.
+
+---
+
+### Dependencies
+
+**Required:**
+- Python 3.8+
+- NumPy
+- Numba (for JIT compilation)
+- tqdm (progress bars)
+- joblib (parallelization)
+
+**Optional:**
+- matplotlib (visualization)
+- scipy (additional analysis)
+
+---
+
+## References