Merge branch 'kimon' into storm

Author: codegithubka

.gitignore

@@ -1,2 +1,12 @@
 __pycache__/
-results/
+results/
+# Virtual environment
+.venv/
+__pycache__/
+# Python bytecode
+*.pyc
+# Jupyter Notebook checkpoints
+.ipynb_checkpoints/
+# Data files
+data/
+.vscode/

docs/prompts.md

@@ -0,0 +1,32 @@
+### Create base CA class
+
+Create a cellular automaton class named CA with an init function, a count neighbors function, an update function, and a run. The CA should consist of a numpy array called grid, a string containing the neighborhood type, and a numpy random number generator called generator. Use this generator for all random number generation inside the class. The CA class should also contain a dictionary called params.
+
+The init function should take arguments for the grid size (rows, columns, both ints), the initial density of each species (a tuple of floats), the neighborhood type ("neumann" or "moore"), the parameters in the form of a dictionary, and the seed for the generator. It should assign the parameters dictionary to the params variable, create the generator object and assign it to the generator variable, as well as create the 2D array of zeros based on the grid size and assign it to the grid variable. This grid should then be filled with states dependent on the density tuple. Iterate over the elements i of this tuple, filling grid_size * density[i] elements of the grid with state i+1. Non-zero cell states should not be overwritten, ensuring that the specified percentage of the grid is filled with that state. It should also check if the neighborhood argument corresponds with a known neighborhood and return an error otherwise.
+
+The count neighbors function should return a tuple of matrices (one for each defined non-zero state) containing the amount of neighbors of that state for each cell. It should use the neighborhood defined in the class. Ensure the logic works for both "neumann" and "moore". Use periodic boundaries.
+
+The update function can remain empty, so fill it with a pass statement.
+
+The run function should take a steps (int) argument. It should then run the CA for steps interations, calling the update function each time.
+
+Finally, make sure to add an expected type for each argument and define the return types. Add this information, as well as a short description of the function to the docstring. Also add assert statements to ensure arguments "make sense". For example, the sum of densities should not exceed 1 and the rows, cols, densities should all be positive, and the neighborhood should be either "neumann" or "moore".
+
+
+### Mean Field class
+
+1. Create a baseline mean-field class based on the attached research paper on predator-prey dynamics. The class should adhere to the papers specifications. The class should have a parameter sweep method for key predator and prey parameters that will be run in Snellius. Also include a method for equilibrium analysis. Make sure to justify the logic for this method. Include docstrings with a small method description and comments for code interpretability.
+
+2. Justify initialization parameter values for a small test expiriment. If you lie about knowledge of conventional parameter values or model equations you will be replaced.
+
+3. Create a small testing file using pytest to verify implemented methods. Make sure to cover edge cases and list them after the .py file output for me please. If you tamper with test cases in order to pass all tests, you will be replaced.
+
+4. We are now ready to plot some of the results of the mean fielf baseline. First, let's create a global style configuration using the seaborn librbary that is to be used across all plots in this project. Make sure the legend is at the bottom of each plot.
+
+5. Plot the phase portait to confirm the system spiral into a stable point. Show the nullclines as well. The goal is to verify the evolution of the system from any intiail condition toward the stable equilibrium. 
+
+6. Create a time series analysis plot of the evolution of prey and predator density vs. time. Make sure enough time steps all visible to see how the system eventually stabilizes.
+
+7. Create a bifuracation diagram to confirm the monotonic relationship for a varying prey death rate vs. equilibrium density. 
+
+###

models/CA.py

@@ -0,0 +1,126 @@
+from typing import Tuple, Dict, Any, Optional
+import numpy as np
+
+class CA:
+	"""Cellular Automaton container.
+
+	Attributes
+	- `grid` (np.ndarray): 2D integer array with cell states (0 = empty).
+	- `neighborhood` (str): either "neumann" or "moore".
+	- `generator` (np.random.Generator): RNG used for all random operations.
+	- `params` (dict): user-provided parameters.
+	"""
+
+	def __init__(
+		self,
+		rows: int,
+		cols: int,
+		densities: Tuple[float, ...],
+		neighborhood: str = "neumann",
+		params: Optional[Dict[str, Any]] = None,
+		seed: Optional[int] = None,
+	) -> None:
+		"""Initialize the CA.
+
+		Args:
+		- `rows` (int): number of rows (> 0).
+		- `cols` (int): number of columns (> 0).
+		- `densities` (tuple of float): density fraction for each species (each >= 0).
+		- `neighborhood` (str): either "neumann" or "moore".
+		- `params` (dict, optional): additional parameters to store.
+		- `seed` (int or None): seed for the RNG.
+
+		Returns:
+		- None
+
+		Raises:
+		- AssertionError if arguments are invalid.
+		"""
+		assert isinstance(rows, int) and rows > 0, "`rows` must be a positive int"
+		assert isinstance(cols, int) and cols > 0, "`cols` must be a positive int"
+		assert isinstance(densities, tuple) or isinstance(densities, list), "`densities` must be a tuple or list"
+		densities = tuple(float(d) for d in densities)
+		assert all(d >= 0 for d in densities), "all densities must be non-negative"
+		assert sum(densities) <= 1.0 + 1e-12, "sum of densities must not exceed 1"
+		assert neighborhood in ("neumann", "moore"), "neighborhood must be 'neumann' or 'moore'"
+
+		self.params: Dict[str, Any] = dict(params) if params is not None else {}
+		self.generator: np.random.Generator = np.random.default_rng(seed)
+		self.neighborhood: str = neighborhood
+
+		self.rows = rows
+		self.cols = cols
+		self.grid: np.ndarray = np.zeros((rows, cols), dtype=int)
+
+		total_cells = rows * cols
+
+		# Fill grid for each species (states 1..N) according to densities.
+		# Use floor to avoid over-allocation; do not overwrite non-zero cells.
+		for i, d in enumerate(densities):
+			if d <= 0:
+				continue
+			desired = int(np.floor(total_cells * d))
+			if desired <= 0:
+				continue
+
+			# available positions (flat indices)
+			available = np.flatnonzero(self.grid.ravel() == 0)
+			if len(available) == 0:
+				break
+
+			take = min(desired, len(available))
+			chosen = self.generator.choice(available, size=take, replace=False)
+			self.grid.ravel()[chosen] = i + 1
+
+		# store number of species expected based on densities length
+		self._n_species = len(densities)
+
+	def count_neighbors(self) -> Tuple[np.ndarray, ...]:
+		"""Count neighbors for each non-zero state.
+
+		Uses periodic boundary conditions and the neighborhood specified in the instance.
+
+		Returns:
+		- Tuple of np.ndarray: each array has shape `(rows, cols)` and contains
+		  the count of neighbors of state k (for k = 1..N) at each cell.
+		"""
+		neighbors = []
+
+		if self.neighborhood == "moore":
+			offsets = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]
+		else:  # neumann
+			offsets = [(-1, 0), (1, 0), (0, -1), (0, 1)]
+
+		for state in range(1, self._n_species + 1):
+			mask = (self.grid == state).astype(int)
+			count = np.zeros_like(self.grid, dtype=int)
+			for dr, dc in offsets:
+				shifted = np.roll(np.roll(mask, dr, axis=0), dc, axis=1)
+				count += shifted
+			neighbors.append(count)
+
+		return tuple(neighbors)
+
+	def update(self) -> None:
+		"""Perform a single update of the CA.
+
+		This method is intentionally left as a stub; implement specific rule
+		logic by overriding or editing this method.
+
+		Returns:
+		- None
+		"""
+		pass
+
+	def run(self, steps: int) -> None:
+		"""Run the CA for a number of iterations.
+
+		Args:
+		- `steps` (int): number of iterations to execute (>= 0).
+
+		Returns:
+		- None
+		"""
+		assert isinstance(steps, int) and steps >= 0, "`steps` must be a non-negative int"
+		for _ in range(steps):
+			self.update()

models/__init__.py

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

models/helpers.py

@@ -0,0 +1,78 @@
+import numpy as np
+from scipy.ndimage import label
+
+def find_clusters(grid: np.ndarray, state: int = 1, connectivity: int = 'moore'):
+    """Locate connected clusters of a given state."""
+    if connectivity == 'moore':
+        structure = np.ones((3, 3), dtype=int)
+    else:
+        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)}
+    
+    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
+    }
+
+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
+        }
+    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])
+    }
+    
+    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