import numpy as np
from PIL import Image
from pathlib import Path

SEED = 7
N = 700          # simulation grid
OUT = 1600       # render size
STEPS = 1500
DELAY = 70       # memory delay in steps
MEM = 0.006      # weight of delayed copy of earlier displacement
GAMMA = 0.004    # damping

# Seven official Flexoki swatches. RGB is assigned only after resizing fields.
H = lambda s: tuple(int(s[i:i+2], 16) for i in (0, 2, 4))
NEG = [H(x) for x in ("205EA6", "3AA99F", "87D3C3")]
POS = [H(x) for x in ("D14D41", "ECCB60", "FAEEC6")]
BLACK = H("100F0F")

rng = np.random.default_rng(SEED)
y, x = np.mgrid[0:N, 0:N].astype(np.float32) / (N - 1)
yy = 1 - y                      # up is positive
# irregular boundary: radius-modulated blob stretched along the diagonal
cx, cy = 0.5, 0.5
ang = np.arctan2(yy - cy, x - cx)
r = np.hypot(x - cx, yy - cy)
rb = 0.62 + 0.0 * ang
for k in range(2, 14):
    rb = rb + rng.normal(0, 0.09 / k ** 0.8) * np.cos(k * ang + rng.uniform(0, 6.28))
d = ((x - cx) + (yy - cy)) / 1.414
rb = rb + 0.0 * d
mask = (r < rb * (0.95 + 0.0)) & (x + (1 - yy) * 0.0 < 1.0)
mask &= (r < 0.47 + 0.07 * np.sin(5 * ang + 1.0) + 0.05 * np.cos(11 * ang))
mask = mask.astype(np.float32)
# decay gradient: wave speed falls toward upper-right so structures shrink
t = np.clip(((x - 0.0) + (yy - 0.0)) / 2, 0, 1)
c2 = (0.33 * (1 - 0.8 * t ** 1.1)) ** 2
c2 = c2.astype(np.float32)
# broad initial disturbance lower-left
u = (np.exp(-(((x - 0.28) ** 2 + (yy - 0.30) ** 2) / 0.0005)) * 1.0).astype(np.float32)
u -= 0.5 * np.exp(-(((x - 0.38) ** 2 + (yy - 0.22) ** 2) / 0.0003)).astype(np.float32)
for _ in range(14):
    bx, by = rng.uniform(0.12, 0.45), rng.uniform(0.12, 0.45)
    u += rng.normal(0, 0.35) * np.exp(-(((x - bx) ** 2 + (yy - by) ** 2) / 0.00025)).astype(np.float32)
u *= mask
up = u.copy()
hist = [u.copy()] * (DELAY + 1)
for s in range(STEPS):
    lap = (np.roll(u, 1, 0) + np.roll(u, -1, 0) + np.roll(u, 1, 1) + np.roll(u, -1, 1) - 4 * u)
    delayed = hist[-DELAY]
    un = 2 * u - up + c2 * lap * 4 - GAMMA * (u - up) + MEM * (delayed - u)
    un *= mask
    up, u = u, un
    hist.append(u.copy()); hist.pop(0)

def blur(a, sig):
    k = np.fft.fftfreq(N)
    g = np.exp(-2 * (np.pi * sig) ** 2 * (k[:, None] ** 2 + k[None, :] ** 2))
    return np.real(np.fft.ifft2(np.fft.fft2(a) * g)).astype(np.float32)
hp = (u - blur(u, 4)) * mask
env = np.sqrt(np.maximum(blur(hp * hp, 26), 0)) 
env = env + 0.25 * env.max()
u = (hp / env).astype(np.float32)
im = Image.fromarray(u.astype(np.float32), mode="F").resize((OUT, OUT), Image.BICUBIC)
f = np.asarray(im)
m = np.asarray(Image.fromarray(mask * 255).resize((OUT, OUT), Image.BICUBIC)) > 127
v = f / (np.percentile(np.abs(f[m]), 99.0) + 1e-9)
# Quiet the small upper-right focus and repetitive fringe through amplitude
# and contour thresholds, retaining the same band vocabulary everywhere.
yg, xg = np.mgrid[0:OUT, 0:OUT].astype(np.float32) / (OUT - 1)
diag = np.clip((xg + 1 - yg - 0.72) / 0.58, 0, 1)
notch = np.exp(-(((xg - 0.735) / 0.115) ** 2 + ((yg - 0.435) / 0.10) ** 2))
spacing = 0.67 + 0.21 * np.sin(18*xg + 9*yg) + 0.12 * np.cos(7*xg - 21*yg)
strength = (1 - 0.12 * diag ** 1.7) * (1 - 0.80 * notch)
v *= strength
threshold = 0.075 + 0.030 * diag * spacing
img = np.empty((OUT, OUT, 3), np.uint8)
img[:] = BLACK
a = np.abs(v)
levels = np.where(a < 0.34, 0, np.where(a < 0.64, 1, 2))
posm = (v > threshold) & m
negm = (v < -threshold) & m
img[posm] = np.array(POS, np.uint8)[levels[posm]]
img[negm] = np.array(NEG, np.uint8)[levels[negm]]
# nested cores: a pale paper swatch inside the strongest warm peaks, and a
# black hairline ring inside broad cool bands
core = posm & (a > 1.0)
img[core] = H("FAEEC6")
ring = negm & (a > 0.80) & (a < 0.84)
img[ring] = BLACK
# Cross tactic: retain the champion's broken contour fragments, but let
# selected interiors inherit the runner-up's closed contour. A short,
# southeast-shifted fragment acts as the delayed reflection of that contour.
from scipy import ndimage
labels, count = ndimage.label(ring)
print("Crossed contour components:", count)
for label_id in range(1, count + 1):
    region = labels == label_id
    py, px = np.nonzero(region)
    width, height = np.ptp(px) + 1, np.ptp(py) + 1
    if min(width, height) < 45 or len(px) / (width * height) > .4:
        continue
    center = np.array([px.mean(), py.mean()])
    coords = np.column_stack((px, py)) - center
    eigvals, eigvecs = np.linalg.eigh(np.cov(coords.T))
    transformed = coords @ eigvecs
    scale = np.sqrt(np.maximum(eigvals, 1))
    theta = np.arctan2(transformed[:, 0] / scale[0], transformed[:, 1] / scale[1])
    phase = (center[0] * .007 + center[1] * .011) % (2*np.pi)
    turn = (theta + phase) % (2*np.pi)
    # Alternate between the runner-up's closed loop and the champion's
    # broken arcs, keeping the original contour vocabulary coherent.
    crossed = (int(center[0] / 135) + int(center[1] / 115)) % 2 == 0
    if crossed:
        img[region] = np.array(NEG[2], np.uint8)
        img[region] = BLACK
        # A delayed copy of one arc, reflected down-right, stays on the
        # cool band so the large dark field and palette remain untouched.
        echo = ((turn > .42) & (turn < 1.12))
        fragments = np.zeros_like(region)
        fragments[py[echo], px[echo]] = True
        fragments = ndimage.binary_dilation(fragments, iterations=2)
        shifted = np.zeros_like(region)
        dy, dx = 18, 24
        shifted[dy:, dx:] = fragments[:-dy, :-dx]
        valid = shifted & negm & (a > .73)
        img[valid] = BLACK
    else:
        img[region] = np.array(NEG[2], np.uint8)
        fragments = np.zeros_like(region)
        keep = ((turn > .16) & (turn < 2.0)) | ((turn > 3.0) & (turn < 3.65)) | ((turn > 4.8) & (turn < 5.35))
        fragments[py[keep], px[keep]] = True
        fragments = ndimage.binary_dilation(fragments, iterations=2)
        fragments &= negm & (a > .73)
        img[fragments] = BLACK
# Repair: suppress the loose right-hand scatter and make the delayed return read as one related contour.
from PIL import ImageDraw
# Remove only small, fully detached marks in the outer-right field; leave the interior field intact.
labels, n = ndimage.label(np.any(img != np.array(BLACK, np.uint8), axis=2))
counts = np.bincount(labels.ravel())
for label_id in range(1, n + 1):
    if counts[label_id] >= 6200:
        continue
    py, px = np.where(labels == label_id)
    if len(px) and px.min() > 830 and py.min() > 250:
        img[py, px] = BLACK
# A paired, gently bending trace repeats the main outer sweep as a delayed reflection.
def bezier(points, n=300):
    p = np.asarray(points, dtype=float)
    t = np.linspace(0, 1, n)[:, None]
    q = (1-t)**3*p[0] + 3*(1-t)**2*t*p[1] + 3*(1-t)*t**2*p[2] + t**3*p[3]
    return [tuple(v) for v in q]
canvas = Image.fromarray(img, "RGB")
draw = ImageDraw.Draw(canvas)
# One long contour pair, kept beyond the main sweep and away from the quiet upper-right field.
paths = [
    [(760, 330), (1070, 440), (1000, 790), (1170, 1010)],
    [(790, 365), (1100, 490), (1040, 820), (1205, 1045)],
    [(1170, 1010), (1290, 1160), (1220, 1370), (1050, 1485)],
    [(1205, 1045), (1325, 1190), (1250, 1400), (1080, 1510)],
]
colors = [NEG[0], POS[0], NEG[0], POS[0]]
for pts, color in zip(paths, colors):
    line = bezier(pts)
    draw.line(line, fill=color, width=13, joint="curve")
    inner = [(x, y) for x, y in line]
    draw.line(inner, fill=NEG[1] if color == NEG[0] else POS[1], width=4, joint="curve")
canvas.save(Path(__file__).with_name("render-1.png"))