Percolation in a Porous Medium

GRID_SIZE = 20  # side length of the square grid (cells)
N_TRIALS = 200  # independent trials per probability value
N_PROBS = 41  # number of probability values from 0.0 to 1.0 inclusive
RNG_SEED = 7493418

# Theoretical site-percolation threshold for a 2D square lattice with
# 4-neighbor (von Neumann) connectivity.  From Stauffer & Aharony,
# Introduction to Percolation Theory (2nd ed., 1994), Appendix B.
THRESHOLD_THEORY = 0.5927

def make_grid(p, rng):
    """Return a GRID_SIZE x GRID_SIZE boolean array.

    Each cell is independently True (open) with probability p.
    """
    return rng.random((GRID_SIZE, GRID_SIZE)) < p

def percolates(grid):
    """Return True if an open path connects any top-row cell to any bottom-row cell.

    Uses iterative depth-first search with 4-neighbor (von Neumann) connectivity.
    All open cells in row 0 are added to the stack as starting points.
    The search succeeds as soon as any visited cell lies in the last row.
    """
    nrows, ncols = grid.shape
    visited = np.zeros_like(grid, dtype=bool)

    stack = [(0, c) for c in range(ncols) if grid[0, c]]
    for r, c in stack:
        visited[r, c] = True

    while stack:
        r, c = stack.pop()
        if r == nrows - 1:
            return True
        for dr, dc in ((-1, 0), (1, 0), (0, -1), (0, 1)):
            nr, nc = r + dr, c + dc
            if (
                0 <= nr < nrows
                and 0 <= nc < ncols
                and grid[nr, nc]
                and not visited[nr, nc]
            ):
                visited[nr, nc] = True
                stack.append((nr, nc))
    return False

def sweep():
    """Return a DataFrame with the percolation fraction at each probability.

    For each p in a uniform grid from 0 to 1, runs N_TRIALS independent
    experiments and records the fraction that percolate.  All trials share
    a single RNG to ensure reproducibility from RNG_SEED.
    """
    rng = np.random.default_rng(RNG_SEED)
    probs = np.linspace(0.0, 1.0, N_PROBS)
    fractions = [
        sum(percolates(make_grid(p, rng)) for _ in range(N_TRIALS)) / N_TRIALS
        for p in probs
    ]
    return pl.DataFrame({"probability": probs, "fraction_percolating": fractions})

def estimate_threshold(df):
    """Return the smallest p where the percolation fraction first reaches 0.5.

    On a finite grid the percolation transition is sharp but not a true step
    function.  The crossing of 0.5 is a conventional and reproducible
    estimator of the finite-size threshold.
    """
    above = df.filter(pl.col("fraction_percolating") >= 0.5)
    return float(above["probability"].min()) if not above.is_empty() else float("nan")

def plot(df):
    """Return an Altair line chart of percolation fraction vs. open probability."""
    return (
        alt.Chart(df)
        .mark_line(point=True)
        .encode(
            x=alt.X("probability:Q", title="Open-cell probability p"),
            y=alt.Y(
                "fraction_percolating:Q",
                title="Fraction of trials that percolate",
                scale=alt.Scale(domain=[0.0, 1.0]),
            ),
        )
        .properties(width=400, height=300, title="Percolation threshold sweep")
    )

import numpy as np
from porous import (
    RNG_SEED,
    THRESHOLD_THEORY,
    estimate_threshold,
    make_grid,
    percolates,
    sweep,
)


def test_empty_grid_never_percolates():
    # When p = 0 every cell is closed (False), so the stack is empty from
    # the start and the DFS returns False immediately.
    rng = np.random.default_rng(RNG_SEED)
    grid = make_grid(0.0, rng)
    assert not percolates(grid)


def test_full_grid_always_percolates():
    # When p = 1 every cell is open (True).  A straight vertical path
    # through column 0 connects row 0 to row GRID_SIZE - 1.
    rng = np.random.default_rng(RNG_SEED)
    grid = make_grid(1.0, rng)
    assert percolates(grid)


def test_single_row_grid():
    # A one-row grid percolates if and only if at least one cell in that row
    # is open, because the top row and the bottom row are the same row.
    open_row = np.array([[True, False, True]])
    assert percolates(open_row)
    closed_row = np.array([[False, False, False]])
    assert not percolates(closed_row)


def test_fractions_in_unit_interval():
    # Each trial either percolates or not, so the fraction must lie in [0, 1].
    df = sweep()
    assert (df["fraction_percolating"] >= 0.0).all()
    assert (df["fraction_percolating"] <= 1.0).all()


def test_threshold_near_theory():
    # On a 20 x 20 grid with 200 trials the finite-size threshold should be
    # within 0.07 of the infinite-lattice value 0.5927.  Finite-size effects
    # and Monte Carlo variance both contribute; a tolerance of 0.07 gives a
    # comfortable margin without being so wide that a wrong algorithm passes.
    df = sweep()
    p_c = estimate_threshold(df)
    assert abs(p_c - THRESHOLD_THEORY) < 0.07, (
        f"estimated threshold {p_c:.3f} is too far from theory {THRESHOLD_THEORY:.4f}"
    )