Merge pull request #20 from codegithubka/storm

Directed reproduction and phase 6 config

Author: StormKoper

docs/experiments.md

@@ -32,7 +32,7 @@ These phases should be completed sequentially, deepening our understanding at ea
 	- Look for critical slowing down: perturbations to states closer to the critical point should more slowly return to the steady state
 	- This requires time series data
 ### Phase 6: model extensions
-- Investigate whether hydra effect and SOC still occur with diffusion and directed movement
+- Investigate whether hydra effect and SOC still occur with diffusion and directed reproduction
 
 # Todo
 The main functionality is all complete. Thus, the models folder should be relatively untouched.

models/config.py

@@ -263,20 +263,16 @@ def estimate_runtime(self, n_cores: int = 32) -> str:
     timeseries_subsample=1,  # Full resolution for autocorrelation
 )
 
-# Phase 6: Model extensions (directed hunting); same config as phase 1 but with directed hunting
+# Phase 6: Model extensions (directed reproduction); same config as phase 4 but with directed reproduction
 PHASE6_CONFIG = Config(
-    grid_size=1000,
-    n_prey_death=20,
-    prey_birth=0.2,
-    prey_death_range=(0.09, 0.12),
-    predator_birth=0.8,    
-    predator_death=0.05,    
-    n_replicates=30,
-    warmup_steps=1000,
-    measurement_steps=1000,
-    collect_pcf=True,
-    pcf_sample_rate=0.2,
-    save_timeseries=False,
+    grid_size=250,
+    n_replicates=10,
+    warmup_steps=500,
+    measurement_steps=500,
+    with_evolution=False,
+    collect_pcf=False,
+    save_timeseries=True,
+    timeseries_subsample=10,
     directed_hunting=True,
 )
 

models/numba_optimized.py

@@ -56,6 +56,65 @@ def set_numba_seed(seed: int) -> None:
 # ============================================================================
 # PREDATOR-PREY KERNELS
 # ============================================================================
+@njit(cache=True)
+def _pp_async_kernel_fast(
+    grid: np.ndarray,
+    prey_death_arr: np.ndarray,
+    p_birth_val: float,
+    p_death_val: float,
+    pred_birth_val: float,
+    pred_death_val: float,
+    dr_arr: np.ndarray,
+    dc_arr: np.ndarray,
+    evolve_sd: float,
+    evolve_min: float,
+    evolve_max: float,
+    evolution_stopped: bool,
+    occupied_buffer: np.ndarray,
+) -> np.ndarray:
+    """Partially synchronous predator-prey update kernel."""
+    rows, cols = grid.shape
+    n_shifts = len(dr_arr)
+    grid_copy = grid.copy()
+    prey_death_arr_copy = prey_death_arr.copy()
+
+    prey_death = np.random.random(size=grid.shape)
+    grid_copy[(grid == 1) & (prey_death < prey_death_arr)] = 0
+    prey_death_arr_copy[(grid == 1) & (prey_death < prey_death_arr)] = np.nan
+
+    pred_death = np.random.random(size=grid.shape)
+    grid_copy[(grid == 2) & (pred_death < pred_death_val)] = 0
+
+    count = np.count_nonzero(grid)
+    indices = np.random.permutation(count)
+    rs = indices // cols
+    cs = indices % cols
+
+    nb = np.random.randint(0, n_shifts, size=count)
+    nrs = (rs + dr_arr[nb]) % rows
+    ncs = (cs + dc_arr[nb]) % cols
+
+    for r, c, nr, nc in zip(rs, cs, nrs, ncs):
+        state = grid[r, c]
+        nstate = grid[nr, nc]
+
+        if state == 1 and nstate == 0 and np.random.random() < p_birth_val:
+            grid_copy[nr, nc] = 1
+            parent_val = prey_death_arr[r, c]
+            if not evolution_stopped:
+                child_val = parent_val + np.random.normal(0, evolve_sd)
+                prey_death_arr_copy[nr, nc] = np.clip(child_val, evolve_min, evolve_max)
+            else:
+                prey_death_arr_copy[nr, nc] = parent_val
+
+        elif state == 2 and nstate == 1 and np.random.random() < pred_birth_val:
+            grid_copy[nr, nc] = 2
+            prey_death_arr_copy[nr, nc] = np.nan
+
+    grid = grid_copy
+    prey_death_arr = prey_death_arr_copy
+
+    return grid
 
 @njit(cache=True)
 def _pp_async_kernel_random(
@@ -148,7 +207,7 @@ def _pp_async_kernel_directed(
     evolution_stopped: bool,
     occupied_buffer: np.ndarray,
 ) -> np.ndarray:
-    """Async predator-prey update kernel with directed hunting."""
+    """Async predator-prey update kernel with directed reproduction."""
     rows, cols = grid.shape
     n_shifts = len(dr_arr)
 
@@ -170,19 +229,22 @@ def _pp_async_kernel_directed(
         c = occupied_buffer[i, 1]
 
         state = grid[r, c]
-        if state == 0:
-            continue
 
         if state == 1:  # PREY
-            nbi = np.random.randint(0, n_shifts)
-            nr = (r + dr_arr[nbi]) % rows
-            nc = (c + dc_arr[nbi]) % cols
-            
             if np.random.random() < prey_death_arr[r, c]:
                 grid[r, c] = 0
                 prey_death_arr[r, c] = np.nan
-            elif grid[nr, nc] == 0:
-                if np.random.random() < p_birth_val:
+            elif np.random.random() < p_birth_val:
+                valid = []
+                for k in range(n_shifts):
+                    nr = (r + dr_arr[k]) % rows
+                    nc = (c + dc_arr[k]) % cols
+                    if grid[nr, nc] == 0:
+                        valid.append(k)
+                if len(valid) > 0:
+                    choice = np.random.choice(valid)
+                    nr = (r + dr_arr[choice]) % rows
+                    nc = (c + dc_arr[choice]) % cols
                     grid[nr, nc] = 1
                     parent_val = prey_death_arr[r, c]
                     if not evolution_stopped:
@@ -195,44 +257,22 @@ def _pp_async_kernel_directed(
                     else:
                         prey_death_arr[nr, nc] = parent_val
 
-        elif state == 2:  # PREDATOR - directed hunting
+        elif state == 2:  # PREDATOR
             if np.random.random() < pred_death_val:
                 grid[r, c] = 0
-                continue
-            
-            prey_count = 0
-            for k in range(n_shifts):
-                check_r = (r + dr_arr[k]) % rows
-                check_c = (c + dc_arr[k]) % cols
-                if grid[check_r, check_c] == 1:
-                    prey_count += 1
-            
-            if prey_count > 0:
-                target_idx = np.random.randint(0, prey_count)
-                found = 0
-                nr, nc = 0, 0
+            elif np.random.random() < pred_birth_val:
+                valid = []
                 for k in range(n_shifts):
-                    check_r = (r + dr_arr[k]) % rows
-                    check_c = (c + dc_arr[k]) % cols
-                    if grid[check_r, check_c] == 1:
-                        if found == target_idx:
-                            nr = check_r
-                            nc = check_c
-                            break
-                        found += 1
-                
-                if np.random.random() < pred_birth_val:
+                    nr = (r + dr_arr[k]) % rows
+                    nc = (c + dc_arr[k]) % cols
+                    if grid[nr, nc] == 1:
+                        valid.append(k)
+                if len(valid) > 0:
+                    choice = np.random.choice(valid)
+                    nr = (r + dr_arr[choice]) % rows
+                    nc = (c + dc_arr[choice]) % cols
                     grid[nr, nc] = 2
                     prey_death_arr[nr, nc] = np.nan
-            else:
-                nbi = np.random.randint(0, n_shifts)
-                nr = (r + dr_arr[nbi]) % rows
-                nc = (c + dc_arr[nbi]) % cols
-                
-                if grid[nr, nc] == 1:
-                    if np.random.random() < pred_birth_val:
-                        grid[nr, nc] = 2
-                        prey_death_arr[nr, nc] = np.nan
 
     return grid