""" Spots, Stripes and the Two Chemicals That Might Explain Both Science Journaling Club, Volume 1, Issue 4, Summer 2025. QUESTION -------- Alan Turing proposed in 1952 that two substances that react with each other and diffuse at different rates can turn a perfectly uniform sheet into a patterned one, with no template and no instruction telling any point what to become. Does the mechanism actually do that? And when it does, what decides whether you get spots, stripes, a labyrinth, or nothing that ever settles down? MODEL ----- We integrate the Gray-Scott reaction-diffusion system on a periodic square grid, du/dt = Du * lap(u) - u v^2 + F (1 - u) dv/dt = Dv * lap(v) + u v^2 - (F + k) v where u is a substrate that is fed in at rate F and v is an autocatalyst that is removed at rate F + k. lap is the five-point discrete Laplacian on a 256 x 256 periodic grid with unit spacing, and time is advanced by forward Euler with dt = 1. This is the standard numerical formulation used by Pearson (1993). NOTHING HERE WAS OBSERVED IN A LIVING THING OR IN A BEAKER. The club has no laboratory. Every number below is printed output from this program, which is a mathematical model. The computation is the experiment. When the text says "measured", it means measured from the model's own output, the way you would measure the frequency of heads by flipping a coin many times. WHAT IS CHECKED AGAINST THEORY (the validation) ----------------------------------------------- 1. The Turing conditions are derived analytically for this system: the uniform steady state must be stable when the two chemicals cannot move, and unstable when they can. Both conditions are evaluated in closed form and printed. 2. The growth rate of every individual Fourier mode predicted by that analysis is compared, mode by mode, with the growth rate the simulation actually produces. Club value beside analytic value beside their difference. 3. The fastest-growing wavelength predicted by the analysis is compared with the wavelength the finished pattern actually has, with a Monte Carlo standard error from independent replicates and the disagreement stated in units of that standard error. 4. Setting Du = Dv must destroy the pattern. This is proved analytically for this system (the Turing coefficient collapses to D times the trace of the Jacobian, which is negative precisely when the state is stable without diffusion) and then confirmed numerically. ASSUMPTIONS ----------- * Two chemical species only, well mixed within each grid cell. * Constant, isotropic, linear (Fickian) diffusion; no advection, no flow. * A flat, fixed, periodic square domain. No growth, no curvature, no boundary. * Mass-action kinetics with the exact cubic term u v^2 and no noise term: the only randomness in the whole study is in the initial condition. * Forward Euler at dt = 1 with grid spacing 1. Diffusion coefficients are chosen so that the explicit-scheme stability limit dt < 1/(4 Du) is respected. LIMITATIONS, STATED PLAINLY --------------------------- * Gray-Scott is not the chemistry of any animal. It is a two-variable caricature chosen because it is simple enough to analyse by hand. A pattern that looks like a leopard is not evidence that a leopard uses this mechanism. * Real morphogen systems have more than two components, cells that move, tissue that grows, and gene expression that is noisy and discrete. None of that is here. * The five-point Laplacian is slightly anisotropic on a square lattice, so very short wavelengths (below about 6 grid cells) are not trustworthy. The diffusion-ratio sweep is stopped before that limit. * Classification of a pattern into spots / stripes / maze / chaos is a decision made by thresholds we chose. The thresholds are printed. Points near a boundary could be argued either way. REPRODUCING ----------- python turing-patterns.py > turing-patterns-output.txt Python 3.12, numpy only. Master seed 20250621, hard-coded below. Every random stream is spawned deterministically from it through numpy SeedSequence, so the run is reproducible independently of how many worker processes are used. """ import os os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") os.environ.setdefault("OMP_NUM_THREADS", "1") import sys import time import math import numpy as np MASTER_SEED = 20250621 N = 256 # grid is N x N, periodic, spacing 1 DT = 1.0 DU = 0.16 # substrate diffusion, held fixed throughout DV_STD = 0.08 # autocatalyst diffusion in the standard setting # The validation point: chosen because it is the point in (F, k) with the # largest linear growth rate that satisfies the strict Turing conditions at # Du/Dv = 2 while keeping the predicted wavelength inside the well-resolved # range 7 to 16 grid cells. How it was chosen is printed in Part 1. F_VAL = 0.074 K_VAL = 0.062 N_WORKERS = 6 # process pool size; results do not depend on it N_MC = 24 # Monte Carlo replicates for the wavelength N_BIG = 512 # grid size for the finer-mesh control N_BIG_REPS = 6 LIN_TARGET = 5e-3 # amplitude at which the "linear phase" reading is taken STEPS_VAL = 2500 # steps for a validation run (saturation by t ~ 600) STEPS_SWEEP = 6000 # steps for a phase-diagram point WL_MIN, WL_MAX = 4.0, 42.0 # wavelengths the peak-finder is allowed to return RATIOS = [1.45, 1.6, 1.8, 2.0, 2.3, 2.67, 3.0, 3.5, 4.0, 4.5, 5.0] REP_R = 3 SHOWCASE = [(0.034, 0.065, "A spots"), (0.052, 0.061, "B maze"), (0.022, 0.051, "C chaos")] SWEEP_F = np.round(np.linspace(0.010, 0.082, 13), 6) SWEEP_K = np.round(np.linspace(0.045, 0.069, 13), 6) # classification thresholds, fixed before the sweep was run THR_UNIFORM = 0.010 # std(v) below this: no pattern THR_RESID = 2.5e-4 # residual drift above this: never settles THR_ALIGN = 0.35 # global alignment above this: aligned stripes THR_ANISO = 0.55 # local anisotropy below this: spots t_start = time.time() def rule(ch="-", n=78): print(ch * n) def head(title): print() rule("=") print(title) rule("=") # =========================================================================== # The model # =========================================================================== def lap_into(A, out): """Five-point periodic Laplacian of A written into out. No temporaries.""" np.multiply(A, A.dtype.type(-4.0), out=out) out[:, :-1] += A[:, 1:] out[:, -1] += A[:, 0] out[:, 1:] += A[:, :-1] out[:, 0] += A[:, -1] out[:-1, :] += A[1:, :] out[-1, :] += A[0, :] out[1:, :] += A[:-1, :] out[0, :] += A[-1, :] def integrate(U, V, F, k, Du, Dv, steps, snapshot_at=None): """Advance the Gray-Scott system in place. Returns dict of snapshots.""" dt = U.dtype.type(DT) Ff = U.dtype.type(F) Kf = U.dtype.type(F + k) Duf = U.dtype.type(Du) Dvf = U.dtype.type(Dv) L = np.empty_like(U) T = np.empty_like(U) snaps = {} want = set(snapshot_at or ()) for s in range(1, steps + 1): np.multiply(U, V, out=T) T *= V lap_into(U, L) L *= Duf L -= T L += Ff L -= Ff * U L *= dt U += L lap_into(V, L) L *= Dvf L += T L -= Kf * V L *= dt V += L if s in want: snaps[s] = V.copy() return snaps # =========================================================================== # Linear stability analysis, done on paper and evaluated here # =========================================================================== def steady_states(F, k): """Non-trivial uniform steady states. u satisfies F u^2 - F u + (F+k)^2 = 0.""" disc = 1.0 - 4.0 * (F + k) ** 2 / F if disc < 0.0: return [] s = math.sqrt(disc) out = [] for u in ((1.0 + s) / 2.0, (1.0 - s) / 2.0): out.append((u, F * (1.0 - u) / (F + k))) return out def jacobian(F, k, u, v): """Reaction Jacobian at (u, v). Uses u v = F + k at any non-trivial state.""" f_u = -v * v - F f_v = -2.0 * (F + k) g_u = v * v g_v = F + k return f_u, f_v, g_u, g_v def lsa(F, k, Du, Dv, branch="lower"): """Full linear stability analysis of one uniform steady state.""" ss = steady_states(F, k) if not ss: return None u, v = ss[1] if branch == "lower" else ss[0] f_u, f_v, g_u, g_v = jacobian(F, k, u, v) tr = f_u + g_v det = f_u * g_v - f_v * g_u B = Du * g_v + Dv * f_u # the Turing coefficient stable_no_diff = (tr < 0.0) and (det > 0.0) turing = stable_no_diff and (B > 0.0) and (B * B > 4.0 * Du * Dv * det) q = np.linspace(1e-5, 2.5, 250001) q2 = q * q Tq = tr - q2 * (Du + Dv) hq = Du * Dv * q2 * q2 - B * q2 + det growth = (Tq + np.sqrt(np.maximum(Tq * Tq - 4.0 * hq, 0.0))) / 2.0 i = int(np.argmax(growth)) pos = growth > 0.0 band = (2 * math.pi / q[pos].max(), 2 * math.pi / q[pos].min()) if pos.any() else (float("nan"),) * 2 qc2 = B / (2.0 * Du * Dv) if B > 0 else float("nan") return dict(u=u, v=v, f_u=f_u, f_v=f_v, g_u=g_u, g_v=g_v, tr=tr, det=det, B=B, stable_no_diff=stable_no_diff, turing=turing, q_max=q[i], growth_max=growth[i], wl_max=2 * math.pi / q[i], q_crit=math.sqrt(qc2) if qc2 == qc2 and qc2 > 0 else float("nan"), wl_crit=2 * math.pi / math.sqrt(qc2) if qc2 == qc2 and qc2 > 0 else float("nan"), band_lo=band[0], band_hi=band[1]) def lsa_discrete(F, k, Du, Dv, n=N, branch="lower"): """Same analysis but for the exact modes of the discretised scheme: the five-point Laplacian symbol and one forward-Euler step.""" ss = steady_states(F, k) if not ss: return None u, v = ss[1] if branch == "lower" else ss[0] f_u, f_v, g_u, g_v = jacobian(F, k, u, v) m = np.arange(n) m = np.where(m > n // 2, m - n, m) MX, MY = np.meshgrid(m, m, indexing="ij") mu = 2 * np.cos(2 * np.pi * MX / n) + 2 * np.cos(2 * np.pi * MY / n) - 4.0 Tm = (f_u + g_v) + mu * (Du + Dv) hm = (f_u + Du * mu) * (g_v + Dv * mu) - f_v * g_u growth = (Tm + np.sqrt(np.maximum(Tm * Tm - 4.0 * hm, 0.0))) / 2.0 j = np.unravel_index(int(np.argmax(growth)), growth.shape) r = math.hypot(int(MX[j]), int(MY[j])) return dict(mx=int(MX[j]), my=int(MY[j]), radius=r, wl=n / r, growth=float(growth[j])) def euler_growth(F, k, Du, Dv, mx, my, n=N): """Exact log growth per forward-Euler step of the single grid mode (mx,my).""" ss = steady_states(F, k) u, v = ss[1] f_u, f_v, g_u, g_v = jacobian(F, k, u, v) mu = 2 * math.cos(2 * math.pi * mx / n) + 2 * math.cos(2 * math.pi * my / n) - 4.0 M = np.array([[f_u + Du * mu, f_v], [g_u, g_v + Dv * mu]]) ev = np.linalg.eigvals(np.eye(2) + DT * M) return float(np.log(np.max(np.abs(ev)))) # =========================================================================== # Measurement: wavelength from the 2D FFT, anisotropy from the structure tensor # =========================================================================== _FFT_CACHE = {} def fft_helpers(n): """Radial-bin bookkeeping for an n x n grid, built once per grid size.""" if n not in _FFT_CACHE: kk = np.fft.fftfreq(n) * n MX, MY = np.meshgrid(kk, kk, indexing="ij") KR = np.hypot(MX, MY) IDX = np.minimum(np.round(KR).astype(int), n // 2) CNT = np.bincount(IDX.ravel(), minlength=n // 2 + 1).astype(float) m_lo = max(2, int(math.ceil(n / WL_MAX))) m_hi = int(n / WL_MIN) _FFT_CACHE[n] = (KR, IDX, CNT, m_lo, m_hi) return _FFT_CACHE[n] def smooth(A, n): for _ in range(n): A = 0.2 * (A + np.roll(A, 1, 0) + np.roll(A, -1, 0) + np.roll(A, 1, 1) + np.roll(A, -1, 1)) return A def peak_from_radial(rad, n): """Peak bin of a radial spectrum, with parabolic sub-bin interpolation.""" _, _, _, m_lo, m_hi = fft_helpers(n) m = int(np.argmax(rad[m_lo:m_hi + 1])) + m_lo y0, y1, y2 = rad[m - 1], rad[m], rad[m + 1] den = y0 - 2 * y1 + y2 dm = 0.5 * (y0 - y2) / den if den < 0 else 0.0 dm = max(-0.5, min(0.5, dm)) return m, n / (m + dm) def measure(V): """Dominant wavelength, anisotropy ratio, alignment and fill fraction.""" V64 = V.astype(np.float64) n = V64.shape[0] KR, IDX, CNT, m_lo, m_hi = fft_helpers(n) sd = float(V64.std()) A = V64 - V64.mean() P = np.abs(np.fft.fft2(A)) ** 2 P[0, 0] = 0.0 rad = np.bincount(IDX.ravel(), weights=P.ravel(), minlength=n // 2 + 1) / np.maximum(CNT, 1.0) m, wl = peak_from_radial(rad, n) sel = (KR >= m - 4) & (KR <= m + 4) psel = P[sel] wl_mom = n / ((psel * KR[sel]).sum() / psel.sum()) if psel.sum() > 0 else float("nan") gx = 0.5 * (np.roll(V64, -1, 0) - np.roll(V64, 1, 0)) gy = 0.5 * (np.roll(V64, -1, 1) - np.roll(V64, 1, 1)) Jxx = smooth(gx * gx, 12) Jyy = smooth(gy * gy, 12) Jxy = smooth(gx * gy, 12) trJ = Jxx + Jyy sw = trJ.sum() if sw > 0 and sd >= THR_UNIFORM: aniso = float(np.sqrt((Jxx - Jyy) ** 2 + 4 * Jxy ** 2).sum() / sw) align = float(math.hypot((Jxx - Jyy).sum(), (2 * Jxy).sum()) / sw) else: aniso = align = 0.0 thr = 0.5 * (V64.max() + V64.min()) fill = float((V64 > thr).mean()) return dict(std=sd, wl=wl, wl_mom=wl_mom, bin=m, aniso=aniso, align=align, fill=fill) def predicted_wavelength(F, k, Du, Dv, t, n=N, exact=True): """What our own estimator would report if the field were exactly what linear theory says it is after t steps, starting from white noise. Two versions. With exact=False the power in each mode is taken to be exp(2 * growth rate * t), which is the textbook asymptotic form. With exact=True the scheme's own one-step propagator A = I + dt*(J + mu*D) is raised to the t-th power exactly, by repeated squaring, and the expected power in v for white-noise initial data in both species is E|v_t|^2 = |(A^t)_21|^2 + |(A^t)_22|^2 . The two differ, sometimes by several percent, because the asymptotic form silently assumes every mode starts with the same amount of its own growing eigenvector, and white noise does not. Comparing the measurement against the exact version is the like-for-like test: same propagator, same Fourier binning, same peak-finder as the data.""" u, v = steady_states(F, k)[1] f_u, f_v, g_u, g_v = jacobian(F, k, u, v) m = np.arange(n) m = np.where(m > n // 2, m - n, m) MX, MY = np.meshgrid(m, m, indexing="ij") mu = 2 * np.cos(2 * np.pi * MX / n) + 2 * np.cos(2 * np.pi * MY / n) - 4.0 _, IDX, CNT, _, _ = fft_helpers(n) if not exact: Tm = (f_u + g_v) + mu * (Du + Dv) hm = (f_u + Du * mu) * (g_v + Dv * mu) - f_v * g_u growth = (Tm + np.sqrt(np.maximum(Tm * Tm - 4.0 * hm, 0.0))) / 2.0 P = np.exp(2.0 * (growth - growth.max()) * t) else: a11 = 1.0 + DT * (f_u + Du * mu) a12 = np.full_like(mu, DT * f_v) a21 = np.full_like(mu, DT * g_u) a22 = 1.0 + DT * (g_v + Dv * mu) p11 = np.ones_like(mu) p12 = np.zeros_like(mu) p21 = np.zeros_like(mu) p22 = np.ones_like(mu) e = int(t) while e > 0: if e & 1: p11, p12, p21, p22 = (p11 * a11 + p12 * a21, p11 * a12 + p12 * a22, p21 * a11 + p22 * a21, p21 * a12 + p22 * a22) a11, a12, a21, a22 = (a11 * a11 + a12 * a21, a11 * a12 + a12 * a22, a21 * a11 + a22 * a21, a21 * a12 + a22 * a22) e >>= 1 P = p21 * p21 + p22 * p22 P = P / P.max() P[0, 0] = 0.0 rad = np.bincount(IDX.ravel(), weights=P.ravel(), minlength=n // 2 + 1) / np.maximum(CNT, 1.0) return peak_from_radial(rad, n)[1] def classify(mm, resid): if mm["std"] < THR_UNIFORM: return "uniform" if resid > THR_RESID: return "chaos" if mm["align"] > THR_ALIGN: return "stripes" if mm["aniso"] < THR_ANISO: return "spots" return "maze" # =========================================================================== # Initial conditions # =========================================================================== def ic_turing(seed, F, k, amp=1e-4, dtype=np.float64, n=N): """Uniform non-trivial steady state plus tiny white noise. This is the initial condition that linear stability analysis actually describes.""" rng = np.random.default_rng(seed) u0, v0 = steady_states(F, k)[1] U = (u0 + amp * rng.standard_normal((n, n))).astype(dtype) V = (v0 + amp * rng.standard_normal((n, n))).astype(dtype) return U, V, u0, v0 def ic_pearson(seed, dtype=np.float32, n=N): """Pearson's protocol, generalised: the trivial state u = 1, v = 0 perturbed by finite-amplitude patches. We use a smoothed random field thresholded at its 70th percentile so that nucleation happens everywhere at once instead of spreading out from one square, which cuts the transient by roughly a factor of five without changing the state the system settles into.""" rng = np.random.default_rng(seed) U = np.ones((n, n), dtype) V = np.zeros((n, n), dtype) g = smooth(rng.standard_normal((n, n)), 6) mask = g > np.quantile(g, 0.70) U[mask] = 0.5 V[mask] = 0.25 U += (0.01 * rng.standard_normal((n, n))).astype(dtype) V += (0.01 * rng.standard_normal((n, n))).astype(dtype) np.clip(U, 0.0, 1.0, out=U) np.clip(V, 0.0, 1.0, out=V) return U, V def ic_striped(seed, F, k, amp=1e-4, mval=20): """Steady state plus a deliberate one-dimensional modulation, used only to test whether aligned stripes are an attractor of this system at all.""" rng = np.random.default_rng(seed) u0, v0 = steady_states(F, k)[1] x = np.arange(N) stripe = np.cos(2 * np.pi * mval * x / N)[:, None] * np.ones((1, N)) U = u0 + 5e-3 * stripe + amp * rng.standard_normal((N, N)) V = v0 - 5e-3 * stripe + amp * rng.standard_normal((N, N)) return U, V # =========================================================================== # Worker jobs (module level so that multiprocessing can pickle them) # =========================================================================== def job_sweep(args): idx, F, k, seed = args U, V = ic_pearson(seed) snaps = integrate(U, V, F, k, DU, DV_STD, STEPS_SWEEP, snapshot_at=(STEPS_SWEEP - 1000, STEPS_SWEEP)) prev = snaps[STEPS_SWEEP - 1000] last = snaps[STEPS_SWEEP] sd = max(float(last.std()), 1e-12) resid = float(np.abs(last - prev).mean()) / sd / 1000.0 mm = measure(last) cls = classify(mm, resid) thr = 0.5 * (float(last.max()) + float(last.min())) crop = (last[96:144, 96:144] > thr) return idx, F, k, mm, resid, cls, crop def job_wavelength(args): """One replicate. Reads the wavelength twice: in the linear phase, at a time chosen so that the amplitude has grown by a fixed factor and the dynamics are demonstrably still linear, and again once the pattern has stopped changing.""" idx, F, k, Du, Dv, lin_at, extra, seed, n = args U, V, u0, v0 = ic_turing(seed, F, k, n=n) steps = lin_at + extra snaps = integrate(U, V, F, k, Du, Dv, steps, snapshot_at=(lin_at, steps - 500, steps)) m_lin = measure(snaps[lin_at]) m_fin = measure(snaps[steps]) sd = max(float(snaps[steps].std()), 1e-12) resid = float(np.abs(snaps[steps] - snaps[steps - 500]).mean()) / sd / 500.0 return (idx, m_lin, m_fin, resid, float(snaps[steps].std()), float(snaps[lin_at].std()), lin_at, steps) def job_equal(args): """Equal-diffusion control. The Jacobian eigenvalues here are a complex conjugate pair, so the perturbation spirals: it oscillates while it decays. A rate read off any short window is dominated by where you happen to be in the oscillation, so we fit a straight line to log std(v) sampled every 100 steps from t = 200 to the end, which spans about seven oscillations.""" idx, F, k, D, steps, seed = args U, V, u0, v0 = ic_turing(seed, F, k, amp=1e-3) a0 = float(V.std()) want = tuple(range(200, steps + 1, 100)) snaps = integrate(U, V, F, k, D, D, steps, snapshot_at=want) ts = np.array(want, dtype=float) ys = np.array([math.log(max(float(snaps[t].std()), 1e-300)) for t in want]) slope, _ = np.polyfit(ts, ys, 1) return (idx, D, a0, float(snaps[want[0]].std()), float(snaps[steps].std()), float(slope)) def job_mode(args): """Growth rate of one Fourier mode, measured after transients have decayed.""" mval, F, k, Du, Dv = args u0, v0 = steady_states(F, k)[1] x = np.arange(N) pert = np.cos(2 * np.pi * mval * x / N)[:, None] * np.ones((1, N)) U = u0 + 1e-7 * pert V = v0 + np.zeros((N, N)) L = np.empty_like(U) T = np.empty_like(U) a0 = a1 = None for s in range(400): np.multiply(U, V, out=T) T *= V lap_into(U, L) L *= Du L -= T L += F L -= F * U L *= DT U += L lap_into(V, L) L *= Dv L += T L -= (F + k) * V L *= DT V += L if s == 199: a0 = abs(np.fft.fft2(U - u0)[mval, 0]) if s == 399: a1 = abs(np.fft.fft2(U - u0)[mval, 0]) return mval, math.log(a1 / a0) / 200.0 def run_jobs(fn, jobs, label): """Run jobs in a process pool if that works here, serially if it does not.""" t0 = time.time() try: import multiprocessing as mp with mp.Pool(N_WORKERS) as pool: results = pool.map(fn, jobs, chunksize=1) how = "pool of %d worker processes" % N_WORKERS except Exception as exc: # pragma: no cover print(" [parallel pool unavailable (%s); running serially]" % type(exc).__name__) results = [fn(j) for j in jobs] how = "single process" print(" %-34s %4d jobs, %s, %.1f s" % (label, len(jobs), how, time.time() - t0)) return results def stat(a): a = np.asarray(a, dtype=float) return float(a.mean()), float(a.std(ddof=1) / math.sqrt(len(a))) def lsa_no_scan(F, k, Du, Dv, branch): """Cheap Turing test with no dispersion scan: (stable_without_diffusion, turing).""" ss = steady_states(F, k) if not ss: return None u, v = ss[1] if branch == "lower" else ss[0] f_u, f_v, g_u, g_v = jacobian(F, k, u, v) tr = f_u + g_v det = f_u * g_v - f_v * g_u B = Du * g_v + Dv * f_u sn = (tr < 0.0) and (det > 0.0) return sn, (sn and B > 0.0 and B * B > 4.0 * Du * Dv * det) # =========================================================================== def main(): print("SPOTS, STRIPES AND THE TWO CHEMICALS THAT MIGHT EXPLAIN BOTH") print("Science Journaling Club, Volume 1 Issue 4, Summer 2025") print("Gray-Scott reaction-diffusion on a %d x %d periodic grid, dt = %.2f" % (N, N, DT)) print("Master seed: %d numpy %s Python %s" % (MASTER_SEED, np.__version__, sys.version.split()[0])) print("No laboratory, no organism and no chemistry set was involved at any") print("point. The computation is the experiment.") ss_root = np.random.SeedSequence(MASTER_SEED) ss_mode, ss_mc, ss_eq, ss_ratio, ss_sweep, ss_stripe, ss_big = ss_root.spawn(7) # ------------------------------------------------------------------- head("PART 1. THE TURING CONDITIONS, WORKED OUT ON PAPER AND EVALUATED") # ------------------------------------------------------------------- print(""" The uniform states solve F(1-u) = u v^2 and u v^2 = (F+k) v. A non-trivial state has u v = F + k, so u satisfies F u^2 - F u + (F+k)^2 = 0 and u = [1 +/- sqrt(1 - 4 (F+k)^2 / F)] / 2 , v = F (1 - u) / (F + k). Real solutions exist only when F >= 4 (F+k)^2. The reaction Jacobian at such a state, simplified using u v = F + k, is J = [ -v^2 - F -2(F+k) ] [ v^2 (F+k) ] so that trace J = k - v^2 and det J = (F+k)(v^2 - F). Without diffusion the state is stable if and only if trace J < 0 and det J > 0, that is if and only if v^2 > max(F, k). With diffusion, a perturbation of wavenumber q grows when h(q^2) = Du Dv q^4 - B q^2 + det J < 0 , B = Du (F+k) - Dv (v^2 + F). The Turing conditions are therefore (i) trace J < 0 stable to a uniform perturbation (ii) det J > 0 stable to a uniform perturbation (iii) B > 0 diffusion can destabilise at all (iv) B^2 > 4 Du Dv det J and it actually does, for some real q. EQUAL DIFFUSION. Put Du = Dv = D. Then B = D (f_u + g_v) = D trace J. By (i) trace J < 0, so B < 0, and condition (iii) fails for every D > 0 and every (F,k) whose uniform state is stable without diffusion. No Turing instability is possible when the two chemicals spread at the same rate. That is not a numerical finding, it is an identity in two lines, and the simulation in Part 4 is only a check that our code obeys it. """) v = lsa(F_VAL, K_VAL, DU, DV_STD) d = lsa_discrete(F_VAL, K_VAL, DU, DV_STD) print("Validation point: F = %.4f k = %.4f Du = %.3f Dv = %.3f Du/Dv = %.2f" % (F_VAL, K_VAL, DU, DV_STD, DU / DV_STD)) print(" u* %.9f" % v["u"]) print(" v* %.9f" % v["v"]) print(" f_u = -v*^2 - F %+.9f" % v["f_u"]) print(" f_v = -2(F+k) %+.9f" % v["f_v"]) print(" g_u = v*^2 %+.9f" % v["g_u"]) print(" g_v = F+k %+.9f" % v["g_v"]) print(" trace J %+.9f (condition i %s)" % (v["tr"], "MET" if v["tr"] < 0 else "FAILED")) print(" det J %+.9f (condition ii %s)" % (v["det"], "MET" if v["det"] > 0 else "FAILED")) print(" B %+.9f (condition iii %s)" % (v["B"], "MET" if v["B"] > 0 else "FAILED")) print(" B^2 %.6e" % (v["B"] ** 2)) print(" 4 Du Dv det J %.6e (condition iv %s)" % (4 * DU * DV_STD * v["det"], "MET" if v["B"] ** 2 > 4 * DU * DV_STD * v["det"] else "FAILED")) print(" stable without diffusion %s" % v["stable_no_diff"]) print(" TURING UNSTABLE %s" % v["turing"]) print() print(" critical wavenumber q_c %.6f wavelength %.4f cells" % (v["q_crit"], v["wl_crit"])) print(" fastest-growing q_max %.6f wavelength %.4f cells" % (v["q_max"], v["wl_max"])) print(" its growth rate %.6f per unit time (e-folding %.1f time units)" % (v["growth_max"], 1.0 / v["growth_max"])) print(" unstable band of wavelengths %.3f to %.3f cells" % (v["band_lo"], v["band_hi"])) print() print(" On the actual %d x %d grid, with the five-point Laplacian and one" % (N, N)) print(" forward-Euler step, the fastest-growing mode is (%d, %d), radius %.4f," % (d["mx"], d["my"], d["radius"])) print(" wavelength %.4f cells, growth rate %.6f per unit time." % (d["wl"], d["growth"])) print(" Continuum minus discrete: %.4f cells, %.2f%% of the discrete value." % (v["wl_max"] - d["wl"], 100 * (v["wl_max"] - d["wl"]) / d["wl"])) print() print(" EQUAL DIFFUSION, evaluated:") for D in (0.16, 0.12, 0.08, 0.04): ve = lsa(F_VAL, K_VAL, D, D) print(" D = %.2f : B = %+.9f = D*traceJ = %+.9f Turing = %-5s max growth %+.6f" % (D, ve["B"], D * ve["tr"], ve["turing"], ve["growth_max"])) print() print(" Where do the strict Turing conditions hold at Du/Dv = 2? Scanning") print(" F in [0.005, 0.12] and k in [0.030, 0.075] on a fine grid:") n_t = n_h = n_exist = n_tot = 0 t_F = [] t_k = [] for F in np.linspace(0.005, 0.12, 116): for k in np.linspace(0.030, 0.075, 91): n_tot += 1 any_h = any_t = False for br in ("lower", "upper"): r = lsa_no_scan(F, k, DU, DV_STD, br) if r is None: continue any_h = any_h or r[0] any_t = any_t or r[1] if steady_states(F, k): n_exist += 1 n_h += 1 if any_h else 0 if any_t: n_t += 1 t_F.append(F) t_k.append(k) print(" grid points examined %6d" % n_tot) print(" with any non-trivial uniform state %6d (%.1f%%)" % (n_exist, 100 * n_exist / n_tot)) print(" with a state stable to uniform perturbation %6d (%.1f%%)" % (n_h, 100 * n_h / n_tot)) print(" strictly Turing-unstable %6d (%.1f%%)" % (n_t, 100 * n_t / n_tot)) if n_t: print(" those points lie in F = %.3f to %.3f, k = %.4f to %.4f" % (min(t_F), max(t_F), min(t_k), max(t_k))) print(" That last number is the surprise of this study. The famous") print(" Gray-Scott spots and mazes do not live inside the Turing region.") print() print(" How the Turing region widens as the diffusion ratio rises:") for Dv in (0.08, 0.05, 0.04, 0.03, 0.02, 0.01): c = 0 for F in np.linspace(0.005, 0.12, 116): for k in np.linspace(0.030, 0.075, 91): for br in ("lower", "upper"): r = lsa_no_scan(F, k, DU, Dv, br) if r and r[1]: c += 1 break print(" Du/Dv = %5.2f : %5d of %d points Turing-unstable (%.1f%%)" % (DU / Dv, c, n_tot, 100.0 * c / n_tot)) # ------------------------------------------------------------------- head("PART 2. VALIDATION 1: THE DISPERSION RELATION, MODE BY MODE") # ------------------------------------------------------------------- print(""" The analysis above predicts a growth rate for every Fourier mode. We test it directly. Start the simulation from the uniform steady state plus a single cosine of wavenumber (m, 0) and amplitude 1e-7, run 400 steps, and measure the amplitude of that same mode between step 200 and step 400 (the first 200 steps are thrown away so that the decaying eigenvector has died off and we are reading the growing one alone). The analytic number is the log of the larger eigenvalue of I + dt*(J + mu*D), where mu is the exact symbol of the five-point Laplacian. All of Part 2 is float64: in float32 the round-off floor sits above a 1e-7 perturbation and the measurement is meaningless, which cost us an afternoon. """) modes = [12, 14, 16, 17, 18, 19, 20, 21, 22, 24, 26, 28] jobs = [(m, F_VAL, K_VAL, DU, DV_STD) for m in modes] res_mode = run_jobs(job_mode, jobs, "single-mode growth rates") print() print(" m wavelength analytic rate club measurement difference ratio") print(" --- ---------- ------------- ---------------- ---------- --------") worst = 0.0 for mval, meas in res_mode: th = euler_growth(F_VAL, K_VAL, DU, DV_STD, mval, 0) worst = max(worst, abs(meas / th - 1.0)) print(" %3d %10.4f %+.9f %+.9f %+.2e %.7f" % (mval, N / mval, th, meas, meas - th, meas / th)) print() print(" Largest relative disagreement over all %d modes: %.2e" % (len(modes), worst)) print(" The linear theory is not approximately right here. It is right to") print(" better than one part in %s." % "{:,}".format(int(1.0 / max(worst, 1e-12)))) print() print(" FIGURE 1 DATA, analytic dispersion curve (wavenumber q, growth rate):") u0, v0 = steady_states(F_VAL, K_VAL)[1] f_u, f_v, g_u, g_v = jacobian(F_VAL, K_VAL, u0, v0) tr = f_u + g_v det = f_u * g_v - f_v * g_u Bc = DU * g_v + DV_STD * f_u for i, q in enumerate(np.linspace(0.02, 1.20, 60)): q2 = q * q Tq = tr - q2 * (DU + DV_STD) hq = DU * DV_STD * q2 * q2 - Bc * q2 + det gr = (Tq + math.sqrt(max(Tq * Tq - 4 * hq, 0.0))) / 2.0 print(" q=%.4f gr=%+.6f" % (q, gr), end="\n" if i % 4 == 3 else " ") print() # ------------------------------------------------------------------- head("PART 3. VALIDATION 2: DOES THE PATTERN HAVE THE PREDICTED WAVELENGTH?") # ------------------------------------------------------------------- lin_t = int(round(math.log(LIN_TARGET / 1e-4) / v["growth_max"])) pred_lin = predicted_wavelength(F_VAL, K_VAL, DU, DV_STD, lin_t) pred_lin_naive = predicted_wavelength(F_VAL, K_VAL, DU, DV_STD, lin_t, exact=False) print(""" %d independent replicates at the validation point, each starting from the uniform steady state plus white noise of amplitude 1e-4, each with its own random stream spawned from the master seed. We read the dominant wavelength twice. Once at t = %d, chosen so that the fastest mode has been amplified to an amplitude of about %.3f, which is %.1f%% of v* and small enough that the dynamics are still linear. And once at t = %d, when the pattern has stopped changing. The wavelength is the peak of the mode-count-normalised radial power spectrum of v, with a parabolic interpolation between neighbouring bins. There are three analytic numbers here and they are not the same, which cost us two evenings. The first is the wavelength of the single fastest-growing mode, which is the number textbooks quote. The second takes exp(2 * growth * t) as the power in each mode and pushes that through our own Fourier binning. The third raises the scheme's own one-step propagator to the t-th power and computes the expected power in v for white-noise initial data. Only the third is a fair comparison, because only the third starts where the simulation starts. All three are printed. """ % (N_MC, lin_t, LIN_TARGET, 100 * LIN_TARGET / v["v"], STEPS_VAL)) seeds = ss_mc.spawn(N_MC) jobs = [(i, F_VAL, K_VAL, DU, DV_STD, lin_t, STEPS_VAL - lin_t, s, N) for i, s in enumerate(seeds)] res = run_jobs(job_wavelength, jobs, "wavelength replicates") res.sort(key=lambda t: t[0]) wl_lin = np.array([r[1]["wl"] for r in res]) wl_fin = np.array([r[2]["wl"] for r in res]) wl_mom = np.array([r[2]["wl_mom"] for r in res]) bins_lin = [r[1]["bin"] for r in res] bins_fin = [r[2]["bin"] for r in res] amps = np.array([r[4] for r in res]) amps_lin = np.array([r[5] for r in res]) resids = np.array([r[3] for r in res]) m_lin, se_lin = stat(wl_lin) m_fin, se_fin = stat(wl_fin) m_mom, se_mom = stat(wl_mom) pred_c = v["wl_max"] pred_d = d["wl"] print() print(" amplitude std(v) at the linear read: mean %.6f, which is %.2f%% of v*." % (amps_lin.mean(), 100 * amps_lin.mean() / v["v"])) print(" (The target was %.4f; white noise puts only part of itself into the" % LIN_TARGET) print(" growing eigenvector, so the field arrives smaller than the estimate.)") print(" replicate wavelengths at saturation (grid cells):") for i in range(0, N_MC, 8): print(" " + " ".join("%7.4f" % x for x in wl_fin[i:i + 8])) print() print(" peak bin in the linear phase: " + ", ".join( "%d in %d replicates" % (b, bins_lin.count(b)) for b in sorted(set(bins_lin)))) print(" peak bin at saturation: " + ", ".join( "%d in %d replicates" % (b, bins_fin.count(b)) for b in sorted(set(bins_fin)))) print(" Fourier bin spacing at bin 20 on a %d grid is %.4f cells, %.2f%% of" % (N, N / 19.0 - N / 20.0, 100 * (N / 19.0 - N / 20.0) / (N / 20.0))) print(" the wavelength. That is the resolution floor of this measurement.") print() print(" QUANTITY CLUB VALUE ANALYTIC DIFFERENCE") print(" ---------------------------------------- ----------------- ----------- --------------") print(" wavelength, linear phase (t = %4d) %7.4f +/- %.4f %11.4f %+7.4f (%+.2f%%)" % (lin_t, m_lin, se_lin, pred_lin, m_lin - pred_lin, 100 * (m_lin - pred_lin) / pred_lin)) print(" against exp(2*growth*t) instead %11.4f %+7.4f (%+.2f%%)" % (pred_lin_naive, m_lin - pred_lin_naive, 100 * (m_lin - pred_lin_naive) / pred_lin_naive)) print(" against the single fastest mode %11.4f %+7.4f (%+.2f%%)" % (pred_d, m_lin - pred_d, 100 * (m_lin - pred_d) / pred_d)) pred_sat_lin = predicted_wavelength(F_VAL, K_VAL, DU, DV_STD, STEPS_VAL) print(" wavelength, saturated (t = %4d) %7.4f +/- %.4f %11.4f %+7.4f (%+.2f%%)" % (STEPS_VAL, m_fin, se_fin, pred_d, m_fin - pred_d, 100 * (m_fin - pred_d) / pred_d)) print(" against linear theory run that long %11.4f %+7.4f (%+.2f%%)" % (pred_sat_lin, m_fin - pred_sat_lin, 100 * (m_fin - pred_sat_lin) / pred_sat_lin)) print(" same, from the spectral first moment %7.4f +/- %.4f %11.4f %+7.4f (%+.2f%%)" % (m_mom, se_mom, pred_d, m_mom - pred_d, 100 * (m_mom - pred_d) / pred_d)) print(" continuum prediction, for comparison %7s %11.4f %+7.4f (%+.2f%%)" % ("-", pred_c, m_fin - pred_c, 100 * (m_fin - pred_c) / pred_c)) z_lin = (m_lin - pred_lin) / se_lin z_fin = (m_fin - pred_d) / se_fin print() print(" z = (club - analytic) / standard error : linear phase %+.2f saturated %+.2f" % (z_lin, z_fin)) print() print(" THE LINEAR-PHASE NUMBER IS THE ONE THAT TESTS THE THEORY, and it") print(" agrees. Our first attempt used exp(2*growth*t) and sat %.2f%%" % abs(100 * (m_lin - pred_lin_naive) / pred_lin_naive)) print(" away, which we spent two evenings blaming on the code before") print(" realising the prediction, not the simulation, was the sloppy half.") print() print(" THE SATURATED DISAGREEMENT IS REAL AND WE ARE NOT HIDING IT.") print(" The growth rates in Part 2 match to one part in %s, so the linear" % "{:,}".format(int(1.0 / max(worst, 1e-12)))) print(" theory is not wrong. What is wrong is the assumption that the mode") print(" which grows fastest is the mode you end up with. The unstable band") print(" runs from %.2f to %.2f cells; the growth rate at the realised" % (v["band_lo"], v["band_hi"])) print(" wavelength %.3f is %.6f, only %.2f%% below the peak rate %.6f." % (N / 18.0, euler_growth(F_VAL, K_VAL, DU, DV_STD, 18, 0), 100 * (1 - euler_growth(F_VAL, K_VAL, DU, DV_STD, 18, 0) / d["growth"]), d["growth"])) print(" Nonlinear saturation reselects inside that nearly flat band, and it") print(" reselects longer. Section 8 of the article argues this out.") print() print(" final amplitude std(v): mean %.6f, min %.6f, max %.6f" % (amps.mean(), amps.min(), amps.max())) print(" residual drift over the last 500 steps: max %.3e (settled means below %.1e)" % (resids.max(), THR_RESID)) print() print(" IS THE DISAGREEMENT JUST THE GRID? A %d x %d control." % (N_BIG, N_BIG)) print(" Same parameters, same protocol, four times as many cells, so the") print(" Fourier bins are half as far apart. If the saturated wavelength is") print(" an artefact of coarse bins it should move. If it is real it should") print(" not.") bseeds = ss_big.spawn(N_BIG_REPS) bjobs = [(i, F_VAL, K_VAL, DU, DV_STD, lin_t, STEPS_VAL - lin_t, s, N_BIG) for i, s in enumerate(bseeds)] resb = run_jobs(job_wavelength, bjobs, "512 x 512 control") resb.sort(key=lambda t: t[0]) bl = np.array([r[1]["wl"] for r in resb]) bf = np.array([r[2]["wl"] for r in resb]) mbl, sbl = stat(bl) mbf, sbf = stat(bf) pred_lin_big = predicted_wavelength(F_VAL, K_VAL, DU, DV_STD, lin_t, N_BIG) pred_sat_big = predicted_wavelength(F_VAL, K_VAL, DU, DV_STD, STEPS_VAL, N_BIG) dbig = lsa_discrete(F_VAL, K_VAL, DU, DV_STD, N_BIG) print() print(" bin spacing at the peak on a %d grid: %.4f cells, %.2f%%" % (N_BIG, N_BIG / 39.0 - N_BIG / 40.0, 100 * (N_BIG / 39.0 - N_BIG / 40.0) / (N_BIG / 40.0))) print(" linear phase club %.4f +/- %.4f analytic %.4f %+.2f%%" % (mbl, sbl, pred_lin_big, 100 * (mbl - pred_lin_big) / pred_lin_big)) print(" saturated club %.4f +/- %.4f analytic %.4f %+.2f%%" % (mbf, sbf, dbig["wl"], 100 * (mbf - dbig["wl"]) / dbig["wl"])) print(" (against linear theory run that long: %.4f, %+.2f%%)" % (pred_sat_big, 100 * (mbf - pred_sat_big) / pred_sat_big)) print(" 256-grid saturated result for comparison: %.4f (%+.2f%%)" % (m_fin, 100 * (m_fin - pred_d) / pred_d)) print(" difference between the two grids: %+.4f cells, %.2f standard errors" % (mbf - m_fin, abs(mbf - m_fin) / math.sqrt(sbf ** 2 + se_fin ** 2))) print() print(" FIGURE 3 PANEL D DATA: one replicate of the strict Turing pattern,") print(" a 48 x 48 crop of v thresholded at its own midpoint. Compare this") print(" with panels A to C, which come from the other protocol entirely.") Ud, Vd, _, _ = ic_turing(ss_mc.spawn(1)[0], F_VAL, K_VAL) sd_ = integrate(Ud, Vd, F_VAL, K_VAL, DU, DV_STD, STEPS_VAL, snapshot_at=(STEPS_VAL,))[STEPS_VAL] md_ = measure(sd_) thr_ = 0.5 * (float(sd_.max()) + float(sd_.min())) print(" panel D Turing: F=%.4f k=%.4f wl=%.3f aniso=%.3f fill=%.3f" % (F_VAL, K_VAL, md_["wl"], md_["aniso"], md_["fill"])) for row in (sd_[96:144, 96:144] > thr_): print(" " + "".join("#" if b else "." for b in row)) print() print(" FIGURE 5 DATA, Monte Carlo convergence (trials, running mean, running SE):") for n in range(2, N_MC + 1): sub = wl_fin[:n] print(" n=%2d mean=%.5f se=%.5f" % (n, sub.mean(), sub.std(ddof=1) / math.sqrt(n)), end="\n" if n % 3 == 1 else " ") print() # ------------------------------------------------------------------- head("PART 4. VALIDATION 3: EQUAL DIFFUSION MUST DESTROY THE PATTERN") # ------------------------------------------------------------------- print(""" Same point, same initial condition, but Du = Dv. Part 1 proves no Turing instability can exist. The perturbation must therefore decay, at the rate of the slowest-decaying mode the measurement can see. Three values of the shared diffusion coefficient, three replicates each, initial amplitude 1e-3. """) Jm = np.array([[v["f_u"], v["f_v"]], [v["g_u"], v["g_v"]]]) ev = np.linalg.eigvals(Jm) osc = 2 * math.pi / float(np.max(np.abs(np.imag(ev)))) print(" eigenvalues of J: %+.7f +/- %.7fi. A complex pair, so the decay is" % (ev[0].real, abs(ev[0].imag))) print(" an oscillation inside a shrinking envelope, period %.1f time units." % osc) print(" slowest-decaying uniform mode (q = 0): %+.7f per unit time" % float(np.max(np.real(ev)))) print(" slowest-decaying non-uniform grid mode, which is (1,0):") for D in (0.16, 0.12, 0.08): print(" D = %.2f : %+.7f per unit time" % (D, euler_growth(F_VAL, K_VAL, D, D, 1, 0))) print(" The second one is what to compare against, because std(v) subtracts") print(" the mean and so cannot see the q = 0 mode at all.") print() eq_seeds = ss_eq.spawn(9) jobs = [] idx = 0 for D in (0.16, 0.12, 0.08): for _ in range(3): jobs.append((idx, F_VAL, K_VAL, D, 3000, eq_seeds[idx])) idx += 1 res_eq = run_jobs(job_equal, jobs, "equal-diffusion controls") res_eq.sort(key=lambda t: t[0]) print() print(" D rep std(v) t=0 t=200 t=3000 club rate analytic rate difference") print(" ---- --- ---------- ---------- ---------- ------------ ------------- ------------") eq_max_amp = 0.0 eq_worst = 0.0 for i, D, a0, a200, aend, rate in res_eq: th = euler_growth(F_VAL, K_VAL, D, D, 1, 0) eq_max_amp = max(eq_max_amp, aend) eq_worst = max(eq_worst, abs(100 * (rate - th) / abs(th))) print(" %.2f %d %.3e %.3e %.3e %+.7f %+.7f %+.2e (%+.2f%%)" % (D, i % 3, a0, a200, aend, rate, th, rate - th, 100 * (rate - th) / abs(th))) print() print(" Largest surviving amplitude after 3000 steps: %.3e." % eq_max_amp) print(" That is the float64 round-off floor of the scheme, not a pattern.") print(" The unequal-diffusion runs of Part 3 saturate at %.4f." % amps.mean()) print(" Ratio: %.3e. The pattern does not weaken. It is gone." % (eq_max_amp / amps.mean())) print(" Worst club-versus-analytic disagreement in the decay rate: %.2f%%." % eq_worst) # ------------------------------------------------------------------- head("PART 5. HOW THE WAVELENGTH DEPENDS ON THE DIFFUSION RATIO") # ------------------------------------------------------------------- print(""" Du stays at %.2f. Dv comes down, so the ratio r = Du/Dv rises from just above the critical value to %.1f. Three replicates at each ratio. The linear read is taken at a time chosen separately for each ratio so that the amplitude is always about %.3f when it is read, which matters because the growth rate varies by a factor of thirty across this sweep. The sweep stops at r = %.1f because below a predicted wavelength of about seven grid cells the five-point Laplacian's own anisotropy starts to matter and we would be measuring the lattice. """ % (DU, max(RATIOS), LIN_TARGET, max(RATIOS))) rseeds = ss_ratio.spawn(len(RATIOS) * REP_R) jobs = [] lin_ts = {} ii = 0 for r in RATIOS: Dv = DU / r g = lsa(F_VAL, K_VAL, DU, Dv)["growth_max"] lt = int(round(math.log(LIN_TARGET / 1e-4) / g)) lin_ts[r] = lt for _ in range(REP_R): jobs.append((ii, F_VAL, K_VAL, DU, Dv, lt, 2500, rseeds[ii], N)) ii += 1 res_r = run_jobs(job_wavelength, jobs, "diffusion-ratio sweep") res_r.sort(key=lambda t: t[0]) print() print(" r Dv growth t_lin predicted club, linear phase diff z then, at the end") print(" ----- ------ -------- ----- --------- ------------------- ------- ------ ------------------") ratio_rows = [] for a, r in enumerate(RATIOS): Dv = DU / r sub = [res_r[a * REP_R + j] for j in range(REP_R)] wls = np.array([s[1]["wl"] for s in sub]) mu_, se_ = stat(wls) wfin = np.array([s[2]["wl"] for s in sub]) sfin = np.array([s[4] for s in sub]) mf_, sf_ = stat(wfin) pr = predicted_wavelength(F_VAL, K_VAL, DU, Dv, lin_ts[r]) cc = lsa(F_VAL, K_VAL, DU, Dv) dd = lsa_discrete(F_VAL, K_VAL, DU, Dv) z = (mu_ - pr) / se_ if se_ > 0 else float("nan") if sfin.mean() < THR_UNIFORM: fate = "collapsed to u=1" elif np.mean([s[3] for s in sub]) > THR_RESID: fate = "coarsening, wl %.1f" % mf_ else: fate = "%s, wl %.2f" % (classify(sub[0][2], sub[0][3]), mf_) print(" %5.2f %.4f %.6f %5d %9.4f %9.4f +/- %.4f %+6.3f %+6.2f %s" % (r, Dv, cc["growth_max"], lin_ts[r], pr, mu_, se_, mu_ - pr, z, fate)) ratio_rows.append((r, Dv, cc["growth_max"], pr, dd["wl"], cc["wl_max"], mu_, se_, mf_, sf_, float(sfin.mean()), fate)) diffs = np.array([row[6] - row[3] for row in ratio_rows]) rels = np.array([100 * (row[6] - row[3]) / row[3] for row in ratio_rows]) print() print(" Across all %d ratios, club minus linear theory in the linear phase:" % len(RATIOS)) print(" mean difference %+.4f cells" % diffs.mean()) print(" mean absolute difference %.4f cells" % np.abs(diffs).mean()) print(" mean absolute difference %.2f%% of the predicted wavelength" % np.abs(rels).mean()) print(" largest single difference %.4f cells at r = %.2f" % (np.abs(diffs).max(), ratio_rows[int(np.argmax(np.abs(diffs)))][0])) print(" The z column above is not the right thing to read row by row: with") print(" only %d replicates the standard error is itself poorly determined," % REP_R) print(" so a row can show a large z purely because its three runs happened") print(" to agree. The mean absolute difference is the honest summary.") print() lr = np.log(np.array([row[0] for row in ratio_rows])) lm = np.log(np.array([row[6] for row in ratio_rows])) lp = np.log(np.array([row[3] for row in ratio_rows])) sm, im = np.polyfit(lr, lm, 1) sp, ip = np.polyfit(lr, lp, 1) rm = lm - (sm * lr + im) se_slope = math.sqrt((rm ** 2).sum() / (len(lr) - 2) / ((lr - lr.mean()) ** 2).sum()) print(" Power-law fit of the linear-phase wavelength, wavelength ~ r^p:") print(" club measurement p = %+.4f +/- %.4f" % (sm, se_slope)) print(" linear theory p = %+.4f" % sp) print(" asymptotic theory p = %+.4f (q_c^2 -> r g_v / (2 Du) for large r)" % -0.5) print(" club minus theory: %+.4f, which is %.2f standard errors." % (sm - sp, (sm - sp) / se_slope)) print() print(" The exponent is not -1/2 over this range because r is not large.") print(" The exact critical wavenumber is q_c^2 = (r g_v + f_u) / (2 Du), and") print(" with f_u = %+.5f the offset still matters at r of a few. Club and" % f_u) print(" theory bend by the same amount, which is the thing being tested.") print() print(" WHAT HAPPENS AFTERWARDS IS NOT THE SAME AT EVERY RATIO, and this") print(" surprised us. Below r of about 1.8 the pattern grows to a visible") print(" amplitude and then the whole field collapses onto the other stable") print(" state, u = 1, v = 0, leaving nothing. The Turing bifurcation here is") print(" subcritical: there is no small-amplitude patterned state to land on,") print(" so the overshoot escapes to the far attractor. Linear theory says") print(" nothing about this, and cannot. Final amplitudes:") for row in ratio_rows: print(" r = %5.2f final std(v) = %.6f %s" % (row[0], row[10], row[11])) # ------------------------------------------------------------------- head("PART 6. THE PHASE DIAGRAM") # ------------------------------------------------------------------- print(""" Now the standard setting, Du = %.2f, Dv = %.2f, r = 2, and Pearson's protocol: the trivial state u = 1, v = 0 perturbed by finite-amplitude patches, which is what actually nucleates structure in this part of parameter space. %d values of F crossed with %d values of k, %d points, each integrated %d steps and then classified automatically. Classification, decided before the sweep was run: std(v) < %.3f uniform residual drift > %.1e per step chaos (it never settles) global alignment > %.2f stripes (one direction wins) local anisotropy < %.2f spots (round objects) otherwise maze (locally striped, no direction) Local anisotropy is the gradient structure tensor of v, smoothed over about half a wavelength, reported as |lambda1 - lambda2| / (lambda1 + lambda2) and averaged with gradient-energy weighting. It is 0 for a perfectly round blob and 1 for a perfect stripe. Global alignment is the same quantity computed once for the whole field, so a labyrinth with no preferred direction scores near zero while aligned stripes score near one. """ % (DU, DV_STD, len(SWEEP_F), len(SWEEP_K), len(SWEEP_F) * len(SWEEP_K), STEPS_SWEEP, THR_UNIFORM, THR_RESID, THR_ALIGN, THR_ANISO)) sseeds = ss_sweep.spawn(len(SWEEP_F) * len(SWEEP_K)) jobs = [] ii = 0 for F in SWEEP_F: for k in SWEEP_K: jobs.append((ii, float(F), float(k), sseeds[ii])) ii += 1 res_s = run_jobs(job_sweep, jobs, "phase-diagram points") res_s.sort(key=lambda t: t[0]) n_points = len(res_s) grid = {} crops = {} print() print(" Every point, in full:") print(" F k std(v) wl aniso align fill residual class") print(" ------ ------ -------- ------- ------ ------ ----- --------- --------") for idx, F, k, mm, resid, cls, crop in res_s: grid[(round(F, 6), round(k, 6))] = (mm, resid, cls) crops[(round(F, 6), round(k, 6))] = crop pat = mm["std"] >= THR_UNIFORM print(" %.4f %.4f %.6f %7s %6s %6s %5.3f %.3e %s" % (F, k, mm["std"], ("%7.3f" % mm["wl"]) if pat else " -", ("%.4f" % mm["aniso"]) if pat else " -", ("%.4f" % mm["align"]) if pat else " -", mm["fill"], resid, cls)) symbols = {"uniform": ".", "spots": "o", "stripes": "=", "maze": "#", "chaos": "~"} print() print(" PHASE MAP (. uniform o spots # maze = stripes ~ chaos)") print(" A trailing T marks a point where the strict Turing conditions hold.") print() print(" k = " + " ".join("%5.3f" % k for k in SWEEP_K)) for F in SWEEP_F: row = [] for k in SWEEP_K: cls = grid[(round(float(F), 6), round(float(k), 6))][2] r = lsa_no_scan(float(F), float(k), DU, DV_STD, "lower") row.append(symbols[cls] + ("T" if (r and r[1]) else " ")) print(" F = %.4f %s" % (F, " ".join(row))) census = {} for row in res_s: census[row[5]] = census.get(row[5], 0) + 1 print() print(" CENSUS of %d parameter points:" % n_points) for cls in ("uniform", "spots", "maze", "stripes", "chaos"): n = census.get(cls, 0) print(" %-8s %4d %5.1f%%" % (cls, n, 100.0 * n / n_points)) print() print(" Wavelength by class, counting only points that settled:") for cls in ("spots", "maze"): w = np.array([row[3]["wl"] for row in res_s if row[5] == cls]) if len(w) > 1: print(" %-6s n = %2d mean %.3f sd %.3f range %.2f to %.2f" % (cls, len(w), w.mean(), w.std(ddof=1), w.min(), w.max())) allw = np.array([row[3]["wl"] for row in res_s if row[5] in ("spots", "maze")]) print(" all settled patterns: n = %d, mean %.3f, sd %.3f" % (len(allw), allw.mean(), allw.std(ddof=1))) print() print(" Those spreads are inflated by a handful of points where the pattern") print(" survives as a few isolated objects in an otherwise empty field: the") print(" spacing between three spots is not a wavelength. Restricting to") print(" points where the minority phase covers at least 10 percent of the area, so") print(" that the pattern genuinely fills the domain:") fillw = [(row[3]["wl"], row[5]) for row in res_s if row[5] in ("spots", "maze") and min(row[3]["fill"], 1 - row[3]["fill"]) >= 0.10] for cls in ("spots", "maze"): w = np.array([f[0] for f in fillw if f[1] == cls]) if len(w) > 1: print(" %-6s n = %2d mean %.3f sd %.3f range %.2f to %.2f" % (cls, len(w), w.mean(), w.std(ddof=1), w.min(), w.max())) fw = np.array([f[0] for f in fillw]) print(" all space-filling patterns: n = %d, mean %.3f, sd %.3f, range %.2f to %.2f" % (len(fw), fw.mean(), fw.std(ddof=1), fw.min(), fw.max())) print(" excluded as too sparse: %d of %d" % (len(allw) - len(fw), len(allw))) print(" F ranges over a factor of %.1f and k over a factor of %.2f inside" % (max(SWEEP_F) / min(SWEEP_F), max(SWEEP_K) / min(SWEEP_K))) print(" this set, and the wavelength moves by %.1f%% either side of its" % (100 * fw.std(ddof=1) / fw.mean())) print(" mean. That is what a diffusion-set length scale should do: F and") print(" k change the shape of the pattern, Du and Dv change its size.") print() fills = [(row[1], row[2], row[3]["fill"], row[3]["aniso"], row[5]) for row in res_s] sp_fill = [f[2] for f in fills if f[4] == "spots"] mz_fill = [f[2] for f in fills if f[4] == "maze"] if sp_fill: print(" spots: high-v area fraction mean %.3f (range %.3f to %.3f)" % (np.mean(sp_fill), min(sp_fill), max(sp_fill))) if mz_fill: print(" maze : high-v area fraction mean %.3f (range %.3f to %.3f)" % (np.mean(mz_fill), min(mz_fill), max(mz_fill))) print() print(" FIGURE 3 DATA: 48 x 48 crops of v, thresholded at its own midpoint,") print(" taken from the centre of the %d x %d field. '#' is high v." % (N, N)) for Fs, ks, name in SHOWCASE: Fn = min(SWEEP_F, key=lambda a: abs(a - Fs)) kn = min(SWEEP_K, key=lambda a: abs(a - ks)) key = (round(float(Fn), 6), round(float(kn), 6)) mm, resid, cls = grid[key] print() print(" panel %s: F=%.4f k=%.4f class=%s wl=%.3f aniso=%.3f fill=%.3f" % (name, Fn, kn, cls, mm["wl"], mm["aniso"], mm["fill"])) for row in crops[key]: print(" " + "".join("#" if b else "." for b in row)) # ------------------------------------------------------------------- head("PART 7. ARE ALIGNED STRIPES AN ATTRACTOR AT ALL?") # ------------------------------------------------------------------- print(""" The sweep almost never produces aligned stripes, because nothing in a periodic square started from isotropic random noise tells the pattern which way to point. What it produces instead is a labyrinth: stripes with no shared direction. To check that aligned stripes are a genuine solution of the model rather than something the model cannot do, we repeat two Turing points with the symmetry broken by hand, seeding a single one-dimensional cosine on the steady state. """) strseeds = ss_stripe.spawn(2) for i, (F, k) in enumerate([(F_VAL, K_VAL), (0.083, 0.0610)]): U, V = ic_striped(strseeds[i], F, k, mval=20) snaps = integrate(U, V, F, k, DU, DV_STD, 4000, snapshot_at=(3500, 4000)) mm = measure(snaps[4000]) resid = float(np.abs(snaps[4000] - snaps[3500]).mean()) \ / max(float(snaps[4000].std()), 1e-12) / 500.0 print(" F = %.4f, k = %.4f, seeded with a 1D cosine of wavelength %.2f" % (F, k, N / 20)) print(" after 4000 steps: wl = %.3f local anisotropy = %.3f" % (mm["wl"], mm["aniso"])) print(" GLOBAL ALIGNMENT = %.3f amplitude %.4f class = %s" % (mm["align"], mm["std"], classify(mm, resid))) print() print(" Aligned stripes are stable once they exist. They are simply never") print(" chosen out of isotropic noise, because no direction is preferred.") print(" A zebra is not a periodic square, and that difference is not") print(" cosmetic: it is the whole of why real animals have oriented stripes") print(" and our model has labyrinths.") # ------------------------------------------------------------------- head("SUMMARY OF THE NUMBERS THE ARTICLE QUOTES") # ------------------------------------------------------------------- print(" validation point F = %.4f, k = %.4f, Du/Dv = %.1f" % (F_VAL, K_VAL, DU / DV_STD)) print(" uniform steady state u* = %.6f, v* = %.6f" % (v["u"], v["v"])) print(" trace J / det J %+.6f / %+.6f" % (v["tr"], v["det"])) print(" Turing coefficient B %+.6f" % v["B"]) print(" fastest-growing wavelength, theory %.4f cells on this grid" % d["wl"]) print(" same, through our own estimator %.4f cells" % pred_lin) print(" measured wavelength, linear phase %.4f +/- %.4f cells (%+.2f%%, z = %+.2f)" % (m_lin, se_lin, 100 * (m_lin - pred_lin) / pred_lin, z_lin)) print(" measured wavelength, saturated %.4f +/- %.4f cells (%+.2f%%, z = %+.2f)" % (m_fin, se_fin, 100 * (m_fin - pred_d) / pred_d, z_fin)) print(" same on a %d x %d grid %.4f +/- %.4f cells (%+.2f%%)" % (N_BIG, N_BIG, mbf, sbf, 100 * (mbf - dbig["wl"]) / dbig["wl"])) print(" dispersion relation, worst mode %.2e relative" % worst) print(" equal diffusion, decay rate error %.2f%% worst of 9 runs" % eq_worst) print(" equal diffusion, surviving amplitude %.2e against %.4f patterned" % (eq_max_amp, amps.mean())) print(" phase diagram %d points: %s" % (n_points, ", ".join("%s %d" % (c, census.get(c, 0)) for c in ("uniform", "spots", "maze", "stripes", "chaos")))) print(" space-filling pattern wavelength %.3f +/- %.3f cells over %d points" % (fw.mean(), fw.std(ddof=1), len(fw))) print(" diffusion-ratio sweep, linear phase mean |club - theory| = %.3f cells (%.2f%%)" % (np.abs(diffs).mean(), np.abs(rels).mean())) print(" wavelength vs diffusion ratio p = %+.4f +/- %.4f (theory %+.4f)" % (sm, se_slope, sp)) print(" strict Turing region at r = 2 %.1f%% of the scanned window" % (100.0 * n_t / n_tot)) print() print(" total wall clock: %.1f s" % (time.time() - t_start)) if __name__ == "__main__": main()