Merge branch 'main' into storm

Author: StormKoper

docs/HPC_GUIDE.md

@@ -0,0 +1,50 @@
+### Snellius Usage Breakdown
+
+```
+ssh <your_username>@snellius.surf.nl
+
+# On a separate terminal run the following
+
+# Upload the entire project directory (including your models/ folder)
+scp -r ~/Documents/CSS_Project <your_username>@snellius.surf.nl:~/
+
+# On the Snellius terminal
+
+module load 2023 Python/3.11.3-GCCcore-12.3.0
+python3 -m venv ~/css_env
+source ~/css_env/bin/activate
+pip install numpy scipy matplotlib joblib
+
+# To do a dry run for testing the entire environment
+
+python3 pp_analysis.py --mode full --dry-run
+
+# For async run
+
+python3 pp_analysis.py --mode full --output results_${SLURM_JOB_ID} --cores $SLURM_CPUS_PER_TASK --async
+
+# To submit a job
+
+sbatch run_analysis.sh
+
+# Check Queue Status
+
+squeue -u $USER
+
+# Cancel a job
+
+scancel <JOBID>
+
+# Monitoring live progress
+
+tail -f logs_<JOBID>.err
+
+# Fetching the results once the job is done
+
+scp -r <your_username>@snellius.surf.nl:~/results_18514601 ~/Downloads/
+```
+
+The jobscript template can be found in ```run_analysis.sh``` (default rome paritition).
+
+
+Snellius Partitions Page: https://servicedesk.surf.nl/wiki/spaces/WIKI/pages/30660209/Snellius+partitions+and+accounting

docs/prompts.md

@@ -691,4 +691,20 @@ Lastly, add functionality to stop evolution after a certain time-step. This shou
 
 7. Create a bifuracation diagram to confirm the monotonic relationship for a varying prey death rate vs. equilibrium density. 
 
-###
+---
+
+### Testing CA class
+
+1. Create a comprehensive testing suite for the CA and PP classes. Test initialization, async update changes, synchronous update changes, prey growth in isolation behavior, predator starvation, parameter evolution and long run dynamics. Also make sur ethe test_viz mehtod works as desired
+
+### Parameter Sweep and PP Class Analysis
+
+2. Create a skeletal version of a .py script that will be subimtted into Snellius for parameter analysis. The purpose of this script should be to identify power law distribution with the consideration of finite size scaling, the hydra effect, and suitable parameter configurtaions for time series analysis for model evolution. Compare the non-evo to the evo model.
+
+3. Create a config class adjustable depending on the CPU budget. We want to run a prey_birth vs. predator_death parameter sweep (2D), quantify the hydra effect using the derivative, search for the critical point (power law relartionship paramameter), quantify evolution sensitivity and analyze finite grid size scaling. Include script options for cost optimal runs as well. Make sure to have a summary of collected data stored for reference and future usage.
+
+
+4. Add configuration option to run the asynchronous version of the CA class. The synchronous functionality should also be preserved. Provide me with a small smoke test to see if the updated file runs as expected.
+
+5. Create a minimal bash script for Snellius. Use the rome configiuration.
+

models/CA.py

@@ -208,6 +208,9 @@ def run(self, steps: int, stop_evolution_at: Optional[int] = None, snapshot_iter
 		# normalize snapshot iteration list
 		snapshot_set = set(snapshot_iters) if snapshot_iters is not None else set()
 
+		# normalize snapshot iteration list
+		snapshot_set = set(snapshot_iters) if snapshot_iters is not None else set()
+
 		for i in range(steps):
 			self.update()
 			# Update visualization if enabled every `interval` iterations
@@ -241,6 +244,30 @@ def run(self, steps: int, stop_evolution_at: Optional[int] = None, snapshot_iter
 				# disable further evolution
 				self._evolve_info = {}
 
+			# create snapshots if requested at this iteration
+			if (i + 1) in snapshot_set:
+				try:
+					# create snapshot folder if not present
+					if not hasattr(self, "_viz_snapshot_dir") or self._viz_snapshot_dir is None:
+						import os, time
+
+						base = "results"
+						ts = int(time.time())
+						run_folder = f"run-{ts}"
+						full = os.path.join(base, run_folder)
+						os.makedirs(full, exist_ok=True)
+						self._viz_snapshot_dir = full
+					self._viz_save_snapshot(i + 1)
+				except Exception:
+					pass
+
+			# stop evolution at specified time-step (disable further evolution)
+			if stop_evolution_at is not None and (i + 1) == int(stop_evolution_at):
+				try:
+					self._evolve_info = {}
+				except Exception:
+					pass
+
 	def visualize(
 		self,
 		interval: int = 1,

models/__init__.py

@@ -1,2 +1,2 @@
 from .mean_field import MeanFieldModel
-from .CA import CA
+from .CA import CA, PP

models/firebreak.py

@@ -0,0 +1,64 @@
+""""
+Altruistic Suicide Firebreak Model
+
+Concept:
+
+- Standard Prey: Dies of old age or being eaten
+- Firebreak prey: Has a sensor to detect if a predator is nearby. If so, then it has
+    the option to protect the colony by commiting suicde before the predator can reproduce
+    into it
+- Result: Creates empty cells around the predator cluster. Prevents percolation through
+the cluster
+
+
+TO DO:
+
+1. Override update_sync and update_async logic.
+    New Logic:
+        1. Calculate predator_neighbors for every cell
+        2. Create a death rate grid where every cell gets assigned a death rate based on neighbors
+        3. Calc the mask against the dynamic grid
+"""
+
+import numpy as np
+from models.CA import PP
+
+
+class Firebreak(PP):
+    """
+    PP CA where prey commit suicide when threatend to creaty empty firebreaks
+    that starve predators
+    """
+
+    def __init__(
+        self,
+        rows,
+        cols,
+        densities,
+        neighborhood="moore",
+        params=None,
+        cell_params=None,
+        seed=None,
+        synchronous=True,
+    ):
+        # get firebreak specific parameter
+        self.alt_dr = 0.5  # FIXME: Random default value
+        clean_params = dict(params) if params else {}
+        if "altruistic_dr" in clean_params:
+            self.alt_dr = float(clean_params.pop("altruistic_dr"))
+
+        super().__init__(
+            rows,
+            cols,
+            densities,
+            neighborhood,
+            clean_params,
+            cell_params,
+            seed,
+            synchronous,
+        )
+        self.params["altruistic_dr"] = self.alt_dr
+
+    def update_sync(self) -> None:
+        """Override syncrhonous update. Similar to PP update except Prey Death Logic."""
+        pass

models/helpers.py

@@ -1,78 +1,73 @@
 import numpy as np
 from scipy.ndimage import label
 
-def find_clusters(grid: np.ndarray, state: int = 1, connectivity: int = 'moore'):
+
+def find_clusters(grid: np.ndarray, state: int = 1, connectivity: int = "moore"):
     """Locate connected clusters of a given state."""
-    if connectivity == 'moore':
+    if connectivity == "moore":
         structure = np.ones((3, 3), dtype=int)
     else:
-        structure =  ([[0, 1, 0],
-                       [1, 1, 1],
-                      [0, 1, 0]])
-        
+        structure = [[0, 1, 0], [1, 1, 1], [0, 1, 0]]
+
     binary_mask = (grid == state).astype(int)
     labeled_array, num_clusters = label(binary_mask, structure=structure)
-    
+
     return labeled_array, num_clusters
 
+
 def get_cluster_sizes(grid: np.ndarray, state: int = 1):
     """Returns an array containing the sizes of all detected clusters."""
     labeled, n = find_clusters(grid, state)
-    
+
     if n == 0:
         return np.array([])
-    
+
     sizes = np.bincount(labeled.ravel())[1:]  # Exclude background count
-    
+
     return sizes
 
+
 def fit_power_law_mle(sizes: np.ndarray, s_min: float = 1.0):
     """
     Maximum Likelihood Estimator for the power-law exponent.
     """
     # Filter sizes below the minimum threshold
     sizes = sizes[sizes >= s_min]
-    
+
     if len(sizes) < 10:
-        return {'tau': np.nan, 'tau_err': np.nan, 'n_samples': len(sizes)}
-    
+        return {"tau": np.nan, "tau_err": np.nan, "n_samples": len(sizes)}
+
     n = len(sizes)
-    
+
     # Standard MLE formula for power-law index
     tau = 1 + n / np.sum(np.log(sizes / s_min))
     tau_err = (tau - 1) / np.sqrt(n)
-    
-    return {
-        'tau': tau, 
-        'tau_err': tau_err, 
-        'n_samples': n
-    }
+
+    return {"tau": tau, "tau_err": tau_err, "n_samples": n}
+
 
 def detect_hydra_effect(d_r_vals: np.ndarray, R_vals: np.ndarray):
     """Detects hydra effect by identifying increase in prey population with increased death rate."""
-    
+
     dR = np.gradient(R_vals, d_r_vals)
     hydra_mask = dR > 0
-    
+
     if not np.any(hydra_mask):
-        return {
-            'detected': False,
-            'region': None,
-            'max_strength': 0.0
-        }
+        return {"detected": False, "region": None, "max_strength": 0.0}
     hydra_indices = np.where(hydra_mask)[0]
-    
+
     results = {
-        'detected': True,
-        'd_r_start': d_r_vals[hydra_indices[0]],
-        'd_r_end': d_r_vals[hydra_indices[-1]],
-        'max_slope': np.max(dR),
-        'd_r_at_max_slope': d_r_vals[np.argmax(dR)],
-        'peak_density': np.max(R_vals[hydra_mask])
+        "detected": True,
+        "d_r_start": d_r_vals[hydra_indices[0]],
+        "d_r_end": d_r_vals[hydra_indices[-1]],
+        "max_slope": np.max(dR),
+        "d_r_at_max_slope": d_r_vals[np.argmax(dR)],
+        "peak_density": np.max(R_vals[hydra_mask]),
     }
-    
+
     return results
 
+
 def find_spatial_correlation(d_r_vals, ca_results):
     """ca_results: List of grids (np.ndarray) at equilibrium of each d_r value."""
-    pass
+    pass