"""DBSCAN and HDBSCAN spatial clustering in Python: finding Glasgow's dining districts
in Overture Maps places (restaurants, cafés, bars and takeaways).

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/dbscan-spatial-clustering-python/
Data: Overture Maps Foundation places (CDLA-Permissive-2.0), ONS local authority boundaries (OGL v3).
"""
import json
import sys
import time
import warnings
from importlib.metadata import version
from pathlib import Path

import duckdb
import geopandas as gpd
import hdbscan
import matplotlib
import numpy as np
import pandas as pd
import requests
import shapely

matplotlib.use("Agg")   # render to files only; no GUI backend
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib.lines import Line2D
from matplotlib.patches import Patch
from sklearn.cluster import DBSCAN, HDBSCAN
from sklearn.metrics import adjusted_rand_score, silhouette_score
from sklearn.neighbors import BallTree, NearestNeighbors

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
CRS = "EPSG:27700"            # British National Grid, metres
R_EARTH = 6_371_008.8         # mean Earth radius in metres (IUGG)
LAD = "S12000049"             # Glasgow City
RELEASE = "2026-10-21.0"      # Overture keeps only recent releases; the script falls back to the latest one
FOOD = ("restaurant", "casual_eatery", "cafe", "coffee_shop", "fast_food_restaurant", "bar")
ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
LAD_URL = f"{ONS}/Local_Authority_Districts_May_2026_Boundaries_UK_BGC/FeatureServer/0/query"

# colour-blind-safe qualitative palette (Paul Tol "muted" + "bright"), cycled for cluster ids
TOL = ["#332288", "#88CCEE", "#44AA99", "#117733", "#999933", "#DDCC77", "#CC6677", "#882255", "#AA4499",
       "#4477AA", "#EE6677", "#228833", "#CCBB44", "#66CCEE", "#AA3377"]
NOISE = "#c8c8c8"


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


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 cluster_colours(labels):
    return np.array([NOISE if l < 0 else TOL[l % len(TOL)] for l in labels])


section("versions")
print(f"python {sys.version.split()[0]} | scikit-learn {version('scikit-learn')} | hdbscan {version('hdbscan')} | "
      f"geopandas {version('geopandas')} | shapely {version('shapely')} | duckdb {version('duckdb')} | "
      f"numpy {version('numpy')}")

# ---------------------------------------------------------------- 1. data: Glasgow boundary + Overture places
section("Glasgow City boundary")
boundary = gpd.read_file(download(LAD_URL, DATA / "glasgow_lad.geojson", where=f"LAD26CD='{LAD}'",
                                  outFields="LAD26CD,LAD26NM", outSR=27700, f="geojson"))
boundary = boundary.set_crs(CRS, allow_override=True)
city = boundary.geometry.iloc[0]
xmin, ymin, xmax, ymax = boundary.to_crs(4326).total_bounds.round(4)
print(f"{boundary.LAD26NM.iloc[0]}: {city.area / 1e6:.1f} km2, lon/lat bbox {xmin}, {ymin}, {xmax}, {ymax}")

section("Overture Maps places")
cache = DATA / "glasgow_places.parquet"
if not cache.exists():
    stac = requests.get("https://stac.overturemaps.org/catalog.json", timeout=60).json()
    releases = [l["href"].split("/")[1] for l in stac["links"] if l.get("rel") == "child"]
    release = RELEASE if RELEASE in releases else stac["latest"]
    con = duckdb.connect()
    con.sql("INSTALL spatial; LOAD spatial; INSTALL httpfs; LOAD httpfs; SET s3_region = 'us-west-2'")
    files = con.sql(f"""
        SELECT assets.aws.alternate.s3.href AS url FROM
        read_parquet('https://stac.overturemaps.org/{release}/collections.parquet')
        WHERE collection = 'place' AND bbox.xmin <= {xmax} AND bbox.xmax >= {xmin}
          AND bbox.ymin <= {ymax} AND bbox.ymax >= {ymin}""").df().url.tolist()
    con.sql(f"""
        COPY (SELECT id, names.primary AS name, basic_category, taxonomy.hierarchy[1] AS top_category,
                     confidence, operating_status, addresses[1].postcode AS postcode,
                     ST_X(geometry) AS lon, ST_Y(geometry) AS lat, '{release}' AS release
              FROM read_parquet({files})
              WHERE bbox.xmin BETWEEN {xmin} AND {xmax} AND bbox.ymin BETWEEN {ymin} AND {ymax})
        TO '{cache.as_posix()}' (FORMAT parquet)""")
places = pd.read_parquet(cache)
print(f"Overture release {places.release.iloc[0]}: {len(places):,} places of all kinds in the bounding box")

g = gpd.GeoDataFrame(places, geometry=gpd.points_from_xy(places.lon, places.lat), crs=4326).to_crs(CRS)
g = g[g.within(city)].copy()
food = g[g.basic_category.isin(FOOD) & (g.confidence >= 0.5) &
         (g.operating_status.fillna("open") != "permanently_closed")].reset_index(drop=True)
print(f"inside Glasgow City: {len(g):,} places; food and drink ({', '.join(FOOD)}), "
      f"confidence >= 0.5, not permanently closed: {len(food):,}")
print(food.basic_category.value_counts().to_string())
dup = food.duplicated(subset=["lon", "lat"]).sum()
print(f"points sharing exact coordinates with another: {dup} ({dup / len(food):.1%})")
food["district"] = food.postcode.str.upper().str.extract(r"^([A-Z]{1,2}\d{1,2}[A-Z]?)\s")[0]
print(f"postcode district available for {food.district.notna().mean():.0%} of places")

X_m = np.column_stack([food.geometry.x, food.geometry.y])          # metres, British National Grid
X_rad = np.radians(food[["lat", "lon"]].to_numpy())                 # haversine wants [lat, lon] in radians

# ---------------------------------------------------------------- 2. choose eps with a k-distance plot
section("k-distance plot")
MIN_SAMPLES = 10
nn = NearestNeighbors(n_neighbors=MIN_SAMPLES).fit(X_m)     # the query point itself is neighbour 0
kdist = np.sort(nn.kneighbors(X_m)[0][:, -1])


def knee(y):
    """Kneedle-style knee: the point of a sorted curve furthest below the chord from first to last value."""
    x = np.linspace(0, 1, len(y)); yn = (y - y[0]) / (y[-1] - y[0])
    return int(np.argmax(x - yn))


# the top 2% are isolated outliers several km from anything; they would stretch the chord
k_i = knee(kdist[: int(0.98 * len(kdist))])
print(f"distance to the {MIN_SAMPLES - 1}th nearest other place: median {np.median(kdist):.0f} m, "
      f"90th percentile {np.percentile(kdist, 90):.0f} m, max {kdist.max():,.0f} m")
print(f"knee of the k-distance curve at rank {k_i} ({k_i / len(kdist):.0%}): {kdist[k_i]:.0f} m")
EPS = 150.0
print(f"chosen eps = {EPS:.0f} m ({(kdist <= EPS).mean():.0%} of places have {MIN_SAMPLES - 1} others within it)")

fig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 4.8), gridspec_kw={"width_ratios": [1.15, 1]})
boundary.boundary.plot(ax=a1, color="#555555", lw=0.6)
cols = {"restaurant": "#CC6677", "casual_eatery": "#DDCC77", "cafe": "#117733", "coffee_shop": "#44AA99",
        "fast_food_restaurant": "#882255", "bar": "#332288"}
for c in FOOD:
    s = food[food.basic_category == c]
    a1.scatter(s.geometry.x, s.geometry.y, s=2.5, color=cols[c], label=f"{c.replace('_', ' ')} ({len(s):,})", lw=0)
a1.legend(loc="lower left", fontsize=7, markerscale=3, frameon=True, title="Overture basic_category", title_fontsize=7)
a1.set_title(f"{len(food):,} food and drink places in Glasgow City")
a1.set_axis_off(); a1.set_aspect("equal")
a2.plot(np.arange(len(kdist)), kdist, color="#332288", lw=1.5)
a2.axhline(EPS, color="#CC6677", ls="--", lw=1)
a2.annotate(f"eps = {EPS:.0f} m", (0.02 * len(kdist), EPS), xytext=(0, 5), textcoords="offset points", color="#CC6677")
a2.scatter([k_i], [kdist[k_i]], color="#117733", zorder=3, s=25)
a2.annotate(f"knee ≈ {kdist[k_i]:.0f} m", (k_i, kdist[k_i]), xytext=(-90, 25), textcoords="offset points",
            color="#117733", arrowprops={"arrowstyle": "->", "color": "#117733"})
a2.set_yscale("log"); a2.set_ylim(5, kdist.max() * 1.2)
a2.set_xlabel("places sorted by k-distance"); a2.set_ylabel(f"distance to {MIN_SAMPLES - 1}th nearest neighbour (m, log)")
a2.set_title(f"k-distance plot (min_samples = {MIN_SAMPLES})"); a2.grid(alpha=0.3)
fig.savefig(OUT / "01-places-k-distance.png"); plt.close(fig)

# ---------------------------------------------------------------- 3. DBSCAN: haversine on radians vs metres
section("DBSCAN: haversine on lat/lon vs Euclidean on metres")
t0 = time.perf_counter()
db_hav = DBSCAN(eps=EPS / R_EARTH, min_samples=MIN_SAMPLES, metric="haversine", algorithm="ball_tree").fit(X_rad)
t_hav = time.perf_counter() - t0
t0 = time.perf_counter()
db = DBSCAN(eps=EPS, min_samples=MIN_SAMPLES).fit(X_m)
t_m = time.perf_counter() - t0
lab = db.labels_
for name, l, t in [("haversine (radians)", db_hav.labels_, t_hav), ("projected metres", lab, t_m)]:
    print(f"{name:20s}: {l.max() + 1} clusters, noise {np.mean(l < 0):.1%}, {t * 1000:.0f} ms")
print(f"adjusted Rand index between the two: {adjusted_rand_score(db_hav.labels_, lab):.4f}; "
      f"points labelled differently as noise/non-noise: {np.sum((db_hav.labels_ < 0) != (lab < 0))}")
bad = DBSCAN(eps=EPS, min_samples=MIN_SAMPLES, metric="haversine").fit(X_rad).labels_
print(f"mistake check: eps={EPS:.0f} passed straight to haversine (radians) -> {bad.max() + 1} cluster(s), "
      f"noise {np.mean(bad < 0):.0%}")
deg = DBSCAN(eps=EPS / 111_320, min_samples=MIN_SAMPLES).fit(food[["lon", "lat"]].to_numpy()).labels_
print(f"mistake check: Euclidean on raw degrees (eps = {EPS:.0f}/111,320) -> {deg.max() + 1} clusters, "
      f"ARI vs metres {adjusted_rand_score(deg, lab):.3f}")
core = np.zeros(len(lab), bool); core[db.core_sample_indices_] = True
print(f"core points {core.sum():,}, border points {np.sum(~core & (lab >= 0)):,}, noise {np.sum(lab < 0):,}")
sizes = pd.Series(lab[lab >= 0]).value_counts()
print(f"largest cluster {sizes.iloc[0]:,} places ({sizes.iloc[0] / len(lab):.0%} of all); median cluster size "
      f"{int(sizes.median())}")

# ---------------------------------------------------------------- 4. sensitivity to eps and min_samples
section("DBSCAN sensitivity")
EPS_GRID = [75, 100, 150, 200, 300, 400]
MS_GRID = [5, 10, 15, 20, 30]
ncl = np.zeros((len(MS_GRID), len(EPS_GRID)), int); noise = np.zeros_like(ncl, float)
big = np.zeros_like(ncl, float)
for i, ms in enumerate(MS_GRID):
    for j, e in enumerate(EPS_GRID):
        l = DBSCAN(eps=e, min_samples=ms).fit(X_m).labels_
        ncl[i, j] = l.max() + 1; noise[i, j] = np.mean(l < 0)
        big[i, j] = pd.Series(l[l >= 0]).value_counts().max() / len(l) if (l >= 0).any() else 0
print("rows min_samples, columns eps (m): number of clusters / % noise / largest cluster as % of places")
print(pd.DataFrame([[f"{ncl[i, j]:>3} / {noise[i, j]:4.0%} / {big[i, j]:4.0%}" for j in range(len(EPS_GRID))]
                    for i in range(len(MS_GRID))], index=MS_GRID, columns=EPS_GRID).to_string())

fig, axes = plt.subplots(1, 2, figsize=(11, 3.8))
for ax, M, title, cmap, fmt in [(axes[0], ncl, "Number of clusters", "viridis", "{:d}"),
                                (axes[1], noise * 100, "Noise (% of places)", "cividis", "{:.0f}")]:
    im = ax.imshow(M, cmap=cmap, aspect="auto")
    for i in range(M.shape[0]):
        for j in range(M.shape[1]):
            v = M[i, j]
            ax.text(j, i, fmt.format(int(v)) if fmt == "{:d}" else fmt.format(v), ha="center", va="center",
                    color="white" if v < (M.max() + M.min()) / 2 else "black", fontsize=8)
    ax.set_xticks(range(len(EPS_GRID)), EPS_GRID); ax.set_yticks(range(len(MS_GRID)), MS_GRID)
    ax.set_xlabel("eps (m)"); ax.set_ylabel("min_samples"); ax.set_title(title)
    jj, ii = EPS_GRID.index(int(EPS)), MS_GRID.index(MIN_SAMPLES)
    ax.add_patch(plt.Rectangle((jj - 0.5, ii - 0.5), 1, 1, fill=False, ec="#CC6677", lw=2.5))
    fig.colorbar(im, ax=ax, shrink=0.85, label="clusters" if M is ncl else "% of places")
fig.suptitle("DBSCAN on Glasgow food and drink places: sensitivity to eps and min_samples (chosen: red box)", y=1.02)
fig.savefig(OUT / "03-dbscan-sensitivity.png"); plt.close(fig)

# ---------------------------------------------------------------- 5. HDBSCAN for variable density
section("HDBSCAN")
MIN_CLUSTER = 15
t0 = time.perf_counter()
hdb = HDBSCAN(min_cluster_size=MIN_CLUSTER, min_samples=MIN_SAMPLES, cluster_selection_method="eom",
              copy=True).fit(X_m)
t_hdb = time.perf_counter() - t0
hl = hdb.labels_
hs = pd.Series(hl[hl >= 0]).value_counts()
print(f"sklearn HDBSCAN(min_cluster_size={MIN_CLUSTER}, min_samples={MIN_SAMPLES}): {hl.max() + 1} clusters, "
      f"noise {np.mean(hl < 0):.1%}, largest {hs.iloc[0]} places, median size {int(hs.median())}, "
      f"{t_hdb * 1000:.0f} ms")
hdb_hav = HDBSCAN(min_cluster_size=MIN_CLUSTER, min_samples=MIN_SAMPLES, metric="haversine").fit(X_rad)
print(f"same on lat/lon radians with metric='haversine': {hdb_hav.labels_.max() + 1} clusters, "
      f"ARI vs metres {adjusted_rand_score(hdb_hav.labels_, hl):.4f}")
leaf = HDBSCAN(min_cluster_size=MIN_CLUSTER, min_samples=MIN_SAMPLES, cluster_selection_method="leaf").fit(X_m).labels_
print(f"cluster_selection_method='leaf': {leaf.max() + 1} clusters, noise {np.mean(leaf < 0):.1%}")
for mcs in (10, 25, 50):
    l = HDBSCAN(min_cluster_size=mcs, min_samples=MIN_SAMPLES).fit(X_m).labels_
    print(f"min_cluster_size={mcs:>2}: {l.max() + 1} clusters, noise {np.mean(l < 0):.1%}")
print(f"ARI DBSCAN vs HDBSCAN: {adjusted_rand_score(lab, hl):.3f}")

# per-cluster density in places per km2 of convex hull: shows why one eps cannot fit every cluster
food["hdb"] = hl
dens = []
for c, grp in food[food.hdb >= 0].groupby("hdb"):
    a = shapely.convex_hull(shapely.MultiPoint(grp.geometry.tolist())).area / 1e6
    dens.append(len(grp) / max(a, 1e-3))
dens = np.array(dens)
print(f"HDBSCAN cluster density (places per km2 of convex hull): min {dens.min():.0f}, median {np.median(dens):.0f}, "
      f"max {dens.max():.0f} ({dens.max() / dens.min():.0f}x range)")

# the hdbscan package gives the condensed tree and DBCV
hp = hdbscan.HDBSCAN(min_cluster_size=MIN_CLUSTER, min_samples=MIN_SAMPLES, gen_min_span_tree=True).fit(X_m)
print(f"hdbscan package: {hp.labels_.max() + 1} clusters, ARI vs sklearn {adjusted_rand_score(hp.labels_, hl):.4f}, "
      f"relative_validity_ (fast DBCV approximation) {hp.relative_validity_:.3f}")
print(f"membership strength: median {np.median(hp.probabilities_[hp.labels_ >= 0]):.2f}; outlier scores > 0.9: "
      f"{np.sum(hp.outlier_scores_ > 0.9)}")

centre = (258_900, 665_300)   # George Square
zoom = (centre[0] - 4500, centre[0] + 4500, centre[1] - 3000, centre[1] + 3200)
fig, axes = plt.subplots(1, 2, figsize=(12, 4.6), layout="constrained")
for ax, l, title in [(axes[0], lab, f"DBSCAN eps={EPS:.0f} m, min_samples={MIN_SAMPLES}: {lab.max() + 1} clusters, "
                                    f"{np.mean(lab < 0):.0%} noise"),
                     (axes[1], hl, f"HDBSCAN min_cluster_size={MIN_CLUSTER}: {hl.max() + 1} clusters, "
                                   f"{np.mean(hl < 0):.0%} noise")]:
    boundary.boundary.plot(ax=ax, color="#555555", lw=0.6)
    n = l < 0
    ax.scatter(X_m[n, 0], X_m[n, 1], s=3, color=NOISE, lw=0, label="noise")
    ax.scatter(X_m[~n, 0], X_m[~n, 1], s=5, c=cluster_colours(l[~n]), lw=0)
    ax.set_xlim(zoom[:2]); ax.set_ylim(zoom[2:]); ax.set_aspect("equal")
    ax.set_xticks([]); ax.set_yticks([]); ax.set_title(title, fontsize=8.5)
    ax.plot([zoom[0] + 300, zoom[0] + 1300], [zoom[2] + 300] * 2, color="black", lw=2)
    ax.text(zoom[0] + 800, zoom[2] + 420, "1 km", ha="center", fontsize=7)
    ax.legend(handles=[Line2D([], [], marker="o", ls="", color=NOISE, label="noise"),
                       Line2D([], [], marker="o", ls="", color=TOL[0], label="cluster (colour = id)")],
              loc="upper right", fontsize=7, frameon=True)
fig.suptitle("Central Glasgow: one global eps chains the centre into a single large cluster; "
             "HDBSCAN splits it by local density")
fig.savefig(OUT / "02-dbscan-vs-hdbscan.png"); plt.close(fig)

# the minimum density DBSCAN demands of a core point: min_samples places in a circle of radius eps
db_floor = MIN_SAMPLES / (np.pi * (EPS / 1000) ** 2)
print(f"DBSCAN density floor: {MIN_SAMPLES} places in a {EPS:.0f} m circle = {db_floor:.0f} places per km2; "
      f"{np.sum(dens < db_floor)} of {len(dens)} HDBSCAN clusters are sparser than that")
fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 4.4), gridspec_kw={"width_ratios": [1.3, 1]})
# select_clusters=True draws ellipses that fail with matplotlib 3.11, so the tree is drawn plain
hp.condensed_tree_.plot(select_clusters=False, axis=a1, colorbar=True, cmap="viridis", log_size=True)
a1.set_title("HDBSCAN condensed tree (λ = 1 / distance in 1/m; width = log number of places)")
order = np.sort(dens)[::-1]
a2.bar(np.arange(len(order)), order, color="#44AA99", width=0.85)
a2.axhline(db_floor, color="#CC6677", ls="--", lw=1.2)
a2.text(0.97, 0.97, f"dashed line: DBSCAN core density floor\n{MIN_SAMPLES} places within {EPS:.0f} m "
         f"= {db_floor:.0f} per km²", transform=a2.transAxes, ha="right", va="top", color="#CC6677", fontsize=8,
         bbox={"fc": "white", "ec": "none", "alpha": 0.9})
a2.set_yscale("log"); a2.grid(alpha=0.3, axis="y")
a2.set_xlabel("HDBSCAN clusters, densest first"); a2.set_ylabel("places per km² of convex hull (log)")
a2.set_title(f"Cluster densities span a {dens.max() / dens.min():.0f}x range")
fig.savefig(OUT / "04-hdbscan-condensed-tree.png"); plt.close(fig)

# ---------------------------------------------------------------- 6. hulls, profiles and noise
section("cluster polygons: convex vs concave hulls")
rows = []
for c, grp in food[food.hdb >= 0].groupby("hdb"):
    mp = shapely.MultiPoint(grp.geometry.tolist())
    convex = shapely.convex_hull(mp)
    concave = shapely.concave_hull(mp, ratio=0.3)
    top = grp.basic_category.value_counts()
    rows.append({"cluster": c, "n": len(grp), "district": grp.district.mode().iloc[0] if grp.district.notna().any()
                 else "", "convex_km2": convex.area / 1e6, "concave_km2": concave.area / 1e6,
                 "top": f"{top.index[0].replace('_', ' ')} {top.iloc[0] / len(grp):.0%}",
                 "bars": (grp.basic_category == "bar").mean(), "geometry": concave, "convex": convex})
hulls = gpd.GeoDataFrame(rows, geometry="geometry", crs=CRS)
hulls["density"] = hulls.n / hulls.concave_km2.clip(lower=0.005)
print(f"{len(hulls)} clusters: convex hulls cover {hulls.convex_km2.sum():.2f} km2, "
      f"concave hulls (ratio 0.3) {hulls.concave_km2.sum():.2f} km2 "
      f"({hulls.concave_km2.sum() / hulls.convex_km2.sum():.0%} of convex)")
print(f"clustered places: {hulls.n.sum():,} of {len(food):,} ({hulls.n.sum() / len(food):.0%}) "
      f"in {hulls.concave_km2.sum() / (city.area / 1e6):.1%} of the city's area")
print("ten largest clusters:")
print(hulls.sort_values("n", ascending=False).head(10)[["cluster", "n", "district", "concave_km2", "convex_km2",
                                                         "top", "bars"]]
      .to_string(index=False, float_format=lambda v: f"{v:.2f}"))
noise_pts = food[food.hdb < 0]
print(f"noise: {len(noise_pts):,} places; most common categories: "
      + ", ".join(f"{k.replace('_', ' ')} {v:.0%}" for k, v in noise_pts.basic_category.value_counts(normalize=True)
                 .head(3).items()))
nn_noise = BallTree(X_m[hl >= 0]).query(X_m[hl < 0], k=1)[0][:, 0]
print(f"distance from each noise place to the nearest clustered place: median {np.median(nn_noise):.0f} m, "
      f"{np.mean(nn_noise > 500):.0%} more than 500 m away")
hulls.drop(columns="convex").to_file(OUT / "glasgow_dining_clusters.gpkg", layer="hdbscan_concave_hulls",
                                     driver="GPKG")
food[["id", "name", "basic_category", "district", "hdb", "geometry"]].assign(dbscan=lab).to_file(
    OUT / "glasgow_dining_clusters.gpkg", layer="places", driver="GPKG")

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), layout="constrained")
for ax, which, title in [(axes[0], "convex", "Convex hulls"), (axes[1], "geometry", "Concave hulls (ratio = 0.3)")]:
    ax.scatter(noise_pts.geometry.x, noise_pts.geometry.y, s=3, color=NOISE, lw=0)
    gs = gpd.GeoSeries(hulls[which], crs=CRS)
    gs.plot(ax=ax, color=[TOL[c % len(TOL)] for c in hulls.cluster], alpha=0.35, edgecolor="black", lw=0.6)
    cl = food[food.hdb >= 0]
    ax.scatter(cl.geometry.x, cl.geometry.y, s=3, c=cluster_colours(cl.hdb), lw=0)
    for _, r in hulls[hulls.n >= 60].iterrows():
        p = r.geometry.representative_point()
        if zoom[0] < p.x < zoom[1] and zoom[2] < p.y < zoom[3]:
            ax.annotate(f"{r.district} ({r.n})", (p.x, p.y), fontsize=7, ha="center",
                        bbox={"boxstyle": "round,pad=0.15", "fc": "white", "ec": "none", "alpha": 0.8})
    ax.set_xlim(zoom[:2]); ax.set_ylim(zoom[2:]); ax.set_aspect("equal"); ax.set_xticks([]); ax.set_yticks([])
    ax.set_title(f"{title}: {gs.area.sum() / 1e6:.2f} km² in total")
    ax.plot([zoom[0] + 300, zoom[0] + 1300], [zoom[2] + 300] * 2, color="black", lw=2)
    ax.text(zoom[0] + 800, zoom[2] + 420, "1 km", ha="center", fontsize=7)
axes[1].legend(handles=[Patch(fc=TOL[1], alpha=0.5, ec="black", label="HDBSCAN cluster polygon"),
                        Line2D([], [], marker="o", ls="", color=NOISE, label="noise place")],
               loc="upper right", fontsize=7)
fig.suptitle("HDBSCAN dining districts as polygons; labels give the most common postcode district and the "
             "number of places (clusters of 60+)")
fig.savefig(OUT / "05-cluster-hulls.png"); plt.close(fig)

# ---------------------------------------------------------------- 7. internal validation (with caveats)
section("internal validation: silhouette and DBCV")
val = {}
for name, l in [("DBSCAN", lab), ("HDBSCAN", hl)]:
    m = l >= 0
    sil = silhouette_score(X_m[m], l[m], random_state=SEED)
    dbcv = hdbscan.validity.validity_index(X_m.astype(np.float64), l.astype(np.intp))
    val[name] = {"silhouette": round(float(sil), 3), "dbcv": round(float(dbcv), 3), "clusters": int(l.max() + 1),
                 "noise": round(float(np.mean(~m)), 3)}
    print(f"{name:8s}: silhouette (non-noise only) {sil:.3f} | DBCV {dbcv:.3f} | {l.max() + 1} clusters | "
          f"noise {np.mean(~m):.1%}")
rng = np.random.default_rng(SEED)
lo, hi = X_m.min(0), X_m.max(0)
U = rng.uniform(lo, hi, size=(len(X_m), 2))
U = U[gpd.GeoSeries(gpd.points_from_xy(U[:, 0], U[:, 1]), crs=CRS).within(city).to_numpy()]
ul = HDBSCAN(min_cluster_size=MIN_CLUSTER, min_samples=MIN_SAMPLES).fit(U).labels_
print(f"null check: {len(U):,} uniformly random points in the city -> HDBSCAN finds {ul.max() + 1} clusters, "
      f"noise {np.mean(ul < 0):.0%}" + (f", DBCV {hdbscan.validity.validity_index(U, ul.astype(np.intp)):.3f}"
                                         if ul.max() >= 1 else ""))
udb = DBSCAN(eps=EPS, min_samples=MIN_SAMPLES).fit(U).labels_
print(f"null check: DBSCAN (eps={EPS:.0f} m, min_samples={MIN_SAMPLES}) on the same random points -> "
      f"{udb.max() + 1} clusters, noise {np.mean(udb < 0):.0%}")

# ---------------------------------------------------------------- 8. scaling with a BallTree
section("scaling: BallTree neighbourhoods")
allpts = np.radians(g[["lat", "lon"]].to_numpy())
print(f"all Overture places in Glasgow City: {len(allpts):,}")
t0 = time.perf_counter(); tree = BallTree(allpts, metric="haversine"); t_build = time.perf_counter() - t0
t0 = time.perf_counter(); counts = tree.query_radius(allpts, r=EPS / R_EARTH, count_only=True)
t_q = time.perf_counter() - t0
print(f"BallTree build {t_build * 1000:.0f} ms; {len(allpts):,} radius queries of {EPS:.0f} m in {t_q * 1000:.0f} ms; "
      f"neighbours per point median {int(np.median(counts))}, max {counts.max():,}")
print(f"a dense haversine distance matrix for all of them would need {len(allpts) ** 2 * 8 / 1e9:.1f} GB")
sizes_n = [2_000, 4_000, 8_000, 16_000, 32_000, len(allpts)]
DBSCAN(eps=EPS / R_EARTH, min_samples=MIN_SAMPLES, metric="haversine", algorithm="brute").fit(allpts[:500])  # warm-up
timing = []
for n in sizes_n:
    idx = rng.choice(len(allpts), n, replace=False)
    Xs = allpts[idx]
    row = {"n": n}
    for algo in ("ball_tree", "brute"):
        if algo == "brute" and n > 16_000:
            continue
        t0 = time.perf_counter()
        DBSCAN(eps=EPS / R_EARTH, min_samples=MIN_SAMPLES, metric="haversine", algorithm=algo).fit(Xs)
        row[algo] = time.perf_counter() - t0
    t0 = time.perf_counter()
    HDBSCAN(min_cluster_size=MIN_CLUSTER, min_samples=MIN_SAMPLES, metric="haversine").fit(Xs)
    row["hdbscan"] = time.perf_counter() - t0
    timing.append(row)
    print(f"n={n:>6,}: DBSCAN ball_tree {row['ball_tree']:.2f} s"
          + (f" | brute {row['brute']:.2f} s" if "brute" in row else "") + f" | HDBSCAN {row['hdbscan']:.2f} s")
tdf = pd.DataFrame(timing)
t0 = time.perf_counter()
HDBSCAN(min_cluster_size=MIN_CLUSTER, min_samples=MIN_SAMPLES).fit(np.column_stack([g.geometry.x, g.geometry.y]))
print(f"HDBSCAN on all {len(g):,} places in projected metres (KD-tree): {time.perf_counter() - t0:.2f} s")
fig, ax = plt.subplots(figsize=(6.5, 4))
for col, c, lbl in [("ball_tree", "#332288", "DBSCAN, BallTree (haversine)"),
                    ("brute", "#CC6677", "DBSCAN, brute force (haversine)"),
                    ("hdbscan", "#117733", "HDBSCAN (haversine)")]:
    d = tdf.dropna(subset=[col])
    ax.plot(d.n, d[col], marker="o", color=c, label=lbl)
ax.set_xscale("log"); ax.set_yscale("log"); ax.grid(alpha=0.3, which="both")
ax.set_xticks(sizes_n, [f"{n / 1000:.0f}k" for n in sizes_n]); ax.minorticks_off()
ax.set_xlabel("number of points (random samples of all Glasgow places)"); ax.set_ylabel("run time (s, log)")
ax.set_title(f"Run time with eps = {EPS:.0f} m, min_samples = {MIN_SAMPLES} (single thread)")
ax.legend(fontsize=8)
fig.savefig(OUT / "06-scaling-balltree.png"); plt.close(fig)

json.dump({"n_places": len(food), "eps_m": EPS, "min_samples": MIN_SAMPLES, "validation": val,
           "timing": timing}, open(OUT / "metrics.json", "w"), indent=2)
print("\nwrote output/glasgow_dining_clusters.gpkg, output/metrics.json and 6 figures")