"""Point pattern analysis in Python: quadrat test, Clark-Evans, and Ripley's G, F, K and L functions
with Monte Carlo envelopes, on London plane trees in the London Borough of Camden.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/point-pattern-analysis-python/
Data: Trees In Camden (London Borough of Camden, OGL v3), ONS local authority boundary (OGL v3).
"""
import sys
import warnings
from importlib.metadata import version
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
import shapely
import requests
import matplotlib
matplotlib.use("Agg")   # render to files only; no GUI backend
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib.patches import Patch
from scipy import stats
from scipy.spatial import cKDTree
import pointpats
from pointpats import random as pp_random

warnings.filterwarnings("ignore")
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})

DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)
SEED = 42
NSIM = 199            # CSR simulations: a pointwise min/max envelope then has alpha = 2 / (199 + 1) = 0.01
CRS = "EPSG:27700"    # British National Grid, metres
SPECIES = "London plane"
LAD = "E09000007"     # Camden

TREES_URL = "https://opendata.camden.gov.uk/api/views/csqp-kdss/rows.csv?accessType=DOWNLOAD"
ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
LAD_URL = f"{ONS}/Local_Authority_Districts_May_2026_Boundaries_UK_BFC/FeatureServer/0/query"

BLUE, ORANGE, VERMILION, GREY = "#0072B2", "#E69F00", "#D55E00", "#9a9a9a"   # Okabe-Ito


def download(url, path, **params):
    if not path.exists():
        r = requests.get(url, params=params or None, timeout=600, headers={"User-Agent": "TechnicalGIS-tutorial"})
        r.raise_for_status()
        path.write_bytes(r.content)
    return path


def section(title):
    print(f"\n=== {title} ===")


def fmt_p(p):
    return "< 1e-300" if p < 1e-300 else f"= {p:.2g}"


section("versions")
print(f"python {sys.version.split()[0]} | pointpats {version('pointpats')} | geopandas {version('geopandas')} | "
      f"shapely {version('shapely')} | scipy {version('scipy')} | numpy {version('numpy')}")

# ---------------------------------------------------------------- Step 1: the data and the study window
section("Camden trees")
df = pd.read_csv(download(TREES_URL, DATA / "trees_in_camden.csv"), low_memory=False)
snapshot = pd.to_datetime(df["Last Uploaded"], format="%d/%m/%Y %I:%M:%S %p").max()
print(f"{len(df):,} tree records (snapshot uploaded {snapshot:%d %B %Y})")
df = df[df.Easting > 0]
print(f"{len(df):,} with coordinates")
trees = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.Easting, df.Northing), crs=CRS)

boundary = gpd.read_file(download(LAD_URL, DATA / "camden_lad.geojson", where=f"LAD26CD='{LAD}'",
                                  outFields="LAD26CD,LAD26NM", outSR=27700, f="geojson")).set_crs(CRS, allow_override=True)
W = boundary.geometry.iloc[0]
shapely.prepare(W)
trees = trees[shapely.contains_xy(W, trees.Easting, trees.Northing)]
print(f"{len(trees):,} inside the borough boundary ({W.area / 1e6:.2f} km2, perimeter {W.length / 1e3:.1f} km)")
print("by contract area:", trees["Contract Area"].value_counts().to_dict())

planes = trees[trees["Common Name"] == SPECIES]
n_raw = len(planes)
planes = planes.drop_duplicates(["Easting", "Northing"])
xy = np.column_stack([planes.Easting, planes.Northing]).astype(float)
n = len(xy)
print(f"{SPECIES}: {n_raw:,} records, {n_raw - n} at duplicate coordinates dropped -> n = {n:,}")
print(f"  by contract area: {planes['Contract Area'].value_counts().to_dict()}")

# all council-tree locations: the population that random labelling (Step 6) draws from
all_xy = np.unique(np.column_stack([trees.Easting, trees.Northing]).astype(float), axis=0)
print(f"all council trees at distinct locations: {len(all_xy):,}")

# A second, data-driven window: the parts of Camden within 100 m of any council tree. It removes
# Hampstead Heath, Regent's Park and other land where the council does not manage trees.
support = shapely.intersection(W, shapely.union_all(shapely.buffer(shapely.points(all_xy), 100, quad_segs=4)))
shapely.prepare(support)
xmin, ymin, xmax, ymax = W.bounds
bbox_area = np.ptp(xy[:, 0]) * np.ptp(xy[:, 1])
print(f"windows: borough {W.area / 1e6:.2f} km2 | within 100 m of a council tree {support.area / 1e6:.2f} km2 | "
      f"bounding box of the planes {bbox_area / 1e6:.2f} km2")
lam = n / W.area
print(f"intensity: {lam * 1e4:.2f} planes per hectare ({lam * 1e6:.0f} per km2)")

fig, axes = plt.subplots(1, 2, figsize=(11, 5.4))
for ax in axes:
    boundary.boundary.plot(ax=ax, color="black", linewidth=0.8)
    ax.set_axis_off(); ax.set_aspect("equal")
axes[0].scatter(all_xy[:, 0], all_xy[:, 1], s=0.4, color=GREY, label="other council trees", rasterized=True)
axes[0].scatter(xy[:, 0], xy[:, 1], s=1.6, color=VERMILION, label=f"{SPECIES} (n = {n:,})", rasterized=True)
axes[0].legend(loc="lower left", markerscale=6, frameon=False)
axes[0].set_title("Council-managed trees in Camden")
axes[0].plot([xmin + 200, xmin + 1200], [ymax - 300, ymax - 300], color="black", lw=2)
axes[0].text(xmin + 700, ymax - 220, "1 km", ha="center", fontsize=8)
gpd.GeoSeries([support], crs=CRS).plot(ax=axes[1], color="#cfe3f2", edgecolor=BLUE, linewidth=0.6)
x0, y0 = xy.min(axis=0); x1, y1 = xy.max(axis=0)
axes[1].plot([x0, x1, x1, x0, x0], [y0, y0, y1, y1, y0], color=ORANGE, linewidth=1.2, linestyle="--")
axes[1].scatter(xy[:, 0], xy[:, 1], s=1.2, color=VERMILION, rasterized=True)
axes[1].text(527300, 186600, "Hampstead\nHeath", ha="center", fontsize=8, color="#444444")
axes[1].text(528400, 182750, "Regent's\nPark", ha="center", fontsize=8, color="#444444")
axes[1].legend(handles=[Line2D([], [], color="black", lw=0.8, label=f"borough boundary ({W.area / 1e6:.1f} km²)"),
                        Patch(facecolor="#cfe3f2", edgecolor=BLUE,
                              label=f"within 100 m of a council tree ({support.area / 1e6:.1f} km²)"),
                        Line2D([], [], color=ORANGE, lw=1.2, ls="--",
                               label=f"bounding box of the planes ({bbox_area / 1e6:.1f} km²)")],
               loc="lower left", frameon=False, fontsize=8)
axes[1].set_title("Three possible study windows")
fig.tight_layout()
fig.savefig(OUT / "01-camden-planes-windows.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 2: quadrat counts and a chi-squared test
section("quadrat test (500 m squares)")


def quadrats(size):
    """Square cells clipped to the borough; slivers under a quarter of a full cell are dropped."""
    cells = [shapely.box(x, y, x + size, y + size) for x in np.arange(xmin, xmax, size) for y in np.arange(ymin, ymax, size)]
    g = gpd.GeoDataFrame(geometry=gpd.GeoSeries(cells, crs=CRS).intersection(W), crs=CRS)
    return g[g.area >= 0.25 * size * size].reset_index(drop=True)


pts = gpd.GeoDataFrame(geometry=gpd.points_from_xy(xy[:, 0], xy[:, 1]), crs=CRS)


def quadrat_test(grid):
    grid["count"] = gpd.sjoin(pts, grid, predicate="within").groupby("index_right").size().reindex(grid.index, fill_value=0)
    grid["expected"] = grid["count"].sum() / grid.area.sum() * grid.area   # proportional to each clipped cell's area
    chi2 = (((grid["count"] - grid["expected"]) ** 2) / grid["expected"]).sum()
    return chi2, len(grid) - 1, stats.chi2.sf(chi2, len(grid) - 1)


Q = 500
grid = quadrats(Q)
chi2, dof, p_chi2 = quadrat_test(grid)
full = grid.area > 0.999 * Q * Q
vmr = grid.loc[full, "count"].var() / grid.loc[full, "count"].mean()
print(f"{len(grid)} cells ({full.sum()} complete squares), {grid['count'].sum():,} planes counted, "
      f"{(grid['count'] == 0).sum()} empty cells")
print(f"chi2 = {chi2:,.1f}, df = {dof}, p {fmt_p(p_chi2)}")
print(f"variance/mean ratio on complete squares: {vmr:.2f} (1 under CSR, > 1 clustered, < 1 regular)")
for qs in (250, 1000):   # the verdict is stable across cell sizes, but the statistic is not
    c2, d2, p2 = quadrat_test(quadrats(qs))
    print(f"  {qs:>4} m squares: {d2 + 1} cells, chi2 = {c2:,.1f}, df = {d2}, p {fmt_p(p2)}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.8), gridspec_kw={"width_ratios": [1.15, 1]})
grid.plot(column="count", ax=axes[0], cmap="viridis", edgecolor="white", linewidth=0.3, legend=True,
          legend_kwds={"label": "London planes per 500 m cell", "shrink": 0.75})
boundary.boundary.plot(ax=axes[0], color="black", linewidth=0.6)
axes[0].set_axis_off()
axes[0].set_title(f"Quadrat counts: χ² = {chi2:,.0f} on {dof} df, p < 0.001")
counts = grid.loc[full, "count"].to_numpy()
mu = counts.mean()
edges = np.arange(0, counts.max() + 10, 5)
axes[1].hist(counts, bins=edges, color=BLUE, alpha=0.8, label=f"observed ({full.sum()} complete cells)")
expected_bins = np.diff(stats.poisson.cdf(edges - 1, mu)) * len(counts)   # P(lo <= X < hi) for each 5-count bin
axes[1].step(edges[:-1], expected_bins, where="post", color=VERMILION, lw=2, label=f"Poisson, mean {mu:.1f} (CSR)")
axes[1].set_xlabel("London planes in a 500 m cell"); axes[1].set_ylabel("number of cells")
axes[1].set_title(f"Too many empty and too many full cells (variance/mean = {vmr:.1f})")
axes[1].legend(frameon=False)
fig.tight_layout()
fig.savefig(OUT / "02-quadrat-counts.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 3: Clark-Evans nearest-neighbour index
section("Clark-Evans nearest-neighbour index")
nnd = cKDTree(xy).query(xy, k=2)[0][:, 1]
A, P = W.area, W.length
exp_naive = 0.5 / np.sqrt(n / A)
se_naive = 0.26136 / np.sqrt(n * n / A)
R = nnd.mean() / exp_naive
z = (nnd.mean() - exp_naive) / se_naive
# Donnelly (1978): expected NND and its standard error corrected for the window's perimeter
exp_don = 0.5 * np.sqrt(A / n) + (0.0514 + 0.041 / np.sqrt(n)) * P / n
se_don = np.sqrt(0.0703 * A / n ** 2 + 0.037 * P * np.sqrt(A / n ** 5))
z_don = (nnd.mean() - exp_don) / se_don
print(f"mean NND {nnd.mean():.1f} m, median {np.median(nnd):.1f} m, "
      f"{(nnd < 5).mean():.1%} under 5 m, {(nnd < 15).mean():.1%} under 15 m")
print(f"CSR expectation 0.5/sqrt(lambda) = {exp_naive:.1f} m -> R = {R:.3f}, z = {z:.1f}, "
      f"p {fmt_p(2 * stats.norm.sf(abs(z)))}")
print(f"Donnelly edge-corrected expectation {exp_don:.1f} m -> R = {nnd.mean() / exp_don:.3f}, z = {z_don:.1f}")

# CSR in the borough window with pointpats: NSIM patterns of n points each
sims = pp_random.poisson(W, size=(n, NSIM), rng=SEED)
sim_mean_nnd = np.array([cKDTree(s).query(s, k=2)[0][:, 1].mean() for s in sims])
p_ce = (1 + (sim_mean_nnd <= nnd.mean()).sum()) / (NSIM + 1)
print(f"Monte Carlo: mean NND under CSR in the borough {sim_mean_nnd.mean():.1f} m "
      f"(range {sim_mean_nnd.min():.1f}-{sim_mean_nnd.max():.1f}); one-sided p = {p_ce:.3f} (floor {1 / (NSIM + 1):.3f})")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.0))
axes[0].hist(nnd, bins=np.arange(0, 121, 2), color=BLUE, alpha=0.85)
axes[0].axvline(exp_naive, color=VERMILION, lw=2, label=f"CSR expectation {exp_naive:.0f} m")
axes[0].axvline(nnd.mean(), color="black", lw=1.5, ls="--", label=f"observed mean {nnd.mean():.0f} m")
axes[0].set_xlabel("distance to nearest other London plane (m)"); axes[0].set_ylabel("trees")
axes[0].set_title(f"Nearest-neighbour distances (n = {n:,}; {(nnd > 120).mean():.1%} beyond 120 m not shown)")
axes[0].legend(frameon=False)
axes[1].hist(sim_mean_nnd, bins=20, color=GREY, label=f"{NSIM} CSR simulations in the borough")
axes[1].axvline(nnd.mean(), color="black", lw=1.5, ls="--", label=f"observed {nnd.mean():.1f} m")
axes[1].set_xlabel("mean nearest-neighbour distance (m)"); axes[1].set_ylabel("simulations")
axes[1].set_title(f"Clark-Evans R = {R:.2f}: far more clustered than CSR")
axes[1].legend(frameon=False, loc="upper center")
fig.tight_layout()
fig.savefig(OUT / "03-nearest-neighbour.png")
plt.close(fig)


# ---------------------------------------------------------------- Step 4: G, F, K and L with border correction
def ratio_border(values, b, r):
    """Border-corrected (reduced-sample) estimate of P(value <= r) using only points at least r from the edge:
    #{value <= r <= b} / #{b >= r}."""
    ok = values <= b
    num = np.searchsorted(np.sort(values[ok]), r, "right") - np.searchsorted(np.sort(b[ok]), r, "left")
    den = len(b) - np.searchsorted(np.sort(b), r, "left")
    return num / np.maximum(den, 1)


def g_func(pts, b, r):
    return ratio_border(cKDTree(pts).query(pts, k=2)[0][:, 1], b, r)


def f_func(pts, grid_pts, grid_b, r):
    return ratio_border(cKDTree(pts).query(grid_pts, k=1)[0], grid_b, r)


def k_func(pts, b, r, area, correct=True):
    """Ripley's K. correct=True: border correction, so each point counts only at distances it can 'see' in full."""
    pairs = cKDTree(pts).query_pairs(r.max(), output_type="ndarray")
    d = np.hypot(*(pts[pairs[:, 0]] - pts[pairs[:, 1]]).T)
    lam_hat = len(pts) / area
    if not correct:
        return 2 * np.searchsorted(np.sort(d), r, "right") / (len(pts) * lam_hat)
    i = np.concatenate([pairs[:, 0], pairs[:, 1]]); dd = np.concatenate([d, d])
    bi = b[i]; ok = dd <= bi
    num = np.searchsorted(np.sort(dd[ok]), r, "right") - np.searchsorted(np.sort(bi[ok]), r, "left")
    m = len(b) - np.searchsorted(np.sort(b), r, "left")
    return num / (np.maximum(m, 1) * lam_hat)


def l_minus_r(K, r):
    return np.sqrt(K / np.pi) - r


def edge_dist(pts, window):
    return shapely.distance(shapely.points(pts), window.boundary)


section("G, F, K and L functions (border correction, CSR envelopes)")
r_g = np.arange(0, 81, 1.0)
r_f = np.arange(0, 401, 5.0)
r_k = np.arange(0, 501, 5.0)
b_obs = edge_dist(xy, W)
gx, gy = np.meshgrid(np.arange(xmin + 10, xmax, 20), np.arange(ymin + 10, ymax, 20))
gpts = np.column_stack([gx.ravel(), gy.ravel()])
gpts = gpts[shapely.contains_xy(W, gpts[:, 0], gpts[:, 1])]
gb = edge_dist(gpts, W)
print(f"{len(gpts):,} F-function sample points on a 20 m grid; "
      f"{(b_obs >= 500).sum():,} of {n:,} planes are at least 500 m from the boundary")

obs = {"G": g_func(xy, b_obs, r_g), "F": f_func(xy, gpts, gb, r_f), "K": k_func(xy, b_obs, r_k, A)}
sim = {"G": [], "F": [], "K": [], "K_raw": []}
for s in sims:
    bs = edge_dist(s, W)
    sim["G"].append(g_func(s, bs, r_g)); sim["F"].append(f_func(s, gpts, gb, r_f))
    sim["K"].append(k_func(s, bs, r_k, A)); sim["K_raw"].append(k_func(s, bs, r_k, A, correct=False))
sim = {key: np.array(v) for key, v in sim.items()}
L_obs, L_sim = l_minus_r(obs["K"], r_k), l_minus_r(sim["K"], r_k)


def report(name, r, o, s):
    above = r[o > s.max(axis=0)]; below = r[o < s.min(axis=0)]
    span = lambda a: f"{a.min():.0f}-{a.max():.0f} m" if len(a) else "none"
    print(f"{name}: above the envelope at {span(above)} ({len(above)} of {len(r)} distances), below at {span(below)}")


for name, r in (("G", r_g), ("F", r_f), ("K", r_k)):
    report(name, r, obs[name], sim[name])
for rr in (10, 50, 100, 200, 400):
    i = np.searchsorted(r_k, rr)
    print(f"  r = {rr:>3} m: L(r) - r = {L_obs[i]:6.1f} m (CSR envelope {L_sim[:, i].min():5.1f} to {L_sim[:, i].max():5.1f})")
i = np.searchsorted(r_g, 15)
print(f"G(15 m) = {obs['G'][i]:.2f} (CSR envelope {sim['G'][:, i].min():.3f}-{sim['G'][:, i].max():.3f})")
i = np.searchsorted(r_f, 100)
print(f"F(100 m) = {obs['F'][i]:.2f} (CSR envelope {sim['F'][:, i].min():.2f}-{sim['F'][:, i].max():.2f})")

fig, axes = plt.subplots(2, 2, figsize=(10.5, 8))
panels = [(r_g, obs["G"], sim["G"], "G(r): nearest-neighbour distances", "G(r)"),
          (r_f, obs["F"], sim["F"], "F(r): empty-space distances", "F(r)"),
          (r_k, obs["K"] / 1e6, sim["K"] / 1e6, "K(r): pairs within r", "K(r)  (km²)"),
          (r_k, L_obs, L_sim, "L(r) − r: variance-stabilised K", "L(r) − r  (m)")]
for ax, (r, o, s, title, ylab) in zip(axes.ravel(), panels):
    ax.fill_between(r, s.min(axis=0), s.max(axis=0), color=GREY, alpha=0.45, lw=0, label=f"CSR envelope ({NSIM} sims)")
    ax.plot(r, s.mean(axis=0), color=GREY, lw=1, ls="--", label="CSR mean")
    ax.plot(r, o, color=VERMILION, lw=2, label=f"London planes (n = {n:,})")
    ax.set_title(title); ax.set_xlabel("distance r (m)"); ax.set_ylabel(ylab)
    ax.grid(alpha=0.3)
axes[0, 0].legend(frameon=False, loc="lower right")
fig.suptitle("Border-corrected summary functions with pointwise CSR envelopes (Camden borough window)", y=1.0)
fig.tight_layout()
fig.savefig(OUT / "04-g-f-k-l-envelopes.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 5: the window and the edge correction matter
section("window shape and edge correction")
_, l_pp = pointpats.l(xy, support=r_k, linearized=True)       # pointpats: bounding-box area, no edge correction
L_raw = l_minus_r(k_func(xy, b_obs, r_k, A, correct=False), r_k)
b_sup = edge_dist(xy, support)
L_sup = l_minus_r(k_func(xy, b_sup, r_k, support.area), r_k)
m_sup = np.array([(b_sup >= rr).sum() for rr in r_k])
L_sup[m_sup < 100] = np.nan   # too few planes far enough from the holes in this window to estimate K
print(f"planes at least r from the window edge: borough {(b_obs >= 250).sum():,} at 250 m, {(b_obs >= 500).sum():,} at 500 m | "
      f"support window {(b_sup >= 250).sum():,} at 250 m, {(b_sup >= 500).sum():,} at 500 m "
      f"(curve dropped beyond {r_k[m_sup >= 100].max():.0f} m)")
L_raw_sim = l_minus_r(sim["K_raw"], r_k)
for rr in (100, 250, 500):
    i = np.searchsorted(r_k, rr)
    print(f"r = {rr} m: L - r = {l_pp[i]:6.1f} (pointpats, bounding box) | {L_raw[i]:6.1f} (borough, no correction) | "
          f"{L_obs[i]:6.1f} (borough, border) | {L_sup[i]:6.1f} (100 m support window, border)")
    print(f"           CSR mean in the borough: {L_raw_sim[:, i].mean():6.1f} without correction, "
          f"{L_sim[:, i].mean():5.1f} with border correction")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2))
axes[0].plot(r_k, l_pp, color=ORANGE, lw=2, label=f"bounding box, {bbox_area / 1e6:.1f} km² (pointpats.l)")
axes[0].plot(r_k, L_obs, color=VERMILION, lw=2, label=f"borough, {A / 1e6:.1f} km²")
axes[0].plot(r_k, L_sup, color=BLUE, lw=2, label=f"near a council tree, {support.area / 1e6:.1f} km²")
axes[0].axhline(0, color="black", lw=0.8)
axes[0].set_title("Same trees, different windows")
axes[0].set_xlabel("distance r (m)"); axes[0].set_ylabel("L(r) − r  (m)")
axes[0].legend(frameon=False, fontsize=8); axes[0].grid(alpha=0.3)
axes[1].fill_between(r_k, L_raw_sim.min(axis=0), L_raw_sim.max(axis=0), color=ORANGE, alpha=0.35, lw=0,
                     label="CSR envelope, no edge correction")
axes[1].fill_between(r_k, L_sim.min(axis=0), L_sim.max(axis=0), color=BLUE, alpha=0.35, lw=0,
                     label="CSR envelope, border correction")
axes[1].axhline(0, color="black", lw=0.8, label="theoretical CSR value")
axes[1].set_title("Ignoring edges biases the CSR benchmark downwards")
axes[1].set_xlabel("distance r (m)"); axes[1].set_ylabel("L(r) − r  (m)")
axes[1].legend(frameon=False, fontsize=8, loc="lower left"); axes[1].grid(alpha=0.3)
fig.tight_layout()
fig.savefig(OUT / "05-window-edge-effects.png")
plt.close(fig)


# ---------------------------------------------------------------- Step 6: a global test, and a better null hypothesis
def mad_test(obs_curve, sim_curves):
    """Maximum absolute deviation test (Diggle 1986; Loosmore and Ford 2006) over all distances at once."""
    t_obs = np.abs(obs_curve - sim_curves.mean(axis=0)).max()
    # each simulation is compared with the mean of the others, so it is treated exactly like the data
    total = sim_curves.sum(axis=0)
    t_sim = np.array([np.abs(s - (total - s) / (len(sim_curves) - 1)).max() for s in sim_curves])
    return t_obs, t_sim, (1 + (t_sim >= t_obs).sum()) / (len(t_sim) + 1)


section("global MAD test against CSR")
t_obs, t_sim, p_mad = mad_test(L_obs, L_sim)
print(f"MAD on L(r) - r, 0-500 m: observed {t_obs:.1f} m; CSR simulations {np.median(t_sim):.1f} m median, "
      f"{t_sim.max():.1f} m max; p = {p_mad:.3f}")
# how often does a pointwise envelope 'reject' for a pattern that really is CSR? (each simulation against the rest)
outside = [((L_sim[j] > np.delete(L_sim, j, 0).max(axis=0)) | (L_sim[j] < np.delete(L_sim, j, 0).min(axis=0))).any()
           for j in range(NSIM)]
print(f"pointwise envelope check: {np.mean(outside):.0%} of CSR patterns leave the envelope of the other "
      f"{NSIM - 1} at one or more of {len(r_k)} distances")

section("random labelling: are planes more clustered than council trees in general?")
rng = np.random.default_rng(SEED)
b_all = edge_dist(all_xy, W)
L_rl = []
for _ in range(NSIM):
    idx = rng.choice(len(all_xy), size=n, replace=False)
    L_rl.append(l_minus_r(k_func(all_xy[idx], b_all[idx], r_k, A), r_k))
L_rl = np.array(L_rl)
L_all = l_minus_r(k_func(all_xy, b_all, r_k, A), r_k)
report("L - r vs random labelling", r_k, L_obs, L_rl)
for rr in (20, 50, 100, 200, 400):
    i = np.searchsorted(r_k, rr)
    print(f"  r = {rr:>3} m: planes {L_obs[i]:6.1f} m | all council trees {L_all[i]:6.1f} m | "
          f"random-label envelope {L_rl[:, i].min():6.1f} to {L_rl[:, i].max():6.1f}")
t_obs_rl, t_sim_rl, p_mad_rl = mad_test(L_obs, L_rl)
print(f"MAD test vs random labelling: observed {t_obs_rl:.1f} m; simulations median {np.median(t_sim_rl):.1f} m, "
      f"max {t_sim_rl.max():.1f} m; p = {p_mad_rl:.3f}")
for lo, hi in ((0, 100), (100, 500)):
    sel = (r_k >= lo) & (r_k <= hi)
    _, _, p_band = mad_test(L_obs[sel], L_rl[:, sel])
    print(f"  restricted to {lo}-{hi} m: p = {p_band:.3f}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2))
axes[0].fill_between(r_k, L_sim.min(axis=0), L_sim.max(axis=0), color=GREY, alpha=0.45, lw=0, label="CSR envelope")
axes[0].fill_between(r_k, L_rl.min(axis=0), L_rl.max(axis=0), color=BLUE, alpha=0.35, lw=0,
                     label=f"random labelling envelope ({NSIM} draws)")
axes[0].plot(r_k, L_all, color=BLUE, lw=1.2, ls="--", label=f"all council trees (n = {len(all_xy):,})")
axes[0].plot(r_k, L_obs, color=VERMILION, lw=2, label="London planes")
axes[0].set_xlabel("distance r (m)"); axes[0].set_ylabel("L(r) − r  (m)")
axes[0].set_title("Planes against two null hypotheses")
axes[0].legend(frameon=False, fontsize=8, loc="lower right", bbox_to_anchor=(1, 0.13)); axes[0].grid(alpha=0.3)
bins = np.linspace(0, max(t_obs, t_obs_rl, t_sim.max(), t_sim_rl.max()) * 1.05, 50)
axes[1].hist(t_sim, bins=bins, color=GREY, alpha=0.8, label="CSR simulations")
axes[1].hist(t_sim_rl, bins=bins, color=BLUE, alpha=0.6, label="random-labelling simulations")
axes[1].axvline(t_obs, color="black", lw=1.5, ls="--", label=f"planes vs CSR: {t_obs:.0f} m (p = {p_mad:.3f})")
axes[1].axvline(t_obs_rl, color=VERMILION, lw=2, label=f"planes vs random labelling: {t_obs_rl:.0f} m (p = {p_mad_rl:.3f})")
axes[1].set_xlabel("maximum absolute deviation of L(r) − r over 0–500 m (m)"); axes[1].set_ylabel("simulations")
axes[1].set_title(f"Global MAD test (p = {1 / (NSIM + 1):.3f} is the smallest possible)")
axes[1].legend(frameon=False, fontsize=8)
fig.tight_layout()
fig.savefig(OUT / "06-random-labelling-mad.png")
plt.close(fig)

planes[["Identifier", "Common Name", "Contract Area", "Height In Metres", "Spread In Metres", "geometry"]].to_file(
    OUT / "camden_london_planes.gpkg")
gpd.GeoDataFrame({"window": ["borough", "support_100m"]}, geometry=[W, support], crs=CRS).to_file(
    OUT / "camden_windows.gpkg")
print(f"\nwrote {OUT / 'camden_london_planes.gpkg'} and {OUT / 'camden_windows.gpkg'}")