This post is not about tiling, but about the periodic behavior of the three colors ("species"), which appears to settle down into the same oscillation throughout, but phase-shifted. Once you have identified the phase at each position, you can compute its entire periodic history. I imagine this is well known to those who study this, but it was news to me, and inspired me to use this fact to creating a looping animation, which I'll post next assuming animated GIFs work here.
Since Claude generated the code and the empirical results, I had it produce a detailed write-up. The rest of this post is from AI, but it was the result of a long discussion with a lot of trial and error.
Here is the complete code to reproduce the images below.
Code: Select all
import argparse
import numpy as np
from scipy.ndimage import convolve
from scipy.optimize import curve_fit, brentq
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from PIL import Image
# Standalone BZ reaction-diffusion simulator with PERIODIC boundary
# conditions and a SQUARE (Cartesian) neighborhood, instead of the rhombus
# tile + hexagonal-footprint + rotation/reflection-equivalence machinery
# used in zb3_padded.py. That machinery existed to solve a real but
# separate problem (correct wraparound for a non-rectangular rhombus
# tile); on a plain square grid, scipy.ndimage.convolve's built-in
# mode='wrap' already gives exact periodic boundary conditions for free,
# so none of that is needed here.
#
# Workflow:
# 1. Run the simulation for --iterations steps (random initial
# condition), discarding --skip as transient, on an n x n periodic
# grid with a (2*radius+1) x (2*radius+1) square averaging
# neighborhood.
# 2. Measure the oscillation period directly from the CENTER pixel's
# species-A peak spacing (this grid has translational symmetry, so
# "the center" is not special dynamically, but it's a convenient,
# unambiguous reference point).
# 3. Fit the single periodic pulse-shape model (see bz_empirical_fit.py)
# to the center pixel's trajectory, exactly as before.
# 4. For every pixel in the grid, infer its PHASE OFFSET relative to the
# center pixel by numerically inverting the fitted curve against
# that pixel's own (most recent) measured value -- not by re-fitting
# shape parameters per pixel, since the locality findings (see
# bz_locality.py) showed neighboring/all pixels share the same
# attractor shape and differ only in phase, once alpha=beta=gamma.
# 5. Render frames at ANY requested time purely by evaluating the
# fitted formula at (t + pixel_phase) for every pixel -- no further
# simulation steps needed.
parser = argparse.ArgumentParser(
description="Run a periodic-boundary, square-neighborhood BZ "
"simulation, infer per-pixel phase from a single fitted "
"waveform, and render frames at arbitrary requested times "
"without further simulation.")
parser.add_argument('-p', '--param', type=float, default=1.0,
help='Value for alpha=beta=gamma (default: 1.0).')
parser.add_argument('-n', '--size', type=int, default=200,
help='Side length of the (square, periodic) simulation grid '
'(default: 200).')
parser.add_argument('-r', '--radius', type=int, default=3,
help='Neighborhood radius; the averaging filter is a '
'(2*radius+1) x (2*radius+1) square block (default: 3).')
parser.add_argument('-i', '--iterations', type=int, default=700,
help='Total iterations to run (default: 700).')
parser.add_argument('--skip', type=int, default=450,
help='Iterations discarded as transient before period-fitting/phase '
'inference (default: 450).')
parser.add_argument('--seed', type=int, default=42,
help='Random seed for the initial condition (default: 42).')
parser.add_argument('--render-times', type=float, nargs='+', default=None,
help='One or more continuous-time values (relative to the END of the '
'simulated/skip window) at which to render a static frame. Default: '
'render one frame per third of the fitted period (0, period/3, '
'2*period/3), to show all three species taking a turn at being '
'dominant. Ignored if --gif is given.')
parser.add_argument('--gif', action='store_true',
help='Instead of a few static frames, render --gif-frames evenly spaced '
'frames across exactly one period and assemble them into a seamless '
'looping animated GIF.')
parser.add_argument('--gif-frames', type=int, default=60,
help='Number of frames per loop for --gif (default: 60).')
parser.add_argument('--gif-fps', type=float, default=20.0,
help='Playback frames-per-second for --gif (default: 20.0).')
parser.add_argument('--gif-size', type=int, default=400,
help='Output GIF width/height in pixels (default: 400; the simulation '
'grid is resized to this for the GIF, independent of --size).')
parser.add_argument('--phase-blur-sigma', type=float, default=1.5,
help='Gaussian blur sigma (in pixels) applied to the inferred phase map, '
'done correctly for a circular/periodic quantity via cos/sin '
'components. Set to 0 to disable. Smooths out pixel-level phase-'
'inference noise; near genuine spiral-core defects, blurring cannot '
'recover a well-defined phase (none exists there), but it keeps those '
'pixels from looking like sharp, arbitrary noise (default: 1.5).')
parser.add_argument('-o', '--output-prefix', default='bz_phase',
help='Filename prefix for output images (default: bz_phase). Produces '
'<prefix>_fit.png, <prefix>_phasemap.png, and <prefix>_frame_N.png '
'for each render time.')
args = parser.parse_args()
n = args.size
r = args.radius
niter = args.iterations
skip = args.skip
alpha = beta = gamma = args.param
np.random.seed(args.seed)
# Square (2r+1)x(2r+1) neighborhood, uniformly weighted, normalized -- the
# direct Cartesian analog of the original's truncated-corner hex filter,
# but with no corners cut, since there's no hex-tiling constraint here.
filt_size = 2 * r + 1
FILTER = np.ones((filt_size, filt_size)) / (filt_size * filt_size)
# Random initial condition, n x n, no padding needed -- mode='wrap' in
# convolve handles periodic boundaries directly.
arr = np.zeros((2, 3, n, n))
arr[0] = np.random.random(size=(3, n, n))
center = n // 2
trajectory = np.zeros((niter - skip, 3)) # center pixel's A,B,C
growth = np.zeros((niter - skip, 3)) # center pixel's g_a,g_b,g_c
full_at_skip_end = None # full-grid snapshot at the last iteration
for i in range(niter):
p = i % 2
q = (p + 1) % 2
s = np.zeros((3, n, n))
for k in range(3):
s[k] = convolve(arr[p, k], FILTER, mode='wrap')
g_a = alpha * s[1] - gamma * s[2]
g_b = beta * s[2] - alpha * s[0]
g_c = gamma * s[0] - beta * s[1]
arr[q, 0] = s[0] + s[0] * g_a
arr[q, 1] = s[1] + s[1] * g_b
arr[q, 2] = s[2] + s[2] * g_c
np.clip(arr[q], 0, 1, arr[q])
if i >= skip:
trajectory[i - skip] = arr[q, :, center, center]
growth[i - skip] = (g_a[center, center], g_b[center, center], g_c[center, center])
if i == niter - 1:
full_at_skip_end = arr[q].copy() # (3, n, n), the final settled state
print(f"Simulation complete: n={n}, radius={r} (filter {filt_size}x{filt_size}), "
f"param={args.param}, {niter} iterations ({skip} discarded as transient).")
# --- Measure period from the center pixel's peak spacing. ---
y_a = trajectory[:, 0]
peaks = [i for i in range(1, len(y_a) - 1)
if y_a[i] > y_a[i - 1] and y_a[i] > y_a[i + 1] and y_a[i] > 0.5]
if len(peaks) < 3:
raise SystemExit("Not enough peaks found at center pixel to measure period -- "
"try more iterations or a different seed.")
period_est = float(np.mean(np.diff(peaks)))
print(f"Measured period at center pixel: {period_est:.4f} iterations "
f"({len(peaks)} peaks found)")
def pulse(phase, c_rise, k_rise, c_fall, k_fall, period):
"""Smooth, genuinely periodic plateau-rise-plateau-fall pulse (see
bz_empirical_fit.py for the derivation of why sin() is used here
instead of a raw mod-wrapped linear argument -- this avoids a kink
that the simpler version develops when a center lands near the wrap)."""
omega = 2 * np.pi / period
rise = 1.0 / (1.0 + np.exp(-k_rise * period / (2 * np.pi) * np.sin(omega * (phase - c_rise))))
fall = 1.0 / (1.0 + np.exp(k_fall * period / (2 * np.pi) * np.sin(omega * (phase - c_fall))))
return rise * fall
# --- Fit the pulse shape to the center pixel's species-A trace. ---
search_len = int(round(period_est)) + 1
trough_idx = int(np.argmin(y_a[:search_len]))
start_offset = (trough_idx - int(round(period_est * 0.1))) % len(y_a)
n_periods_to_fit = max(1, min(4, (len(y_a) - start_offset) // max(1, int(period_est))))
fit_len = int(round(period_est * n_periods_to_fit))
fit_len = min(fit_len, len(y_a) - start_offset)
y_fit = y_a[start_offset:start_offset + fit_len]
t_fit_local = np.arange(fit_len)
first_period = y_fit[:int(round(period_est))]
rise_guess = next((i for i in range(1, len(first_period))
if first_period[i - 1] < 0.5 <= first_period[i]), period_est * 0.25)
fall_guess = next((i for i in range(int(rise_guess) + 1, len(first_period))
if first_period[i - 1] >= 0.5 > first_period[i]), period_est * 0.75)
def model(t, c_rise, k_rise, c_fall, k_fall):
return pulse(t, c_rise, k_rise, c_fall, k_fall, period_est)
p0 = [rise_guess, 1.0, fall_guess, 1.0]
bounds = ([0, 0.05, 0, 0.05], [period_est, 10, period_est, 10])
popt, _ = curve_fit(model, t_fit_local, y_fit, p0=p0, bounds=bounds, maxfev=20000)
c_rise_global = float(np.mod(popt[0] + start_offset, period_est))
c_fall_global = float(np.mod(popt[2] + start_offset, period_est))
shape_params = np.array([c_rise_global, popt[1], c_fall_global, popt[3]])
y_pred = model(t_fit_local, *popt)
ss_res = np.sum((y_fit - y_pred) ** 2)
ss_tot = np.sum((y_fit - y_fit.mean()) ** 2)
r2 = 1 - ss_res / ss_tot if ss_tot > 0 else float('nan')
print(f"Fitted shape: c_rise={shape_params[0]:.3f}, k_rise={shape_params[1]:.3f}, "
f"c_fall={shape_params[2]:.3f}, k_fall={shape_params[3]:.3f} (R^2={r2:.4f})")
shift = period_est / 3.0
def species_curve(phase, species_shift):
"""The fitted curve, evaluated for species A (shift=0), B (shift=shift),
or C (shift=2*shift)."""
return pulse(phase + species_shift, *shape_params, period_est)
# --- Fit-quality plot at the center pixel. ---
names = ['A', 'B', 'C']
colors = ['red', 'green', 'blue']
fig_fit, axes_fit = plt.subplots(3, 1, figsize=(11, 8), sharex=True)
t_check = np.arange(fit_len)
for k in range(3):
y_k = trajectory[start_offset:start_offset + fit_len, k]
sp_shift = k * shift
y_pred_k = species_curve(t_check + start_offset, sp_shift)
ss_res_k = np.sum((y_k - y_pred_k) ** 2)
ss_tot_k = np.sum((y_k - y_k.mean()) ** 2)
r2_k = 1 - ss_res_k / ss_tot_k if ss_tot_k > 0 else float('nan')
print(f" {names[k]} (shift={sp_shift:.3f}): R^2 vs actual = {r2_k:.4f}")
ax = axes_fit[k]
ax.plot(t_check, y_k, color=colors[k], alpha=0.5, linewidth=1.2, label='simulation')
ax.plot(t_check, y_pred_k, color='black', linewidth=1.3, linestyle='--', label='fitted model')
ax.set_ylabel(names[k])
ax.set_ylim(-0.05, 1.05)
ax.legend(loc='upper right', fontsize=8)
axes_fit[-1].set_xlabel('iteration (within fit window)')
fig_fit.suptitle(f'Center-pixel fit, n={n}, radius={r}, param={args.param}, period~{period_est:.2f}')
fig_fit.tight_layout()
fig_fit.savefig(f'{args.output_prefix}_fit.png', dpi=120)
print(f"Saved {args.output_prefix}_fit.png")
# --- Turnaround vs growth-rate zero-crossing alignment plot. Confirms
# that each species' peak/trough coincides with its own growth-rate term
# (g_a, g_b, g_c) crossing zero -- i.e. the direction of each species'
# rise/decay is set by the OTHER two species, not by itself. ---
def find_turnarounds(y):
return [i for i in range(1, len(y) - 1)
if (y[i] > y[i - 1] and y[i] > y[i + 1]) or (y[i] < y[i - 1] and y[i] < y[i + 1])]
def find_zero_crossings(g):
sign = np.sign(g)
return [i for i in range(len(g) - 1)
if sign[i] != 0 and sign[i + 1] != 0 and sign[i] != sign[i + 1]]
fig_turn, axes_turn = plt.subplots(3, 1, figsize=(12, 9), sharex=True)
t_full = np.arange(skip, niter)
for k in range(3):
y_k = trajectory[:, k]
g_k = growth[:, k]
turns = find_turnarounds(y_k)
crossings = find_zero_crossings(g_k)
if turns and crossings:
offsets = np.array([min(crossings, key=lambda c: abs(c - t)) - t for t in turns])
print(f" {names[k]}: {len(turns)} turnarounds, offset vs growth-rate zero-crossing: "
f"mean={offsets.mean():.2f}, max|offset|={np.abs(offsets).max()}, "
f"median={np.median(offsets):.1f}")
ax = axes_turn[k]
ax.plot(t_full, y_k, color=colors[k], label=f'{names[k]}(t)', linewidth=1.3)
ax.set_ylabel(names[k], color=colors[k])
ax.set_ylim(-0.05, 1.05)
ax2 = ax.twinx()
ax2.plot(t_full, g_k, color='black', alpha=0.4, linewidth=1, linestyle='--',
label=f'growth rate g_{names[k]}')
ax2.axhline(0, color='gray', linewidth=0.6)
ax2.set_ylabel(f'g_{names[k]}', color='black')
for ti in turns:
ax.axvline(skip + ti, color=colors[k], alpha=0.25, linewidth=0.8)
for ci in crossings:
ax2.axvline(skip + ci, color='black', alpha=0.25, linewidth=0.8, linestyle=':')
ax.set_title(f'{names[k]}(t) [solid, colored axis] vs growth rate g_{names[k]} '
f'[dashed, black axis]')
axes_turn[-1].set_xlabel('iteration')
fig_turn.suptitle(f'Turnaround vs growth-rate zero-crossing alignment, center pixel, '
f'param={args.param}, n={n}, radius={r}')
fig_turn.tight_layout()
fig_turn.savefig(f'{args.output_prefix}_turnaround.png', dpi=120)
print(f"Saved {args.output_prefix}_turnaround.png")
# --- Infer phase at every pixel by numerically inverting the fitted curve
# against that pixel's own final (A, B, C) snapshot. ---
#
# Strategy: search over phase in [0, period) for the value that makes the
# predicted (A, B, C) triple closest (least squares over all 3 species) to
# the pixel's actual final snapshot, using a coarse grid scan followed by
# a local refinement (brentq on the derivative is overkill given the
# function is cheap to evaluate; a fine-enough grid scan is simpler and
# robust to the function being non-monotonic over a full period).
N_PHASE_SAMPLES = 720 # 0.5-iteration resolution at period~17; adjust if needed
phase_candidates = np.linspace(0, period_est, N_PHASE_SAMPLES, endpoint=False)
pred_a = species_curve(phase_candidates, 0.0)
pred_b = species_curve(phase_candidates, shift)
pred_c = species_curve(phase_candidates, 2 * shift)
pred_table = np.stack([pred_a, pred_b, pred_c], axis=1) # (N_PHASE_SAMPLES, 3)
# Flatten the grid's final snapshot to (3, n*n), compare against the table.
flat_vals = full_at_skip_end.reshape(3, n * n) # (3, n*n)
# For each pixel, find the phase_candidate index minimizing squared error
# across all 3 species. Vectorized via broadcasting: (N_PHASE_SAMPLES, 3, 1)
# vs (1, 3, n*n) -> (N_PHASE_SAMPLES, n*n) squared-error sums.
diff = pred_table[:, :, None] - flat_vals[None, :, :] # (N_PHASE_SAMPLES, 3, n*n)
sq_err = np.sum(diff ** 2, axis=1) # (N_PHASE_SAMPLES, n*n)
best_idx = np.argmin(sq_err, axis=0) # (n*n,)
phase_map_raw = phase_candidates[best_idx].reshape(n, n)
residual_map = np.sqrt(np.min(sq_err, axis=0)).reshape(n, n)
print(f"\nPhase inference: mean residual (RMS error in A,B,C space) = {residual_map.mean():.4f}, "
f"max = {residual_map.max():.4f}")
print(f"(For reference, species values range over [0,1]; a residual well under ~0.05-0.1 "
f"indicates the single-shape-fits-all-pixels assumption is holding up well.)")
# --- Smooth the phase map with a Gaussian blur, done CORRECTLY for a
# circular/periodic quantity. Phase wraps at period_est -- naively
# averaging raw phase values would be wrong right where it matters most
# (e.g. averaging phase~0.3 and phase~17.0 at period~17.38 should give
# ~17.15 i.e. near the wrap, not ~8.65, the wrong-way "average" you'd get
# from blurring the raw numbers). The standard correct fix: represent
# phase as a point on the unit circle (cos, sin), blur those two
# ordinary (non-circular) channels with a normal Gaussian filter, then
# recover phase via atan2. Right at a spiral core -- a genuine
# topological defect where neighboring phases point in every direction at
# once -- the blurred (cos, sin) vector naturally shrinks toward zero
# length, which is a meaningful signal (low confidence) rather than an
# error: there ISN'T a well-defined phase there for any method to recover.
from scipy.ndimage import gaussian_filter
omega_blur = 2 * np.pi / period_est
cos_map = np.cos(omega_blur * phase_map_raw)
sin_map = np.sin(omega_blur * phase_map_raw)
# mode='wrap' matches the simulation's own periodic boundary conditions.
cos_blurred = gaussian_filter(cos_map, sigma=args.phase_blur_sigma, mode='wrap')
sin_blurred = gaussian_filter(sin_map, sigma=args.phase_blur_sigma, mode='wrap')
confidence_map = np.sqrt(cos_blurred ** 2 + sin_blurred ** 2) # 1 = confident, ->0 = defect
phase_map_smoothed = np.mod(np.arctan2(sin_blurred, cos_blurred) / omega_blur, period_est)
if args.phase_blur_sigma > 0:
phase_map = phase_map_smoothed
print(f"Applied circular Gaussian blur to phase map, sigma={args.phase_blur_sigma} pixels. "
f"Confidence (post-blur vector length) mean={confidence_map.mean():.3f}, "
f"min={confidence_map.min():.3f} (near 0 = spiral core / genuine defect, not fixable "
f"by smoothing -- this is where phase truly isn't well-defined).")
else:
phase_map = phase_map_raw
print("--phase-blur-sigma=0: using raw (unsmoothed) phase map.")
# --- Plot the inferred phase map, residual map, and (if blurred) the
# confidence map. ---
n_panels = 3 if args.phase_blur_sigma > 0 else 2
fig_map, axes_map = plt.subplots(1, n_panels, figsize=(6.2 * n_panels, 5.5))
im0 = axes_map[0].imshow(phase_map, cmap='twilight', vmin=0, vmax=period_est)
axes_map[0].set_title('Inferred phase per pixel' +
(' (blurred)' if args.phase_blur_sigma > 0 else ''))
fig_map.colorbar(im0, ax=axes_map[0], shrink=0.8, label='phase (iterations)')
im1 = axes_map[1].imshow(residual_map, cmap='viridis')
axes_map[1].set_title('Phase-fit residual (RMS error, pre-blur)')
fig_map.colorbar(im1, ax=axes_map[1], shrink=0.8, label='RMS error')
if args.phase_blur_sigma > 0:
im2 = axes_map[2].imshow(confidence_map, cmap='magma', vmin=0, vmax=1)
axes_map[2].set_title('Post-blur confidence\n(low = spiral core / defect)')
fig_map.colorbar(im2, ax=axes_map[2], shrink=0.8, label='vector length')
fig_map.suptitle(f'n={n}, radius={r}, param={args.param}')
fig_map.tight_layout()
fig_map.savefig(f'{args.output_prefix}_phasemap.png', dpi=120)
print(f"Saved {args.output_prefix}_phasemap.png")
# --- Render frames, either as a few static PNGs or as a full animated
# GIF, purely from the phase map + fitted formula, no further simulation. ---
if args.gif:
n_frames = args.gif_frames
frame_times = np.linspace(0, period_est, n_frames, endpoint=False)
rgb_frames = []
for t_render in frame_times:
frame = np.zeros((n, n, 3))
for k, sp_shift in enumerate([0.0, shift, 2 * shift]):
frame[:, :, k] = species_curve(phase_map + t_render, sp_shift)
frame_u8 = np.clip(frame * 255, 0, 255).astype(np.uint8)
rgb_frames.append(frame_u8)
# Build ONE shared color palette from a sample of pixels across ALL
# frames combined, then quantize every frame against that same
# palette. Quantizing each frame independently (PIL's default when
# you just .save() a list of 'RGB' images as a GIF) lets the palette
# drift slightly frame to frame, which can cause a subtle flicker/
# color-shift even though the underlying RGB data loops perfectly --
# and it also produces a much larger file, since GIF compresses
# better when consecutive frames share a palette.
stacked = np.concatenate([f.reshape(-1, 3) for f in rgb_frames], axis=0)
palette_source = Image.fromarray(
stacked.reshape(1, -1, 3) if stacked.shape[0] < 1_000_000 else stacked[:1_000_000].reshape(1, -1, 3),
mode='RGB')
palette_img = palette_source.convert('P', palette=Image.ADAPTIVE, colors=256)
shared_palette = palette_img.getpalette()
pil_frames = []
for frame_u8 in rgb_frames:
img = Image.fromarray(frame_u8, mode='RGB')
if args.gif_size != n:
img = img.resize((args.gif_size, args.gif_size), Image.LANCZOS)
img_p = img.quantize(palette=palette_img, dither=Image.NONE)
pil_frames.append(img_p)
gif_path = f'{args.output_prefix}_loop.gif'
duration_ms = int(round(1000.0 / args.gif_fps))
pil_frames[0].save(
gif_path, save_all=True, append_images=pil_frames[1:],
duration=duration_ms, loop=0, optimize=True, disposal=2)
print(f"\nSaved animated GIF: {gif_path} ({n_frames} frames, {duration_ms}ms/frame, "
f"{args.gif_size}x{args.gif_size}px, shared 256-color palette, period={period_est:.3f} "
f"iterations covered exactly once per loop -- seamless by construction since frame N "
f"is sampled 1/{n_frames} of a period before frame 0 wraps, the same spacing as any "
f"other consecutive frame pair).")
else:
render_times = args.render_times if args.render_times is not None else \
[0.0, period_est / 3.0, 2.0 * period_est / 3.0]
for t_render in render_times:
frame = np.zeros((n, n, 3))
for k, sp_shift in enumerate([0.0, shift, 2 * shift]):
frame[:, :, k] = species_curve(phase_map + t_render, sp_shift)
fig_frame, ax_frame = plt.subplots(figsize=(6.5, 6.5))
ax_frame.imshow(frame)
ax_frame.set_title(f't = {t_render:.2f}')
ax_frame.axis('off')
fig_frame.tight_layout()
fname = f'{args.output_prefix}_frame_{t_render:.2f}.png'
fig_frame.savefig(fname, dpi=120)
plt.close(fig_frame)
print(f"Saved {fname}")
print(f"\nDone. Model summary: shape params {shape_params}, period={period_est:.4f}, "
f"shift={shift:.4f}. Any frame can now be rendered as "
f"species_curve(phase_map + t, k*shift) for k in (0,1,2), for any t, with no "
f"further simulation.")
This is a writeup of some analysis on a Belousov-Zhabotinsky-style reaction-diffusion simulation (three species A/B/C, rock-paper-scissors growth/decay, local averaging via convolution, hard-clipped to [0,1], periodic boundary conditions). Nearly every pixel in the grid runs the exact same oscillation, just phase-shifted, and that oscillation itself collapses to a small empirical formula. This is enough to render or animate the whole system without simulating it frame by frame.
All three images below, and the numbers quoted alongside them, come from a single run of the same script:
Code: Select all
python3 bz_phase_render.py -n 200 -p 1.0 -r 3 -i 700 --skip 450 --seed 42 --phase-blur-sigma 1.5 -o bz_phase
The system
Each cell holds three concentrations A, B, C -- referred to here as "species," following the usual chemical-kinetics term for the reacting quantities in this kind of system. In the rendered images, A, B, and C are mapped directly to the R, G, B channels, so they are also literally the three colors. Every iteration:
Code: Select all
s = locally_averaged(A, B, C) # convolution with a square neighborhood, periodic boundaries
A_new = s_A + s_A * (alpha * s_B - gamma * s_C)
B_new = s_B + s_B * (beta * s_C - alpha * s_A)
C_new = s_C + s_C * (gamma * s_A - beta * s_B)
clip(A_new, B_new, C_new, 0, 1)
Finding 1: every pixel runs the same cycle, just phase-shifted
After the initial transient, nearly all pixels lock into the same periodic oscillation -- same waveform shape, same period -- differing only in phase. This was checked directly: fit a single empirical model to one pixel, use it (only shifted in phase) to predict every other pixel's value at every point in time, and compare against the actual simulated data. The fit holds at R^2 > 0.99 almost everywhere on the grid (the exceptions are discussed under Finding 5).
Finding 2: A, B, and C are themselves the same function, phase-shifted by 1/3 period
Because alpha=beta=gamma, the cycle has an exact 3-fold symmetry: B(t) = A(t - period/3), C(t) = A(t - 2*period/3). This was checked directly rather than assumed -- fit a curve to species A only, then predict B and C purely by shifting that same curve, and compare against the actual simulated B/C data:
Code: Select all
A: fit directly R^2 = 0.997
B: predicted as A shifted by period/3 (5.79 iter) R^2 = 0.997
C: predicted as A shifted by 2*period/3 (11.59 iter) R^2 = 0.997
Finding 3: the mechanism behind the curve's shape
The time series look like sharp narrow spikes sitting on long near-zero plateaus, not smooth sinusoids. Two things explain this.
First, summing the three reaction terms above shows A+B+C is exactly conserved by the reaction step -- the cross terms cancel algebraically, and this checks out numerically to float64 precision in the actual simulation. Clipping is the only thing that breaks this conservation, and clipping never fully stops, even deep in the "locked in" periodic regime -- it is a continuous low-level leak, not just a transient cleanup. This is part of why the curve has the shape it does rather than being a clean sinusoid.
Second, and more directly: near a species' minimum, its own update is approximately multiplicative:
A_new ~ s_A * (1 + alpha*s_B - gamma*s_C)
-- A's value times a growth factor set entirely by the other two species. Fitting log(A) vs. iteration on the rising/falling segments gives R^2 ~ 0.97-0.98, versus R^2 ~ 0.75-0.82 for a plain linear fit on the same segments -- so each segment is exponential growth/decay, not linear.
That growth factor's sign, g_A = alpha*s_B - gamma*s_C, decides whether A is currently growing or shrinking, and has nothing to do with A's own value. The natural check is whether A's turnaround (peak or trough) lines up with g_A crossing zero. The same run above also tracks g_A, g_B, g_C at the center pixel and reports this directly:
Code: Select all
A: 29 turnarounds, offset vs growth-rate zero-crossing: mean=0.45, max|offset|=1, median=0.0
B: 29 turnarounds, offset vs growth-rate zero-crossing: mean=0.62, max|offset|=1, median=1.0
C: 28 turnarounds, offset vs growth-rate zero-crossing: mean=0.43, max|offset|=1, median=0.0
- bz_phase_turnaround.png (417.71 KiB) Viewed 176 times
This is the actual mechanism: each species keeps accelerating in whichever direction it is currently going right up until the other two swap dominance, then flips abruptly -- which is why the peaks are sharp rather than gently rounded.
Finding 4: it all reduces to 7 numbers
Given the above, the whole oscillation -- at every pixel -- can be written as one smooth, exactly periodic function:
Code: Select all
def pulse(phase, c_rise, k_rise, c_fall, k_fall, period):
omega = 2 * np.pi / period
rise = 1 / (1 + np.exp(-k_rise * period / (2*np.pi) * np.sin(omega * (phase - c_rise))))
fall = 1 / (1 + np.exp( k_fall * period / (2*np.pi) * np.sin(omega * (phase - c_fall))))
return rise * fall
Fit once to one pixel's species-A trace (same run as above), the result is 7 numbers:
Code: Select all
c_rise=7.367, k_rise=1.150, c_fall=13.500, k_fall=1.169, period=17.385, shift=period/3=5.795
- bz_phase_fit.png (155.12 KiB) Viewed 176 times
Since every pixel runs the same curve, the fitted formula can be inverted against each pixel's final simulated (A,B,C) snapshot to back out a phase value per pixel -- a full phase map of the grid, built without any further simulation:
Code: Select all
N_PHASE_SAMPLES = 720
phase_candidates = np.linspace(0, period_est, N_PHASE_SAMPLES, endpoint=False)
pred_a = pulse(phase_candidates, *shape_params, period_est)
pred_b = pulse(phase_candidates + shift, *shape_params, period_est)
pred_c = pulse(phase_candidates + 2*shift, *shape_params, period_est)
pred_table = np.stack([pred_a, pred_b, pred_c], axis=1) # (720, 3)
flat_vals = final_snapshot.reshape(3, n*n) # (3, n*n)
diff = pred_table[:, :, None] - flat_vals[None, :, :] # (720, 3, n*n)
sq_err = np.sum(diff**2, axis=1) # (720, n*n)
phase_map = phase_candidates[np.argmin(sq_err, axis=0)].reshape(n, n)
residual_map = np.sqrt(np.min(sq_err, axis=0)).reshape(n, n)
Smoothing the phase map must be done carefully, since phase is circular -- naively averaging raw phase values near a wraparound point gives the wrong answer. The correct approach is to convert phase to (cos, sin) of the phase angle, Gaussian-blur those two ordinary (non-circular) components, then convert back via atan2:
Code: Select all
from scipy.ndimage import gaussian_filter
omega = 2 * np.pi / period_est
cos_map = np.cos(omega * phase_map)
sin_map = np.sin(omega * phase_map)
cos_blur = gaussian_filter(cos_map, sigma=1.5, mode='wrap')
sin_blur = gaussian_filter(sin_map, sigma=1.5, mode='wrap')
confidence_map = np.sqrt(cos_blur**2 + sin_blur**2) # 1 = confident, -> 0 = defect
phase_map_smoothed = np.mod(np.arctan2(sin_blur, cos_blur) / omega, period_est)
- bz_phase_phasemap.png (803.98 KiB) Viewed 176 times
Application: seamless loop animation, no further simulation
Once a per-pixel phase map and the single fitted curve are available, a frame at any time t is just:
Code: Select all
frame = np.zeros((n, n, 3))
for k, sp_shift in enumerate([0.0, shift, 2 * shift]):
frame[:, :, k] = pulse(phase_map + t + sp_shift, *shape_params, period_est)
Same script, GIF mode:
Code: Select all
python3 bz_phase_render.py -n 200 -p 1.0 -r 3 -i 700 --skip 450 --seed 42 \
--phase-blur-sigma 1.5 --gif --gif-frames 60 --gif-fps 20 --gif-size 400 -o bz_phase
Code: Select all
from PIL import Image
n_frames = 60
# n_frames evenly spaced samples across exactly one period:
frame_times = np.linspace(0, period_est, n_frames, endpoint=False)
rgb_frames = []
for t_render in frame_times:
frame = np.zeros((n, n, 3))
for k, sp_shift in enumerate([0.0, shift, 2 * shift]):
frame[:, :, k] = pulse(phase_map + t_render + sp_shift, *shape_params, period_est)
frame_u8 = np.clip(frame * 255, 0, 255).astype(np.uint8)
rgb_frames.append(frame_u8)
# Build one shared color palette from all frames combined, then quantize
# every frame against that same palette. Quantizing each frame
# independently lets the palette drift slightly from frame to frame,
# which can cause a subtle flicker even though the underlying RGB data
# loops perfectly, and it also produces a larger file, since GIF
# compresses better when consecutive frames share a palette.
stacked = np.concatenate([f.reshape(-1, 3) for f in rgb_frames], axis=0)
palette_source = Image.fromarray(stacked.reshape(1, -1, 3), mode='RGB')
palette_img = palette_source.convert('P', palette=Image.ADAPTIVE, colors=256)
pil_frames = []
for frame_u8 in rgb_frames:
img = Image.fromarray(frame_u8, mode='RGB')
img_p = img.quantize(palette=palette_img, dither=Image.NONE)
pil_frames.append(img_p)
duration_ms = int(round(1000.0 / 20.0)) # 20 fps
pil_frames[0].save(
'bz_loop.gif', save_all=True, append_images=pil_frames[1:],
duration=duration_ms, loop=0, optimize=True, disposal=2)