import warnings
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.ticker import FuncFormatter, PercentFormatter
from scipy.stats import spearmanr, mannwhitneyu

try:
    import geopandas as gpd
    import shapely
except ImportError:
    raise ImportError("Install geopandas first:  pip install geopandas  (it is pre-installed on Kaggle)")

warnings.filterwarnings("ignore")
pd.set_option("display.width", 160)
pd.set_option("display.max_columns", 12)
pd.set_option("display.float_format", "{:,.2f}".format)

# 🎨 palette and chart style
BLUE, ORANGE, AQUA, VIOLET = "#2a78d6", "#eb6834", "#1baf7a", "#6250d6"
INK, MUTED, GRID, SURFACE, NODATA = "#0b0b0b", "#52514e", "#e5e4df", "#fcfcfb", "#d9d8d3"
plt.rcParams.update({
    "figure.dpi": 110, "figure.facecolor": SURFACE, "axes.facecolor": SURFACE,
    "axes.edgecolor": GRID, "axes.labelcolor": MUTED, "axes.titlecolor": INK,
    "axes.titlesize": 11, "axes.titleweight": "bold", "axes.titlelocation": "left",
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.color": GRID, "grid.linewidth": .8,
    "xtick.color": MUTED, "ytick.color": MUTED, "font.size": 10,
    "legend.frameon": False, "lines.linewidth": 2, "hatch.color": "#9a9993",
})

# 🏷️ energy labels, best to worst, in the colours used on Dutch label certificates
LABELS = ["A++", "A+", "A", "B", "C", "D", "E", "F", "G"]
LABEL_COLORS = ["#00843d", "#2aa13c", "#55b947", "#bed630", "#ffdd00", "#fdb913", "#f37021", "#ed1c24", "#b5121b"]
LABEL_SCORE = dict(zip(LABELS, range(9, 0, -1)))          # A++ = 9 ... G = 1

# 📁 paths: works locally (./data) and on Kaggle (/kaggle/input/...)
def find_file(*patterns):
    for pattern in patterns:
        for base in [Path("data"), Path("."), Path("/kaggle/input")]:
            hits = sorted(base.rglob(pattern)) if base.exists() else []
            if hits:
                return hits[0]
    raise FileNotFoundError(patterns)

OUT_DIR = Path("/kaggle/working") if Path("/kaggle/working").exists() else Path("output")
FIG_DIR = OUT_DIR / "figures"
FIG_DIR.mkdir(parents=True, exist_ok=True)

# 🖼️ numbered figures, saved as PNG
_FIG = {"n": 0}
def figure_title(fig, text, y=1.0):
    _FIG["n"] += 1
    fig.suptitle(f"Figure {_FIG['n']} · {text}", x=0.01, y=y, ha="left", fontsize=13, fontweight="bold", color=INK)
    return _FIG["n"]

def save(fig, name):
    fig.savefig(FIG_DIR / f"fig{_FIG['n']:02d}_{name}.png", dpi=150, bbox_inches="tight", facecolor=SURFACE)
    plt.show()

eur = FuncFormatter(lambda v, _: f"€{v/1000:,.0f}k")
pct = PercentFormatter(100, decimals=0)

# 🗺️ map helpers
def map_axis(ax):
    ax.set_axis_off()
    ax.set_aspect("equal")

def classed_colors(values, bins, cmap):
    """Colour each value by the class it falls in (bins = class edges); NaN -> no-data grey."""
    from matplotlib.colors import ListedColormap
    k = len(bins) - 1
    cm = ListedColormap(plt.get_cmap(cmap)(np.linspace(0.2, 0.95, k)))   # skip the near-white end
    cls = np.clip(np.digitize(values, bins[1:-1]), 0, len(bins) - 2)
    return [NODATA if pd.isna(v) else cm(c) for v, c in zip(values, cls)], cm

def choropleth(gdf, col, ax, cmap, k=5, fmt="{:.0f}", title=None, bins=None):
    """Quantile choropleth with a compact class legend; missing values hatched."""
    vals = gdf[col]
    bins = np.nanquantile(vals, np.linspace(0, 1, k + 1)) if bins is None else np.asarray(bins)
    colors, cm = classed_colors(vals, bins, cmap)
    gdf.plot(ax=ax, color=colors, edgecolor="white", linewidth=0.4)
    miss = gdf[vals.isna()]
    if len(miss):
        miss.plot(ax=ax, color=NODATA, edgecolor="white", linewidth=0.4, hatch="////")
    handles = [Patch(color=cm(i), label=f"{fmt.format(bins[i])} – {fmt.format(bins[i+1])}") for i in range(len(bins) - 1)]
    if len(miss):
        handles.append(Patch(facecolor=NODATA, hatch="////", edgecolor="#9a9993", label="no data"))
    ax.legend(handles=handles, loc="upper left", fontsize=8, title=title, title_fontsize=8, handlelength=1.2)
    map_axis(ax)

# 📐 statistics helpers
def fmt_p(p):
    return "< 0.001" if p < 0.001 else f"{p:.3f}"

def spearman_ci(x, y, n_boot=2000, seed=42):
    """Spearman rho with a percentile bootstrap 95% CI."""
    x, y = np.asarray(x), np.asarray(y)
    rho, p = spearmanr(x, y)
    rng = np.random.default_rng(seed)
    idx = rng.integers(0, len(x), (n_boot, len(x)))
    boots = [spearmanr(x[i], y[i])[0] for i in idx]
    lo, hi = np.nanpercentile(boots, [2.5, 97.5])
    return rho, p, lo, hi

def queen_weights(gdf, decimals=7):
    """Row-standardised contiguity matrix: two areas are neighbours if they share a boundary point."""
    xy, idx = shapely.get_coordinates(gdf.geometry.to_numpy(), return_index=True)
    key = pd.Series(idx).groupby([np.round(xy[:, 0], decimals), np.round(xy[:, 1], decimals)]).unique()
    n = len(gdf)
    W = np.zeros((n, n))
    for ids in key:
        if len(ids) > 1:
            for i in ids:
                W[i, ids] = 1
    np.fill_diagonal(W, 0)
    return W / W.sum(1, keepdims=True).clip(min=1)

def morans_i(y, W, n_perm=999, seed=42):
    """Global Moran's I with a one-sided permutation p-value (clustering)."""
    z = np.asarray(y, float) - np.mean(y)
    I = len(z) / W.sum() * (z @ W @ z) / (z @ z)
    rng = np.random.default_rng(seed)
    sims = np.array([(lambda s: len(s) / W.sum() * (s @ W @ s) / (s @ s))(rng.permutation(z)) for _ in range(n_perm)])
    p = (1 + (sims >= I).sum()) / (n_perm + 1)
    return I, p, sims

print("✅ setup done · figures go to", FIG_DIR)

CSV_FILE = find_file("neighborhoods.csv", "*neighborhoods*.csv")
raw = pd.read_csv(CSV_FILE, index_col=0)
raw.index.name = "neighbourhood"
print(f"{CSV_FILE.name}: {raw.shape[0]} rows × {raw.shape[1]} columns")
print("Columns:", [c.replace("Energielabel ", "") for c in raw.columns[:-1]], "+ income")
raw.tail(3)

# 🗺️ boundaries: GeoJSON preferred (one file); shapefile as a fallback
GEO_FILE = find_file("neighborhoods.json", "*neighborhoods*.json", "neighborhoods.shp", "*neighborhoods*.shp")
geo_raw = gpd.read_file(GEO_FILE)
# the code+name column is 'neigb_cijfers' in the GeoJSON and 'neigb_cijf' in the shapefile (10-character limit)
key_col = [c for c in geo_raw.columns if c.startswith("neigb_cijf")][0]
geo_raw = geo_raw.rename(columns={key_col: "neighbourhood", "neighb_cbs": "name_cbs"})
print(f"{GEO_FILE.name}: {len(geo_raw)} shapes · CRS {geo_raw.crs}")

# 🔗 does every table row find a shape, and vice versa?
in_table = set(raw.index) - {"Unknown"}
in_map = set(geo_raw["neighbourhood"])
print(f"Table rows without a shape : {sorted(in_table - in_map) or 'none'}")
print(f"Shapes without a table row : {sorted(in_map - in_table) or 'none'}")
geo_raw.head(3)[["neighbourhood", "name_cbs"]]

df = raw.drop(index="Unknown").copy()
df.columns = [c.replace("Energielabel ", "") for c in df.columns[:-1]] + ["income"]
df = df[LABELS + ["income"]]                      # reorder by name: A++ … G

df["labelled_homes"] = df[LABELS].sum(axis=1)
shares = df[LABELS].div(df["labelled_homes"], axis=0) * 100
df["pct_A_or_better"] = shares[["A++", "A+", "A"]].sum(axis=1)
df["pct_B_to_D"] = shares[["B", "C", "D"]].sum(axis=1)
df["pct_E_to_G"] = shares[["E", "F", "G"]].sum(axis=1)
df["label_score"] = (df[LABELS] * pd.Series(LABEL_SCORE)).sum(axis=1) / df["labelled_homes"]

df["code"] = df.index.str.split(" ").str[0].astype(int)
df["area"] = np.where(df["code"] >= 101, "New-build (Leidschenveen-Ypenburg)", "Rest of the city")

print(f"{len(df)} neighbourhoods · {df['labelled_homes'].sum():,.0f} homes with a registered label")
print(f"Income missing for {df['income'].isna().sum()} neighbourhoods")
df.sort_values("labelled_homes")[["labelled_homes", "income", "pct_A_or_better"]].head(8)

MIN_HOMES = 100
df["in_analysis"] = (df["labelled_homes"] >= MIN_HOMES) & df["income"].notna()
HIGH_INCOME = df.loc[df["in_analysis"], "income"].quantile(0.75)
df["high_income"] = df["income"] > HIGH_INCOME

# 🔗 join to the boundaries and project to the Dutch national grid (RD New, metres) for mapping
nb = geo_raw.merge(df, left_on="neighbourhood", right_index=True, how="left").to_crs(28992)
an = df[df["in_analysis"]].copy()

print(f"In the analysis: {len(an)} of {len(df)} neighbourhoods "
      f"({an['labelled_homes'].sum() / df['labelled_homes'].sum():.1%} of labelled homes)")
print(f"High-income threshold (75th percentile): €{HIGH_INCOME:,.0f} → {an['high_income'].sum()} neighbourhoods")
print(f"New-build neighbourhoods in the analysis: {(an['area'] != 'Rest of the city').sum()}")

fig, ax = plt.subplots(figsize=(8, 7.5))
choropleth(nb.assign(labelled_homes=nb["labelled_homes"].where(nb["in_analysis"])), "labelled_homes", ax,
           "Greys", k=4, fmt="{:,.0f}", title="Homes with a label")
pts = nb[~nb["in_analysis"]].geometry.representative_point()
for (x, y), name in zip(zip(pts.x, pts.y), nb.loc[~nb["in_analysis"], "neighbourhood"]):
    ax.annotate(name.split(" ", 1)[1], (x, y), fontsize=7, color=ORANGE, ha="center", fontweight="bold")
figure_title(fig, f"{(~nb['in_analysis']).sum()} areas are left out: almost no homes or no income figure", y=0.97)
save(fig, "coverage_map")

order = an.sort_values("label_score", ascending=False)
sh = order[LABELS].div(order["labelled_homes"], axis=0) * 100

fig, (ax, ax2) = plt.subplots(2, 1, figsize=(13, 5.6), height_ratios=[4, 1], sharex=True)
bottom = np.zeros(len(sh))
x = np.arange(len(sh))
for lab, col in zip(LABELS, LABEL_COLORS):
    ax.bar(x, sh[lab], bottom=bottom, color=col, width=1.0, label=lab, linewidth=0)
    bottom += sh[lab].to_numpy()
ax.set_ylim(0, 100); ax.yaxis.set_major_formatter(pct); ax.grid(False)
ax.set_ylabel("Share of labelled homes")
ax.legend(ncol=9, loc="upper left", bbox_to_anchor=(0, 1.13), fontsize=9, handlelength=1, columnspacing=1)
new = (order["area"] != "Rest of the city").to_numpy()
ax2.bar(x, order["income"], color=np.where(new, AQUA, "#b9b8b2"), width=0.85)
ax2.yaxis.set_major_formatter(eur); ax2.set_ylabel("Income"); ax2.set_xticks([])
ax2.set_xlabel("Neighbourhoods, best average label → worst     (aqua = new-build district)")
figure_title(fig, "A few neighbourhoods are almost all A-label; most of the city is D or worse", y=1.04)
save(fig, "label_mix_strip")

fig, axes = plt.subplots(1, 3, figsize=(17, 6.4))
view = nb.copy()
for c in ["income", "pct_A_or_better", "pct_E_to_G"]:
    view[c] = view[c].where(view["in_analysis"])
choropleth(view, "income", axes[0], "Blues", fmt="€{:,.0f}", title="Income (quintiles)")
choropleth(view, "pct_A_or_better", axes[1], "Greens", fmt="{:.0f}%", title="% A-label or better")
choropleth(view, "pct_E_to_G", axes[2], "Reds", fmt="{:.0f}%", title="% E, F or G")
for ax, t in zip(axes, ["Household income", "Homes with label A or better", "Homes with label E–G"]):
    ax.set_title(t, fontsize=11)
figure_title(fig, "The richest areas, along the dunes in the north-west, have the most E–G homes", y=1.0)
fig.tight_layout()
save(fig, "maps_side_by_side")

BIV = np.array([["#e8e8e8", "#dfb0d6", "#be64ac"],     # rows: income tercile low→high
                ["#ace4e4", "#a5add3", "#8c62aa"],     # cols: % A-or-better tercile low→high
                ["#5ac8c8", "#5698b9", "#3b4994"]])
view = nb[nb["in_analysis"]].copy()
inc_t = pd.qcut(view["income"], 3, labels=False)
lab_t = pd.qcut(view["pct_A_or_better"], 3, labels=False)
view["biv"] = [BIV[i, j] for i, j in zip(inc_t, lab_t)]

fig, ax = plt.subplots(figsize=(9, 8))
nb[~nb["in_analysis"]].plot(ax=ax, color=NODATA, edgecolor="white", linewidth=0.4, hatch="////")
view.plot(ax=ax, color=view["biv"].tolist(), edgecolor="white", linewidth=0.4)
map_axis(ax)

# 3×3 legend
lax = fig.add_axes([0.16, 0.66, 0.15, 0.15])
for i in range(3):
    for j in range(3):
        lax.add_patch(plt.Rectangle((j, i), 1, 1, color=BIV[i, j]))
lax.set_xlim(0, 3); lax.set_ylim(0, 3); lax.set_aspect("equal"); lax.grid(False)
lax.set_xticks([]); lax.set_yticks([])
lax.set_xlabel("% A or better →", fontsize=8); lax.set_ylabel("Income →", fontsize=8)
for s in lax.spines.values():
    s.set_visible(False)

counts = pd.crosstab(inc_t, lab_t)
figure_title(fig, "No diagonal: rich areas are as likely to have few A-labels as many", y=0.95)
save(fig, "bivariate_map")
print("Neighbourhoods per cell (rows: income tercile low→high, cols: % A-or-better tercile low→high)")
print(counts.to_string())

fig, axes = plt.subplots(1, 2, figsize=(14, 5.6), sharex=True)
for ax, (col, ylab) in zip(axes, [("pct_A_or_better", "Homes with label A or better"), ("pct_E_to_G", "Homes with label E, F or G")]):
    for area, color in [("Rest of the city", "#7d7c77"), ("New-build (Leidschenveen-Ypenburg)", AQUA)]:
        s = an[an["area"] == area]
        rho, p = spearmanr(s["income"], s[col])
        ax.scatter(s["income"], s[col], s=np.sqrt(s["labelled_homes"]) * 2.2, color=color, alpha=.75,
                   edgecolor="white", linewidth=.6, label=f"{area} (n = {len(s)}, ρ = {rho:+.2f})")
    rho, p = spearmanr(an["income"], an[col])
    ax.axvline(HIGH_INCOME, color=ORANGE, ls="--", lw=1.2)
    ax.text(HIGH_INCOME, 101, " high income →", color=ORANGE, fontsize=8, va="bottom")
    ax.set_ylim(-3, 108); ax.yaxis.set_major_formatter(pct); ax.xaxis.set_major_formatter(eur)
    ax.set_xlabel("Average standardised disposable household income")
    ax.set_title(f"{ylab}   ·   all neighbourhoods ρ = {rho:+.2f} (p = {fmt_p(p)})", fontsize=10)
    ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.13), fontsize=8, ncol=1)
for name in an.nlargest(5, "pct_A_or_better").index.tolist() + ["74 Duttendel", "68 Marlot"]:
    r = an.loc[name]
    axes[0].annotate(name.split(" ", 1)[1], (r["income"], r["pct_A_or_better"]), xytext=(5, 2),
                     textcoords="offset points", fontsize=7, color=MUTED)
figure_title(fig, "Income does not predict good labels; outside the new-build district it predicts bad ones", y=1.02)
fig.tight_layout()
save(fig, "scatter_income_labels")

def quartile_mix(frame):
    q = pd.qcut(frame["income"], 4, labels=["Q1\nlowest", "Q2", "Q3", "Q4\nhighest"])
    tot = frame.groupby(q, observed=True)[LABELS].sum()          # pool homes within each quartile
    rng = frame.groupby(q, observed=True)["income"].agg(["min", "max"])
    return tot.div(tot.sum(axis=1), axis=0) * 100, rng

fig, axes = plt.subplots(1, 2, figsize=(14, 4.6), sharey=True)
for ax, (title, frame) in zip(axes, [("All neighbourhoods", an), ("Without the new-build district", an[an["area"] == "Rest of the city"])]):
    mix, rng = quartile_mix(frame)
    left = np.zeros(len(mix))
    for lab, col in zip(LABELS, LABEL_COLORS):
        ax.barh(range(len(mix)), mix[lab], left=left, color=col, height=.7, label=lab)
        left += mix[lab].to_numpy()
    for i, (a, e) in enumerate(zip(mix[["A++", "A+", "A"]].sum(axis=1), mix[["E", "F", "G"]].sum(axis=1))):
        ax.text(1, i, f"≥A: {a:.0f}%", va="center", fontsize=8, color="white", fontweight="bold")
        ax.text(99, i, f"E–G: {e:.0f}%", va="center", ha="right", fontsize=8, color="white", fontweight="bold")
    ax.set_yticks(range(len(mix)), [f"{q}  €{lo/1000:.0f}–{hi/1000:.0f}k".replace("\n", " ") for q, (lo, hi) in zip(mix.index, rng.to_numpy())])
    ax.invert_yaxis(); ax.set_xlim(0, 100); ax.xaxis.set_major_formatter(pct); ax.grid(False)
    ax.set_title(f"{title} (n = {len(frame)})")
axes[1].legend(ncol=9, loc="upper right", bbox_to_anchor=(1, -0.1), fontsize=8, handlelength=1)
figure_title(fig, "The richest quarter has the worst label mix, with or without the new-build district", y=1.03)
fig.tight_layout()
save(fig, "quartile_label_mix")

rows = []
for col, label in [("pct_A_or_better", "% A-label or better"), ("pct_E_to_G", "% E–G"), ("label_score", "Average label score (1–9)")]:
    hi, lo = an.loc[an["high_income"], col], an.loc[~an["high_income"], col]
    u, p = mannwhitneyu(hi, lo, alternative="two-sided")
    rows.append({"measure": label, "high income (median)": hi.median(), "rest (median)": lo.median(),
                 "direction": "high income better" if (hi.median() > lo.median()) == (col != "pct_E_to_G") else "high income worse",
                 "Mann-Whitney p": fmt_p(p)})
h1_table = pd.DataFrame(rows).set_index("measure")
print(f"High income: n = {an['high_income'].sum()} · rest: n = {(~an['high_income']).sum()}")
h1_table

groups = {"All neighbourhoods": an,
          "Rest of the city": an[an["area"] == "Rest of the city"],
          "New-build district": an[an["area"] != "Rest of the city"]}
measures = {"pct_A_or_better": "% A or better", "label_score": "Average label score", "pct_E_to_G": "% E–G"}
res = []
for g, frame in groups.items():
    for col, mlabel in measures.items():
        rho, p, lo, hi = spearman_ci(frame["income"], frame[col])
        res.append({"group": g, "measure": mlabel, "n": len(frame), "rho": rho, "ci_low": lo, "ci_high": hi, "p": fmt_p(p)})
corr = pd.DataFrame(res)

fig, ax = plt.subplots(figsize=(10, 5.2))
colors = {"All neighbourhoods": INK, "Rest of the city": "#7d7c77", "New-build district": AQUA}
ypos, yt = 0, []
for m in measures.values():
    for g in groups:
        r = corr[(corr["group"] == g) & (corr["measure"] == m)].iloc[0]
        ax.plot([r["ci_low"], r["ci_high"]], [ypos, ypos], color=colors[g], lw=2.5, solid_capstyle="round")
        ax.scatter(r["rho"], ypos, color=colors[g], s=60, zorder=3, edgecolor="white")
        ax.text(1.04, ypos, f"{r['rho']:+.2f}  (n = {r['n']})", va="center", fontsize=8, color=colors[g])
        yt.append((ypos, f"{m} · {g}"))
        ypos -= 1
    ypos -= 0.6
ax.axvline(0, color=MUTED, lw=1)
ax.set_yticks([y for y, _ in yt], [t for _, t in yt], fontsize=8.5)
ax.set_xlim(-1, 1.0); ax.set_xlabel("Spearman ρ with household income (95% bootstrap CI)")
ax.grid(axis="y", visible=False)
figure_title(fig, "In the existing city, richer means worse labels; the citywide 'no effect' is two effects cancelling", y=1.0)
save(fig, "forest_correlations")
corr.round(2).set_index(["measure", "group"])

geo_an = geo_raw.merge(an, left_on="neighbourhood", right_index=True, how="inner").reset_index(drop=True)
W = queen_weights(geo_an)
n_nb = (W > 0).sum(axis=1)
print(f"Neighbours per area: mean {n_nb.mean():.1f}, min {n_nb.min()}, max {n_nb.max()} · areas with none: {(n_nb == 0).sum()}")

moran = {}
fig, axes = plt.subplots(1, 3, figsize=(15, 3.8), sharey=True)
for ax, (col, label, color) in zip(axes, [("income", "Household income", BLUE),
                                          ("pct_A_or_better", "% A-label or better", AQUA),
                                          ("pct_E_to_G", "% E–G", ORANGE)]):
    I, p, sims = morans_i(geo_an[col], W)
    moran[label] = {"Moran's I": I, "p (999 permutations)": fmt_p(p), "expected if random": -1 / (len(W) - 1)}
    ax.hist(sims, bins=30, color="#c9c8c2")
    ax.axvline(I, color=color, lw=3)
    ax.text(I, ax.get_ylim()[1] * .9, f"  I = {I:.2f}", color=color, fontweight="bold", fontsize=10)
    ax.set_title(label); ax.set_xlabel("Moran's I")
    ax.set_xlim(-0.35, max(0.8, I + 0.1))
axes[0].set_ylabel("Shuffled maps")
figure_title(fig, "All three measures are strongly clustered on the map", y=1.06)
fig.tight_layout()
save(fig, "morans_i")
pd.DataFrame(moran).T

top_A = an.nlargest(5, "pct_A_or_better")
retrofit = an[an["high_income"]].nlargest(5, "pct_E_to_G")
cols = ["income", "labelled_homes", "pct_A_or_better", "pct_E_to_G", "label_score", "area"]
print("🟢 Highest share of A-label homes")
display(top_A[cols])
print("🔴 High income, most E–G homes (retrofit candidates)")
display(retrofit[cols])

fig, (ax, tx) = plt.subplots(1, 2, figsize=(13, 8), width_ratios=[2.2, 1])
nb.plot(ax=ax, color="#ecebe6", edgecolor="white", linewidth=0.5)
groups5 = [("A", top_A, LABEL_COLORS[1], LABEL_COLORS[0], "🟢 Highest share of A-label homes", "pct_A_or_better", "≥A"),
           ("R", retrofit, LABEL_COLORS[7], LABEL_COLORS[8], "🔴 Retrofit candidates: high income, most E–G", "pct_E_to_G", "E–G")]
tx.set_axis_off(); ty = 0.95
for tag, frame, fill, dark, head, col, unit in groups5:
    sel = nb[nb["neighbourhood"].isin(frame.index)]
    sel.plot(ax=ax, color=fill, edgecolor="white", linewidth=0.8)
    pts = sel.geometry.representative_point()
    pos = dict(zip(sel["neighbourhood"], zip(pts.x, pts.y)))
    tx.text(0, ty, head.split(" ", 1)[1], fontsize=11, fontweight="bold", color=dark, transform=tx.transAxes); ty -= 0.06
    for k, (name, r) in enumerate(frame.iterrows(), 1):
        x, y = pos[name]
        ax.scatter(x, y, s=260, color="white", edgecolor=dark, linewidth=1.8, zorder=4)
        ax.text(x, y, f"{tag}{k}", ha="center", va="center", fontsize=7, fontweight="bold", color=dark, zorder=5)
        tx.text(0, ty, f"{tag}{k}", fontsize=9, fontweight="bold", color=dark, va="top", transform=tx.transAxes)
        tx.text(0.12, ty, f"{name.split(' ', 1)[1]}\n€{r['income']:,.0f} · {r[col]:.0f}% {unit} · {r['labelled_homes']:,.0f} homes",
                fontsize=9, color=INK, va="top", transform=tx.transAxes, linespacing=1.4)
        ty -= 0.075
    ty -= 0.05
map_axis(ax)
figure_title(fig, "The best labels are in new-build areas; the richest dune-side areas have the worst", y=0.97)
fig.tight_layout()
save(fig, "five_neighbourhoods")

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, silhouette_samples, adjusted_rand_score
from scipy.stats import kruskal

FEATURES = ["≥A", "B", "C", "D", "E", "F", "G"]
FEATURE_COLORS = [LABEL_COLORS[2]] + LABEL_COLORS[3:]

cl = df[df["labelled_homes"] >= MIN_HOMES].copy()
X = cl[LABELS].div(cl["labelled_homes"], axis=0) * 100
X["≥A"] = X[["A++", "A+", "A"]].sum(axis=1)
X = X[FEATURES]
print(f"{len(X)} neighbourhoods × {X.shape[1]} features")

K_RANGE = range(2, 11)
scan = []
for k in K_RANGE:
    base = KMeans(n_clusters=k, n_init=50, random_state=42).fit(X)
    ari = [adjusted_rand_score(base.labels_, KMeans(n_clusters=k, n_init=50, random_state=s).fit(X).labels_)
           for s in range(1, 11)]
    scan.append({"k": k, "inertia": base.inertia_, "silhouette": silhouette_score(X, base.labels_),
                 "stability (min ARI)": min(ari), "smallest cluster": np.bincount(base.labels_).min()})
scan = pd.DataFrame(scan).set_index("k")

fig, axes = plt.subplots(1, 3, figsize=(15, 3.8))
for ax, col, color, title in [(axes[0], "inertia", MUTED, "Inertia (elbow)"),
                              (axes[1], "silhouette", BLUE, "Silhouette score ↑"),
                              (axes[2], "stability (min ARI)", AQUA, "Stability across 10 seeds ↑")]:
    ax.plot(scan.index, scan[col], marker="o", color=color)
    ax.axvline(5, color=ORANGE, ls="--", lw=1.2)
    ax.set_title(title); ax.set_xlabel("Number of clusters k"); ax.set_xticks(list(K_RANGE))
axes[0].text(5.15, scan["inertia"].max() * .95, "k = 5", color=ORANGE, fontsize=9)
axes[2].set_ylim(0.5, 1.05)
figure_title(fig, "k = 5 is stable and close to the best silhouette; k = 6 and k ≥ 8 depend on the random start", y=1.06)
fig.tight_layout()
save(fig, "choose_k")
scan.round(3)

K = 5
km = KMeans(n_clusters=K, n_init=50, random_state=42).fit(X)
cl["cluster_raw"] = km.labels_

# order clusters from best to worst average label and give them readable names
score_by_cluster = cl.groupby("cluster_raw")["label_score"].mean().sort_values(ascending=False)
centres = pd.DataFrame(km.cluster_centers_, columns=FEATURES).loc[score_by_cluster.index]

def name_type(c):
    if c["≥A"] >= 60:  return "Mostly A"
    if c["≥A"] + c["B"] >= 60: return "Mostly A and B"
    if c["G"] >= 35:  return "G-dominated"
    if c["E"] >= 35:  return "E-dominated"
    return "Mixed, centred on C"

TYPE_NAMES = [name_type(centres.loc[i]) for i in score_by_cluster.index]
assert len(set(TYPE_NAMES)) == K, TYPE_NAMES        # every cluster gets its own name
cl["type"] = pd.Categorical(cl["cluster_raw"].map(dict(zip(score_by_cluster.index, TYPE_NAMES))),
                            categories=TYPE_NAMES, ordered=True)
centres.index = TYPE_NAMES
TYPE_COLORS = dict(zip(TYPE_NAMES, [LABEL_COLORS[0], LABEL_COLORS[3], LABEL_COLORS[4], LABEL_COLORS[6], LABEL_COLORS[8]]))

# 📊 profile chart: the cluster centre = the typical label mix of each type
sizes = cl["type"].value_counts().reindex(TYPE_NAMES)
fig, ax = plt.subplots(figsize=(12, 3.8))
left = np.zeros(K)
for f, col in zip(FEATURES, FEATURE_COLORS):
    ax.barh(range(K), centres[f], left=left, color=col, height=.68, label=f)
    for i, (l, v) in enumerate(zip(left, centres[f])):
        if v >= 7:
            ax.text(l + v / 2, i, f"{v:.0f}%", ha="center", va="center", fontsize=8,
                    color="white" if f in ["≥A", "E", "F", "G"] else INK, fontweight="bold")
    left += centres[f].to_numpy()
ax.set_yticks(range(K), [f"{t}  (n = {n})" for t, n in zip(TYPE_NAMES, sizes)])
ax.invert_yaxis(); ax.set_xlim(0, 100); ax.xaxis.set_major_formatter(pct); ax.grid(False)
ax.legend(ncol=7, loc="upper left", bbox_to_anchor=(0, 1.16), fontsize=9, handlelength=1)
figure_title(fig, "Five types: two A/B types (mostly new-build), one mixed type, and two dominated by one poor label", y=1.08)
save(fig, "cluster_profiles")
centres.round(1)

geo_cl = nb.merge(cl[["type"]], left_on="neighbourhood", right_index=True, how="left")
fig, ax = plt.subplots(figsize=(9.5, 8.5))
geo_cl[geo_cl["type"].isna()].plot(ax=ax, color=NODATA, edgecolor="white", linewidth=0.4, hatch="////")
known = geo_cl[geo_cl["type"].notna()]
known.plot(ax=ax, color=[TYPE_COLORS[t] for t in known["type"]], edgecolor="white", linewidth=0.5)
map_axis(ax)
ax.legend(handles=[Patch(color=TYPE_COLORS[t], label=f"{t} ({sizes[t]})") for t in TYPE_NAMES]
                  + [Patch(facecolor=NODATA, hatch="////", edgecolor="#9a9993", label="fewer than 100 labelled homes")],
          loc="lower left", bbox_to_anchor=(0.98, 0.02), fontsize=9, title="Neighbourhood type", title_fontsize=9)
figure_title(fig, "The types form zones: G in the pre-war centre and north-west, A/B on the new-build edges", y=0.95)
save(fig, "cluster_map")

cl["silhouette"] = silhouette_samples(X, km.labels_)
cl["dominant_label"] = pd.Categorical(X.idxmax(axis=1), categories=FEATURES, ordered=True)

fig, (ax, ax2) = plt.subplots(1, 2, figsize=(15, 5.4), width_ratios=[1.35, 1])
y0 = 0
for t in TYPE_NAMES:
    s = cl.loc[cl["type"] == t, "silhouette"].sort_values()
    ax.barh(np.arange(y0, y0 + len(s)), s, color=TYPE_COLORS[t], height=1.0)
    ax.text(-0.02, y0 + len(s) / 2, t, ha="right", va="center", fontsize=9, color=INK)
    y0 += len(s) + 3
avg = cl["silhouette"].mean()
ax.axvline(avg, color=INK, ls="--", lw=1); ax.text(avg, y0, f" mean {avg:.2f}", fontsize=9, va="bottom")
ax.axvline(0, color=MUTED, lw=.8)
ax.set_yticks([]); ax.grid(axis="y", visible=False); ax.set_xlabel("Silhouette (−1 wrong cluster · 0 on the border · 1 clear fit)")
ax.set_title("Fit of each neighbourhood to its type")

ct = pd.crosstab(cl["type"], cl["dominant_label"]).reindex(columns=FEATURES, fill_value=0)
ax2.imshow(ct.to_numpy(), cmap="Greys", aspect="auto", vmin=0, vmax=ct.to_numpy().max() * 1.6)
for i in range(ct.shape[0]):
    for j in range(ct.shape[1]):
        v = ct.iat[i, j]
        if v:
            ax2.text(j, i, v, ha="center", va="center", fontsize=10, fontweight="bold", color=INK)
ax2.set_xticks(range(len(FEATURES)), FEATURES); ax2.set_yticks(range(K), TYPE_NAMES)
ax2.set_xlabel("Most common label in the neighbourhood"); ax2.grid(False)
ax2.set_title("Type vs most common label")
figure_title(fig, "Most neighbourhoods fit their type; the border cases are mixed areas", y=1.02)
fig.tight_layout()
save(fig, "cluster_evaluation")
print("Five weakest fits (lowest silhouette):")
cl.nsmallest(5, "silhouette")[["type", "dominant_label", "silhouette"]].join(X.round(0))