#!/usr/bin/env python3
"""Stubborn residue: a crossed deposit mat, then a lift that fails where claims stack."""
import math
import os
import sys
from collections import deque
import numpy as np
from PIL import Image, ImageDraw

N = 1200
BLACK = (16, 15, 15)      # #100F0F
CYAN = (58, 169, 159)     # cyan-400 #3AA99F
YELLOW = (208, 162, 21)   # yellow-400 #D0A215
ORANGE = (218, 112, 44)   # orange-400 #DA702C
RED = (209, 77, 65)       # red-400 #D14D41
PAL = [BLACK, CYAN, YELLOW, ORANGE, RED]


def hash2(x, y, s):
    n = (int(x) * 374761393 + int(y) * 668265263 + int(s) * 1442695041) & 0xFFFFFFFF
    n = (n ^ (n >> 13)) * 1274126177 & 0xFFFFFFFF
    return ((n ^ (n >> 16)) & 0xFFFFFFFF) / 4294967295.0


def flow_angle(x, y, seed, phase):
    ix = int(x) // 48
    iy = int(y) // 48
    a = hash2(ix, iy, seed + phase)
    return a * math.tau * 0.35 + phase  # keep a family near its bias, with local drift


def draw_family(dr, rng, seed, phase, gap, steps, width):
    for gy in range(0, N, gap):
        for gx in range(0, N, gap):
            x = gx + rng.random() * gap
            y = gy + rng.random() * gap
            ang = flow_angle(x, y, seed, phase) + (rng.random() - 0.5) * 0.55
            pts = []
            for t in range(steps):
                pts.append((x, y))
                if t % 6 == 0:
                    ang = 0.75 * ang + 0.25 * flow_angle(x, y, seed, phase)
                ang += (rng.random() - 0.5) * 0.4
                x += math.cos(ang) * 2.0
                y += math.sin(ang) * 2.0
                if x < -20 or y < -20 or x > N + 20 or y > N + 20:
                    break
            if len(pts) > 2:
                dr.line(pts, fill=1, width=width)


def deposit(seed, gap, steps):
    claim = np.zeros((N, N), dtype=np.uint16)
    col = np.zeros((N, N), dtype=np.uint8)
    families = [
        (0.15, 1, 2),
        (1.65, 2, 2),
        (2.9, 3, 2),
    ]
    for phase, ci, width in families:
        rng = np.random.default_rng(seed * 10 + ci * 17)
        mask = Image.new("L", (N, N), 0)
        dr = ImageDraw.Draw(mask)
        draw_family(dr, rng, seed, phase, gap, steps, width)
        layer = np.asarray(mask, dtype=np.uint16)
        claim += layer
        col[layer > 0] = ci
    return np.minimum(claim, 18).astype(np.uint8), col


def reinforce(claim, col, seed, origin):
    rng = np.random.default_rng(seed + 404)
    mask = Image.new("L", (N, N), 0)
    dr = ImageDraw.Draw(mask)
    x = float(origin[0]) + rng.uniform(-40, 40)
    y = float(origin[1]) + rng.uniform(-40, 40)
    ang = rng.random() * math.tau
    for i in range(1100):
        ang += (rng.random() - 0.5) * 0.5
        x += math.cos(ang) * 5.5
        y += math.sin(ang) * 5.5
        if x < 20 or y < 20 or x > N - 20 or y > N - 20:
            ang += math.pi * 0.65
            x = min(max(x, 20), N - 21)
            y = min(max(y, 20), N - 21)
        pts = []
        a = ang + (rng.random() - 0.5) * 1.0
        px, py = x, y
        for t in range(12):
            pts.append((px, py))
            a += (rng.random() - 0.5) * 0.5
            px += math.cos(a) * 1.7
            py += math.sin(a) * 1.7
        dr.line(pts, fill=1, width=3)
    layer = np.asarray(mask, dtype=np.uint16)
    claim = claim.astype(np.uint16) + layer * 2
    col = col.copy()
    col[layer > 0] = 4
    return np.minimum(claim, 18).astype(np.uint8), col


def erase(claim, seed, origin, p_table, max_frac, wander):
    rng = np.random.default_rng(seed + 97)
    lifted = np.zeros((N, N), dtype=np.bool_)
    seen = np.zeros((N, N), dtype=np.bool_)
    ox, oy = origin
    ang = rng.random() * math.tau
    x, y = float(ox), float(oy)
    q = deque()
    for _ in range(wander):
        x += math.cos(ang) * 2.0
        y += math.sin(ang) * 2.0
        ang += (rng.random() - 0.5) * 1.1
        xi, yi = int(x) % N, int(y) % N
        if not seen[yi, xi]:
            seen[yi, xi] = True
            q.append((yi, xi))
    target = int(N * N * max_frac)
    n = 0
    dirs = [(-1, 0), (1, 0), (0, -1), (0, 1), (-1, -1), (1, 1), (-1, 1), (1, -1)]
    while q and n < target:
        yi, xi = q.popleft()
        c = int(claim[yi, xi])
        p = p_table[c] if c < len(p_table) else p_table[-1]
        if rng.random() < p:
            lifted[yi, xi] = True
            n += 1
            if n % 29 == 0:
                rng.shuffle(dirs)
            for dy, dx in dirs:
                ny, nx_ = yi + dy, xi + dx
                if 0 <= ny < N and 0 <= nx_ < N and not seen[ny, nx_]:
                    seen[ny, nx_] = True
                    q.append((ny, nx_))
    return lifted


def paint(claim, col, lifted):
    img = np.empty((N, N, 3), dtype=np.uint8)
    img[:] = BLACK
    remain = (~lifted) & (claim > 0)
    for k in range(1, len(PAL)):
        img[remain & (col == k)] = PAL[k]
    return img


def mono(claim, lifted):
    img = np.full((N, N), 16, dtype=np.uint8)
    remain = (~lifted) & (claim > 0)
    img[remain] = 80
    img[remain & (claim >= 2)] = 140
    img[remain & (claim >= 4)] = 200
    img[remain & (claim >= 6)] = 245
    return img


def measure(claim, lifted):
    area = int(lifted.sum())
    ys, xs = np.where(lifted)
    cy = float(ys.mean()) / N if len(ys) else 0.5
    cx = float(xs.mean()) / N if len(xs) else 0.5
    # boundary roughness: std of radius from centroid
    if len(ys) > 100:
        # sample edge
        edge = lifted.copy()
        edge[1:-1, 1:-1] &= ~(
            lifted[:-2, 1:-1] & lifted[2:, 1:-1] & lifted[1:-1, :-2] & lifted[1:-1, 2:]
        )
        ey, ex = np.where(edge)
        if len(ey) > 50:
            rad = np.hypot(ex - cx * N, ey - cy * N)
            rough = float(rad.std() / (rad.mean() + 1))
        else:
            rough = 0
    else:
        rough = 0
    dil = lifted.copy()
    for _ in range(5):
        n = dil.copy()
        n[1:] |= dil[:-1]
        n[:-1] |= dil[1:]
        n[:, 1:] |= dil[:, :-1]
        n[:, :-1] |= dil[:, 1:]
        dil = n
    islands = dil & ~lifted & (claim > 0)
    return {
        "cover": round(float((claim > 0).mean()), 3),
        "void": round(area / (N * N), 3),
        "remain": round(float(((~lifted) & (claim > 0)).mean()), 3),
        "islands_frac": round(float(islands.mean()), 3),
        "cx": round(cx, 3),
        "cy": round(cy, 3),
        "off": round(math.hypot(cx - 0.5, cy - 0.5), 3),
        "rough": round(rough, 3),
    }


def pick_origin(claim, seed):
    rng = np.random.default_rng(seed + 3)
    best = None
    for _ in range(120):
        x = int(rng.integers(80, N - 80))
        y = int(rng.integers(80, N - 80))
        if (x - N / 2) ** 2 + (y - N / 2) ** 2 < (N * 0.15) ** 2:
            continue
        local = float(claim[y - 14:y + 14, x - 14:x + 14].mean())
        score = abs(local - 1.5)
        if best is None or score < best[0]:
            best = (score, x, y, local)
    return best[1], best[2], round(best[3], 2)


def run(seed, out_png=None, out_mono=None, gap=14, steps=26, max_frac=0.20):
    claim, col = deposit(seed, gap, steps)
    ox, oy, local = pick_origin(claim, seed)
    claim, col = reinforce(claim, col, seed, (ox, oy))
    p_table = [0.99, 0.62, 0.22, 0.08, 0.03, 0.01, 0.0]
    lifted = erase(claim, seed, (ox, oy), p_table, max_frac, wander=380)
    # fray the contour: singly claimed pixels on the edge may still lift; stacked claims hold
    rng = np.random.default_rng(seed + 515)
    edge = np.zeros_like(lifted)
    edge[1:-1, 1:-1] = (~lifted[1:-1, 1:-1]) & (
        lifted[:-2, 1:-1] | lifted[2:, 1:-1] | lifted[1:-1, :-2] | lifted[1:-1, 2:]
    )
    ys, xs = np.where(edge & (claim <= 2))
    for y, x in zip(ys.tolist(), xs.tolist()):
        c = int(claim[y, x])
        p = 0.62 if c <= 1 else 0.28
        if rng.random() < p:
            lifted[y, x] = True
    st = measure(claim, lifted)
    st.update(origin=(ox, oy), local=local, seed=seed)
    if out_png:
        Image.fromarray(paint(claim, col, lifted), "RGB").save(out_png, optimize=True)
    if out_mono:
        Image.fromarray(mono(claim, lifted), "L").save(out_mono, optimize=True)
    return st


if __name__ == "__main__":
    here = os.path.dirname(os.path.abspath(__file__))
    out = os.path.join(here, "render-1.png")
    st = run(4104, out, None)
    print(st)