Replicating Sandpile Research: Nature Fig 1a) From Greedy Control Of Cascading Failures In Interdependent Networks

Mu Controlled Sandpile Model 

#I switched to python from C to make the sandpile model to make my life easier

#replicating Nature fig 1a: https://www.nature.com/articles/s41598-021-82843-8/figures/1

import numpy as np

import matplotlib.pyplot as plt

import random

def make_area(x, y):

    return [[0]*x for _ in range(y)]

def find_at_capacity_sites(height, XMAX, YMAX):

    return [(x,y) for x in range(XMAX) for y in range(YMAX) if height[y][x] >= 3] 

def find_not_at_capacity_sites(height, XMAX, YMAX):

    return [(x,y) for x in range(XMAX) for y in range(YMAX) if height[y][x] < 3] 

def flip(mu):

    return random.random() < mu

def add_grain(area, mu, XMAX, YMAX):

    at_capacity = find_at_capacity_sites(area, XMAX, YMAX)

    not_at_capacity = find_not_at_capacity_sites(area, XMAX, YMAX)

    if not at_capacity:

        rx, ry = random.choice(not_at_capacity)

    else:

        if flip(mu):

            rx, ry = random.choice(at_capacity)

        else:

            rx, ry = random.choice(not_at_capacity)

    area[ry][rx] += 1

    return (rx, ry)

def topple(area, XMAX, YMAX): 

    avalanche_size = 0

    unstable_sites = []

    # Find initial unstable sites

    for y in range(YMAX):

        for x in range(XMAX):

            if area[y][x] >= 4:

                unstable_sites.append((x, y))

    

    while unstable_sites:

        x, y = unstable_sites.pop()

        if area[y][x] < 4:

            continue

        avalanche_size += 1

        area[y][x] -= 4

        # Update neighbors

        for dx, dy in [(1,0), (-1,0), (0,1), (0,-1)]:

            nx, ny = x + dx, y + dy

            if 0 <= nx < XMAX and 0 <= ny < YMAX:

                area[ny][nx] += 1

                if area[ny][nx] >= 4:  # Check if neighbor becomes unstable

                    unstable_sites.append((nx, ny))

    return avalanche_size

def run_sandpile(mu, grains_to_add, XMAX, YMAX):

    area = make_area(XMAX, YMAX)

    avalanche_sizes = []

    for _ in range(grains_to_add):

        add_grain(area, mu, XMAX, YMAX)

        aval_size = topple(area, XMAX, YMAX)

        if aval_size > 0:

            avalanche_sizes.append(aval_size)

    return avalanche_sizes

def plot_loglog_cascade_size(avalanche_data, mu_values):

    plt.figure(figsize=(8,6))

    

    # Use markers/colors from paper

    markers = {0.20: 'o', 0.37: 's', 0.90: '^'}  # circle, square, triangle

    colors = {0.20: 'blue', 0.37: 'green', 0.90: 'red'}

    

    for mu, avalanches in zip(mu_values, avalanche_data):

        if len(avalanches) == 0:

            continue

            

        max_s = max(avalanches)

        bins = np.logspace(0, np.log10(max_s), 30)

        hist, bin_edges = np.histogram(avalanches, bins=bins, density=True)

        bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

        

        # Plot only non-zero points

        non_zero = hist > 0

        plt.loglog(bin_centers[non_zero], hist[non_zero], 

                   marker=markers[mu], linestyle='-', lw=1.5, ms=6,

                   color=colors[mu], label=f'μ = {mu}')

    

    plt.xlabel("Cascade size, s", fontsize=14)

    plt.ylabel("Probability, P(s)", fontsize=14)

    plt.title("Avalanche Size Distribution", fontsize=16)

    plt.legend(fontsize=12)

    plt.grid(True, which="both", ls="--", alpha=0.7)

    plt.xlim(1, 1e5)  # Match paper's axis limits

    plt.ylim(1e-6, 1)

    plt.tight_layout()

    plt.show()

if __name__ == '__main__':

    XMAX, YMAX = 100, 100

    grains_to_add = 50000

    mu_values = [0.20, 0.37, 0.90]

    avalanche_data = []

    for mu in mu_values:

        print(f"Running sandpile for μ = {mu}...")

        avalanches = run_sandpile(mu, grains_to_add, XMAX, YMAX)

        avalanche_data.append(avalanches)

    

    plot_loglog_cascade_size(avalanche_data, mu_values)