# /// script
# requires-python = ">=3.12"
# dependencies = ["marimo>=0.25.1", "numpy", "matplotlib"]
# ///

import marimo

__generated_with = "0.25.1"
app = marimo.App(width="medium")


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        # Spots and stripes from two chemicals

        In 1952 Alan Turing proposed that two substances which react and spread
        at different rates can turn a uniform tissue into a pattern: the spots
        of a leopard, the stripes of a zebrafish. Forty years later John
        Pearson showed how rich such patterns can be with a two-line model, the
        Gray-Scott reaction:

        $$\frac{\partial u}{\partial t} = D_u \nabla^2 u - u v^2 + F(1 - u), \qquad
        \frac{\partial v}{\partial t} = D_v \nabla^2 v + u v^2 - (F + k)\, v$$

        $u$ is fed in at rate $F$, $v$ consumes it to make more of itself and is
        removed at rate $F + k$, and $u$ spreads twice as fast as $v$. Small
        changes in $F$ and $k$ give spots that divide, mazes, or coral.
        """
    )
    return


@app.cell
def _():
    import marimo as mo
    import numpy as np
    import matplotlib.pyplot as plt
    return mo, np, plt


@app.cell
def _(mo):
    PRESETS = {
        "Dividing spots": (0.0367, 0.0649),
        "Spots": (0.035, 0.065),
        "Coral": (0.0545, 0.062),
        "Maze": (0.029, 0.057),
        "Worms": (0.078, 0.061),
    }
    preset = mo.ui.dropdown(list(PRESETS), value="Dividing spots", label="Pattern")
    preset
    return PRESETS, preset


@app.cell
def _(PRESETS, mo, preset):
    feed = mo.ui.slider(0.01, 0.09, step=0.0005, value=PRESETS[preset.value][0], label="Feed rate F")
    kill = mo.ui.slider(0.045, 0.07, step=0.0005, value=PRESETS[preset.value][1], label="Removal rate k")
    steps = mo.ui.slider(500, 8000, step=500, value=5000, label="Time steps")
    mo.vstack([feed, kill, steps])
    return feed, kill, steps


@app.cell
def _(np):
    def laplacian(X):
        return np.roll(X, 1, -2) + np.roll(X, -1, -2) + np.roll(X, 1, -1) + np.roll(X, -1, -1) - 4 * X

    def gray_scott(F, k, n=128, steps=5000, Du=0.16, Dv=0.08, seed=1):
        """Explicit Euler, five-point Laplacian, periodic boundaries, unit grid and time step.

        F and k may be arrays: each pair then runs as its own dish, all at once.
        """
        F = np.asarray(F, dtype=float)[..., None, None]
        k = np.asarray(k, dtype=float)[..., None, None]
        shape = np.broadcast_shapes(F.shape, k.shape)[:-2] + (n, n)
        rng = np.random.default_rng(seed)
        U = np.ones(shape)
        V = np.zeros(shape)
        for _ in range(10):
            x, y = rng.integers(0, n - 8, 2)
            U[..., x : x + 8, y : y + 8] = 0.5
            V[..., x : x + 8, y : y + 8] = 0.25
        U += 0.02 * rng.random((n, n))
        V += 0.02 * rng.random((n, n))
        for _ in range(steps):
            Lu, Lv = laplacian(U), laplacian(V)
            uvv = U * V * V
            U += Du * Lu - uvv + F * (1 - U)
            V += Dv * Lv + uvv - (F + k) * V
        return U, V

    return (gray_scott,)


@app.cell
def _(feed, gray_scott, kill, mo, plt, steps):
    with mo.status.spinner(title="Reacting and diffusing"):
        _, V = gray_scott(feed.value, kill.value, steps=steps.value)
    fig, ax = plt.subplots(figsize=(5.5, 5.5))
    ax.imshow(V, cmap="magma", interpolation="bicubic")
    ax.set_axis_off()
    ax.set_title(f"v after {steps.value:,} steps, F = {feed.value:.4f}, k = {kill.value:.4f}")
    fig.tight_layout()
    fig
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        ## Where the patterns live

        Pearson mapped which pattern appears for each pair $(F, k)$. Most of the
        plane is dull: either $v$ dies out and the dish stays uniform, or it
        takes over everywhere. Patterns live in a narrow band between the two,
        roughly $0.03 < F < 0.08$ and $0.055 < k < 0.068$. Move a slider a few
        thousandths outside it and the dish goes blank; that sensitivity is
        part of why biologists still argue about which real patterns are
        Turing patterns.
        """
    )
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
        ---
        **Method and limits.** 128 × 128 periodic grid, explicit Euler with
        unit time step and grid spacing, $D_u = 0.16$, $D_v = 0.08$; ten seeded
        squares of $v$ plus 2 percent noise start the reaction. Patterns depend
        on the seed and the grid size, the regimes do not.
        **References.** A. M. Turing, *Phil. Trans. R. Soc. B* 237 (1952) 37;
        P. Gray and S. K. Scott, *Chem. Eng. Sci.* 38 (1983) 29;
        J. E. Pearson, *Science* 261 (1993) 189.
        Fork it to grow a pattern on a sphere or to add a third species.
        """
    )
    return


if __name__ == "__main__":
    app.run()