Merge branch 'main' into kimon

Author: codegithubka

SOC_ANALYSIS_README.md

@@ -0,0 +1,100 @@
+# Self-Organized Criticality (SOC) Analysis - Summary
+
+## File Created
+**Location:** `scripts/soc_analysis.py`
+
+## Overview
+This comprehensive Python analysis script tests whether your prey-predator cellular automaton exhibits **self-organized criticality** (SOC), with focus on perturbations from initial configurations and diverse parameter sampling.
+
+## Key Features
+
+### 1. **Four SOC Properties Analyzed**
+   - **Slow Drive:** Gradual parameter drift without immediate release
+   - **Stress Build-up:** Interface-based metric tracking potential energy accumulation
+   - **Intermittent Release:** Detection of avalanche cascades in population dynamics
+   - **Self-Organization:** Robustness across diverse parameter combinations
+
+### 2. **Parameter Variations** (Beyond just death/birth rates)
+   - Grid sizes: 16×16 to 64×64
+   - Initial densities: prey (0.1–0.4), predator (0.02–0.15)
+   - Neighborhood types: Neumann & Moore
+   - Update modes: Synchronous & Asynchronous
+   - Rate parameters: randomly varied across valid ranges
+
+### 3. **Metrics Computed**
+   - **Stress Metric:** Normalized count of (predator/prey)↔empty adjacent pairs
+     - Represents friction and interface gradient (potential energy)
+   - **Avalanche Detection:** Population change magnitude thresholds
+   - **Population Variance:** Rolling window variance of prey/predator counts
+   - **Robustness Metrics:** 
+     - Avalanche count mean/std across configurations
+     - Magnitude consistency
+     - Coefficient of variation (measures criticality robustness)
+
+### 4. **Perturbation Experiment Design**
+Each experiment runs 230 total steps:
+- **Equilibration phase (0–80 steps):** System reaches quasi-steady state
+  - Stress accumulates during slow drive
+  - No parameter perturbation
+- **Observation phase (80–230 steps):** Gradual parameter drift
+  - Predator death rate increases by +0.05 (slow drive)
+  - System responds with cascading events if critical
+  - Stress release and avalanche events detected
+
+## Visualization Output
+
+The script generates `soc_analysis_results.png` with a **2×2 grid displaying the 4 core SOC properties**:
+
+1. **Panel 1 (Top-Left) - Slow Drive:** Gradual parameter drift over time with equilibration and perturbation phases marked
+2. **Panel 2 (Top-Right) - Build-up of Stress:** Stress accumulation with avalanche event thresholds marked as orange stars
+3. **Panel 3 (Bottom-Left) - Intermittent Release:** Prey and predator population dynamics showing cascade events during perturbation
+4. **Panel 4 (Bottom-Right) - Self-Organization:** Stress-density relation across diverse configurations, colored by avalanche activity
+
+## Usage
+
+```bash
+python scripts/soc_analysis.py
+```
+
+Output:
+- Console report with findings
+- PNG visualization saved to workspace root: `soc_analysis_results.png`
+
+## Key Observations from Test Run
+
+- **8 diverse configurations** sampled with varied grid sizes, densities, neighborhoods
+- **Avalanche detection:** 1/8 experiments showed clear cascade events
+- **Stress persistence:** Mean stress ~0.1529 across all configurations
+- **Robustness metric:** Coefficient of Variation = 2.646 (indicates some parameter-dependence; lower values →  more robust SOC)
+- **Population variance:** Consistent across runs (signature of intermittent release mechanism)
+
+## Code Structure
+
+### Main Functions
+- `compute_grid_stress()` – Interface-based stress metric
+- `compute_population_variance()` – Rolling window variance calculation
+- `detect_avalanche_events()` – Identify cascading population changes
+- `sample_parameter_configurations()` – Generate diverse parameter sets
+- `run_soc_perturbation_experiment()` – Single experiment with slow drive
+- `analyze_soc_robustness()` – Cross-configuration robustness metrics
+- `visualize_soc_properties()` – Comprehensive 8-panel figure
+- `main()` – Orchestrates full analysis pipeline
+
+### Configuration Space
+- **Grid size:** Affects stability and relaxation dynamics
+- **Densities:** Controls predator-prey interaction frequency
+- **Neighborhood:** Changes spatial coupling strength
+- **Rates:** Direct influence on birth/death thresholds
+
+## Scientific Interpretation
+
+The analysis tests the hypothesis:
+> *"Does the prey-predator CA exhibit self-organized criticality independent of specific parameter choices?"*
+
+If coefficient of variation is **low** (< 1.0) → SOC is **robust** (self-organized)
+If coefficient of variation is **high** (> 1.0) → Behavior is **parameter-dependent** (requires tuning)
+
+---
+
+**Created:** January 2026
+**Framework:** NumPy, Matplotlib, custom CA simulation

docs/sary_prompts.md

@@ -0,0 +1,135 @@
+# CA Stochastic Bifurcation Diagram:
+Mutations (evolutions) parameter OFF
+Control parameter: prey death rate
+
+Possible statistical observables:
+- Fraction of prey cells at equilibrium
+- Measure of entropy of the generated pattern.
+- Prey population count
+- Predator population count
+
+Run simulation:
+- Let the system run until a steady state is observed
+- For each death rate value, let the CA run for a specified number of iterations after warmp up, show distribution (scatters) for each sim run at a given prey death rate, and the average line
+
+
+# Phase 1: finding the critical point
+- Create bifurcation diagram of mean population count, varying prey death rate
+	- Look for critical transition
+- Create log-log plot of cluster size distribution, varying prey death rate
+	- Look for power-law
+
+# Experiment Phase: CA Stochastic Bifurcation Diagram:
+
+1) Write a Config Object specific to that experiment
+2) Make sure the experiment running on the cluster is running 15 reps of each runs at all sweeped values.
+3) Make sure the outputs of the experiment are a 1D and 2D array (explained below)
+
+# Bifurcation Diagram Prompts:
+1) Help me write a function for creating a stochastic bifurcation diagram, of the population count at equilibrium, varying the prey death rate (as the control variable). 
+2) At each sweeped value of the prey death control variable, we should be measuring the population count at equilibrium for at least 15 simulation runs. 
+3) Which means that the two inputs for my function should be a 1D Array for the sweep parameter, and a 2D array for the experiment results at each sweep for the rows, and the results for each iteration for the columns.
+4) When running my function, using the argparse module, my command-line arguments specifies which analysis to do, in this case the analysis is the bifurcation diagram.
+
+
+# Output: 
+def load_bifurcation_results(results_dir: Path) -> Tuple[np.ndarray, np.ndarray]:
+    """
+    Load bifurcation analysis results.
+    
+    Returns
+    -------
+    sweep_params : np.ndarray
+        1D array of control parameter values (prey death rates).
+    results : np.ndarray
+        2D array of shape (n_sweep, n_replicates) with population counts
+        at equilibrium.
+    """
+    npz_file = results_dir / "bifurcation_results.npz"
+    json_file = results_dir / "bifurcation_results.json"
+    
+    if npz_file.exists():
+        logging.info(f"Loading bifurcation results from {npz_file}")
+        data = np.load(npz_file)
+        return data['sweep_params'], data['results']
+    elif json_file.exists():
+        logging.info(f"Loading bifurcation results from {json_file}")
+        with open(json_file, 'r') as f:
+            data = json.load(f)
+        return np.array(data['sweep_params']), np.array(data['results'])
+    else:
+        raise FileNotFoundError(f"Bifurcation results not found in {results_dir}")
+
+
+def plot_bifurcation_diagram(sweep_params: np.ndarray, results: np.ndarray,
+                             output_dir: Path, dpi: int = 150,
+                             control_label: str = "Prey Death Rate",
+                             population_label: str = "Population at Equilibrium"):
+    """
+    Generate a stochastic bifurcation diagram.
+    
+    Shows the distribution of equilibrium population counts as a function of
+    a control parameter (e.g., prey death rate), with scatter points for each
+    replicate run overlaid on summary statistics.
+    
+    Parameters
+    ----------
+    sweep_params : np.ndarray
+        1D array of control parameter values (e.g., prey death rates).
+        Shape: (n_sweep,)
+    results : np.ndarray
+        2D array of population counts at equilibrium.
+        Shape: (n_sweep, n_replicates) where rows correspond to sweep_params
+        and columns are replicate simulation runs.
+    output_dir : Path
+        Directory to save the output figure.
+    dpi : int
+        Output resolution (default: 150).
+    control_label : str
+        Label for x-axis (control parameter).
+    population_label : str
+        Label for y-axis (population count).
+    """
+    n_sweep, n_replicates = results.shape
+    
+    fig, ax = plt.subplots(figsize=(12, 7))
+    
+    # Scatter all individual replicates with transparency
+    for i, param in enumerate(sweep_params):
+        ax.scatter(
+            np.full(n_replicates, param),
+            results[i, :],
+            alpha=0.3, s=15, c='steelblue', edgecolors='none'
+        )
+    
+    # Compute summary statistics
+    means = np.mean(results, axis=1)
+    medians = np.median(results, axis=1)
+    q25 = np.percentile(results, 25, axis=1)
+    q75 = np.percentile(results, 75, axis=1)
+    
+    # Plot median line and IQR envelope
+    ax.fill_between(sweep_params, q25, q75, alpha=0.25, color='coral',
+                    label='IQR (25th-75th percentile)')
+    ax.plot(sweep_params, medians, 'o-', color='darkred', linewidth=2,
+            markersize=5, label='Median')
+    ax.plot(sweep_params, means, 's--', color='black', linewidth=1.5,
+            markersize=4, alpha=0.7, label='Mean')
+    
+    ax.set_xlabel(control_label)
+    ax.set_ylabel(population_label)
+    ax.set_title(f"Stochastic Bifurcation Diagram\n({n_replicates} replicates per parameter value)")
+    ax.legend(loc='best')
+    ax.grid(True, alpha=0.3)
+    
+    # Add rug plot at bottom showing parameter sampling density
+    ax.plot(sweep_params, np.zeros_like(sweep_params), '|', color='gray',
+            markersize=10, alpha=0.5)
+    
+    plt.tight_layout()
+    output_file = output_dir / "bifurcation_diagram.png"
+    plt.savefig(output_file, dpi=dpi)
+    plt.close()
+    logging.info(f"Saved {output_file}")
+    
+    return output_file

models/CA.py

@@ -8,7 +8,14 @@
 
 import numpy as np
 import logging
+import sys
+from pathlib import Path
+
+# Add parent directory to path for imports
+sys.path.insert(0, str(Path(__file__).parent.parent))
+
 from models.numba_optimized import PPKernel, set_numba_seed
+from models.cluster_analysis import ClusterAnalyzer
 
 # Module logger
 logger = logging.getLogger(__name__)