diff --git a/notebooks/03_analysis.ipynb b/notebooks/03_analysis.ipynb index a03ce92..2ee26e3 100644 --- a/notebooks/03_analysis.ipynb +++ b/notebooks/03_analysis.ipynb @@ -10,7 +10,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 1, "id": "237bf30d", "metadata": {}, "outputs": [], @@ -24,7 +24,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 2, "id": "14500544", "metadata": {}, "outputs": [], @@ -43,12 +43,12 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 5, "id": "db6528cc", "metadata": {}, "outputs": [], "source": [ - "size = 100 # width and height of the grid\n", + "size = 500 # width and height of the grid\n", "p = 0.5 # starting fraction of vegetation\n", "update_rule = CA.update_Scanlon2007 # function containing update rule\n", "true_frac=0.1 # 'natural' (equilibrium) fraction of vegetation\n", @@ -68,7 +68,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 4, "id": "285dcba6", "metadata": {}, "outputs": [], @@ -85,6 +85,25 @@ ")" ] }, + { + "cell_type": "code", + "execution_count": 5, + "id": "764891c9", + "metadata": {}, + "outputs": [], + "source": [ + "grids = CA.evolve_CA_fast(\n", + " size=size,\n", + " p=p,\n", + " update_rule=CA.update_Scanlon2007_fast,\n", + " true_frac=true_frac,\n", + " k=k,\n", + " M=M,\n", + " N_steps=N_steps,\n", + " skip=skip,\n", + ")" + ] + }, { "cell_type": "markdown", "id": "f831e03a", @@ -155,7 +174,7 @@ ], "metadata": { "kernelspec": { - "display_name": "CLS", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -169,7 +188,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.7" + "version": "3.11.4" } }, "nbformat": 4, diff --git a/src/CA_model.py b/src/CA_model.py index 2f7e221..b925053 100644 --- a/src/CA_model.py +++ b/src/CA_model.py @@ -12,6 +12,7 @@ # Import necessary modules import numpy as np +from scipy.ndimage import convolve import matplotlib.pyplot as plt from numba import jit @@ -50,6 +51,126 @@ def update_basic(grid, cells_to_update): return grid +@jit(nopython=False) +def compute_weight_matrix(M, k, dmin=1): + """ + Precompute the distance-based weights for the neighborhood. + Returns a (2M+1) x (2M+1) matrix with weights (dmin/dist)**k, + and zero at the center (a cell doesn't count itself). + """ + size = 2 * M + 1 + weights = np.zeros((size, size)) + + for dx in range(-M, M + 1): + for dy in range(-M, M + 1): + dist = np.sqrt(dx**2 + dy**2) + if 0 < dist <= M: + weights[dx + M, dy + M] = (dmin / dist) ** k + + return weights + +@jit(nopython=False) +def compute_local_density(grid, weights): + """ + Compute the local weighted tree density for every cell simultaneously. + + For each cell, this calculates: + rho = sum(weights * neighbor_states) / sum(weights) + + Using convolution, this becomes a single operation over the whole grid. + """ + # Numerator: weighted sum of occupied neighbors + weighted_neighbors = convolve( + grid.astype(float), + weights, + mode='wrap' # periodic boundary conditions + ) + + # Denominator: sum of all weights (same for every cell with periodic BC) + total_weight = np.sum(weights) + + # Local density at each cell + rho = weighted_neighbors / total_weight + + return rho + +@jit(nopython=False) +def update_Scanlon2007_fast(grid, cells_to_update, true_frac, weights): + """ + Optimized Scanlon update rule using precomputed weights and convolution. + """ + size = grid.shape[0] + + # Compute local density for ALL cells at once (the expensive part, now fast) + rho = compute_local_density(grid, weights) + + # Global vegetation fraction + frac_occ = np.mean(grid) + + # Avoid division by zero + if frac_occ == 0 or frac_occ == 1: + return grid + + # Update only the selected cells + for i, j in cells_to_update: + rho_local = rho[i, j] + random_nr = np.random.random() + + if grid[i, j] == 0: # currently empty + prob_flip = rho_local + (true_frac - frac_occ) / (1 - frac_occ) + prob_flip = np.clip(prob_flip, 0, 1) # clamp to valid range + if random_nr < prob_flip: + grid[i, j] = 1 + + else: # currently occupied + prob_flip = (1 - rho_local) + (frac_occ - true_frac) / frac_occ + prob_flip = np.clip(prob_flip, 0, 1) + if random_nr < prob_flip: + grid[i, j] = 0 + + return grid + + +def evolve_CA_fast( + size=500, + p=0.5, + update_rule=update_Scanlon2007_fast, + true_frac=0.1, + k=3.0, + M=25, + f_update=0.2, + N_steps=200, + skip=100, + seed=0, +): + """ + Optimized evolution function. + """ + np.random.seed(seed) + grid = initialize_CA(p, size) + grids = [] + + # Precompute weights once (not every iteration!) + weights = compute_weight_matrix(M, k) + + # Fixed: use size**2 for total number of cells + N_update = int(f_update * size**2) + + for n in range(N_steps): + # Randomly select cells to update + cells_to_update = np.column_stack([ + np.random.randint(0, size, N_update), + np.random.randint(0, size, N_update) + ]) + + grid = update_rule(grid, cells_to_update, true_frac, weights) + + if n >= skip: + grids.append(grid.copy()) + + return grids + + @jit(nopython=False) def update_Scanlon2007( grid, cells_to_update, true_frac=0.3, k=3.0, M=25, dmin=1, BC="periodic" diff --git a/src/__pycache__/CA_model.cpython-311.pyc b/src/__pycache__/CA_model.cpython-311.pyc new file mode 100644 index 0000000..ee35e26 Binary files /dev/null and b/src/__pycache__/CA_model.cpython-311.pyc differ diff --git a/src/__pycache__/__init__.cpython-311.pyc b/src/__pycache__/__init__.cpython-311.pyc new file mode 100644 index 0000000..01e0bf6 Binary files /dev/null and b/src/__pycache__/__init__.cpython-311.pyc differ diff --git a/src/__pycache__/analysis.cpython-311.pyc b/src/__pycache__/analysis.cpython-311.pyc new file mode 100644 index 0000000..0a60c6b Binary files /dev/null and b/src/__pycache__/analysis.cpython-311.pyc differ