"""Reproducible spatial diagnostic of London housing affordability.

Source
------
Greater London Authority, Housing in London 2019, Table 3. The outcome is the
published median house-price-to-workplace-earnings ratio for 33 local areas.

Design
------
This is a cross-sectional exploratory spatial data analysis, not a causal
estimate. The script reports distributional contrasts, a deterministic
bootstrap for the Inner/Outer mean difference, Moran's I on official borough
polygons with 99,999 random permutations, binary/row-standardised weight
sensitivity, and a leave-one-out check for the extreme Kensington and Chelsea
value.

Dependencies: Python standard library + Pillow.
"""
import csv
import itertools
import json
import math
import random
import statistics
from pathlib import Path

from PIL import Image, ImageDraw, ImageFont


ROOT = Path(__file__).resolve().parent
OUT = ROOT.parent / "public" / "gis"
OUT.mkdir(parents=True, exist_ok=True)

PAPER = "#f5f4ef"
INK = "#171716"
SOFT = "#6e6b65"
HAIR = "#cfcac0"
COBALT = "#2849a8"
CORAL = "#d56754"
PALE_BLUE = "#dbe3f5"


def font(size, bold=False, serif=False):
    candidates = []
    if serif:
        candidates += [
            "C:/Windows/Fonts/georgiab.ttf" if bold else "C:/Windows/Fonts/georgia.ttf",
            "C:/Windows/Fonts/timesbd.ttf" if bold else "C:/Windows/Fonts/times.ttf",
        ]
    candidates += ["C:/Windows/Fonts/arialbd.ttf" if bold else "C:/Windows/Fonts/arial.ttf"]
    for candidate in candidates:
        if Path(candidate).exists():
            return ImageFont.truetype(candidate, size)
    return ImageFont.load_default()


def percentile(values, q):
    ordered = sorted(values)
    pos = (len(ordered) - 1) * q
    lo, hi = math.floor(pos), math.ceil(pos)
    if lo == hi:
        return ordered[lo]
    return ordered[lo] * (hi - pos) + ordered[hi] * (pos - lo)


def mean(values):
    return statistics.mean(values)


def rings(geometry):
    if geometry["type"] == "Polygon":
        yield from geometry["coordinates"]
    elif geometry["type"] == "MultiPolygon":
        for polygon in geometry["coordinates"]:
            yield from polygon
    else:
        raise ValueError(geometry["type"])


def vertices_and_segments(geometry):
    vertices, segments = set(), set()
    for ring in rings(geometry):
        points = [(round(point[0], 1), round(point[1], 1)) for point in ring]
        vertices.update(points)
        for first, second in zip(points, points[1:]):
            if first != second:
                segments.add(tuple(sorted((first, second))))
    return vertices, segments


def polygon_links(features):
    boundaries = {
        feature["properties"]["name"]: vertices_and_segments(feature["geometry"])
        for feature in features
    }
    links = set()
    for left, right in itertools.combinations(boundaries, 2):
        if boundaries[left][1] & boundaries[right][1]:
            links.add(tuple(sorted((left, right))))
    return links


def morans_i(values_by_name, links):
    names = list(values_by_name)
    centre = mean(values_by_name.values())
    denominator = sum((values_by_name[name] - centre) ** 2 for name in names)
    used_links = [(a, b) for a, b in links if a in values_by_name and b in values_by_name]
    numerator = 2 * sum(
        (values_by_name[a] - centre) * (values_by_name[b] - centre) for a, b in used_links
    )
    return len(names) / (2 * len(used_links)) * numerator / denominator


def row_standardised_morans_i(values_by_name, links):
    names = list(values_by_name)
    centre = mean(values_by_name.values())
    neighbours = {name: set() for name in names}
    for left, right in links:
        if left in neighbours and right in neighbours:
            neighbours[left].add(right)
            neighbours[right].add(left)
    denominator = sum((values_by_name[name] - centre) ** 2 for name in names)
    numerator = sum(
        (values_by_name[name] - centre)
        * sum((values_by_name[other] - centre) / len(neighbours[name]) for other in neighbours[name])
        for name in names
        if neighbours[name]
    )
    return len(names) * numerator / (len(names) * denominator)


def permutation_test(values_by_name, links, iterations=99999, seed=20261010):
    rng = random.Random(seed)
    names = list(values_by_name)
    observed = morans_i(values_by_name, links)
    values = [values_by_name[name] for name in names]
    null = []
    for _ in range(iterations):
        shuffled = values[:]
        rng.shuffle(shuffled)
        null.append(morans_i(dict(zip(names, shuffled)), links))
    p_value = (1 + sum(value >= observed for value in null)) / (iterations + 1)
    return observed, p_value, null


def bootstrap_difference(inner, outer, iterations=10000, seed=2701):
    rng = random.Random(seed)
    draws = []
    for _ in range(iterations):
        a = [rng.choice(inner) for _ in inner]
        b = [rng.choice(outer) for _ in outer]
        draws.append(mean(a) - mean(b))
    return percentile(draws, 0.025), percentile(draws, 0.975)


rows = list(csv.DictReader((ROOT / "data" / "london_housing_2019.csv").open(encoding="utf-8")))
for row in rows:
    row["ratio"] = float(row["price_to_earnings"])
    row["rent"] = int(row["rent_gbp"])
    row["price"] = int(row["house_price_gbp"])

# Deliberately simplified cartogram: west/east and north/south ordering are retained.
XY = {
    "Hillingdon": (0, 3), "Harrow": (2, 1), "Barnet": (4, 0), "Enfield": (6, 0),
    "Haringey": (5, 1), "Waltham Forest": (7, 1), "Redbridge": (8, 2), "Havering": (10, 2),
    "Ealing": (2, 3), "Brent": (3, 2), "Camden": (4, 2), "Islington": (5, 2),
    "Hackney": (6, 2), "Newham": (7, 3), "Barking and Dagenham": (9, 3),
    "Hounslow": (1, 4), "Hammersmith and Fulham": (3, 3), "Kensington and Chelsea": (4, 3),
    "Westminster": (5, 3), "City of London": (6, 3), "Tower Hamlets": (7, 4),
    "Richmond upon Thames": (1, 5), "Wandsworth": (3, 4), "Lambeth": (4, 4),
    "Southwark": (5, 4), "Lewisham": (6, 5), "Greenwich": (7, 5), "Bexley": (9, 5),
    "Kingston upon Thames": (2, 6), "Merton": (3, 5), "Sutton": (3, 7),
    "Croydon": (5, 6), "Bromley": (7, 7),
}
ABBR = {
    "Hillingdon": "HIL", "Harrow": "HAR", "Barnet": "BAR", "Enfield": "ENF",
    "Haringey": "HGY", "Waltham Forest": "WFT", "Redbridge": "RED", "Havering": "HAV",
    "Ealing": "EAL", "Brent": "BRE", "Camden": "CAM", "Islington": "ISL", "Hackney": "HCK",
    "Newham": "NEW", "Barking and Dagenham": "BDG", "Hounslow": "HOU",
    "Hammersmith and Fulham": "HMF", "Kensington and Chelsea": "KNC", "Westminster": "WES",
    "City of London": "COL", "Tower Hamlets": "TWH", "Richmond upon Thames": "RIC",
    "Wandsworth": "WAN", "Lambeth": "LAM", "Southwark": "SOU", "Lewisham": "LEW",
    "Greenwich": "GRE", "Bexley": "BEX", "Kingston upon Thames": "KIN", "Merton": "MER",
    "Sutton": "SUT", "Croydon": "CRO", "Bromley": "BRO",
}

values = {row["borough"]: row["ratio"] for row in rows}
ratios = list(values.values())
inner_values = [row["ratio"] for row in rows if row["sector"] == "Inner"]
outer_values = [row["ratio"] for row in rows if row["sector"] == "Outer"]
median_ratio = statistics.median(ratios)
inner_mean = mean(inner_values)
outer_mean = mean(outer_values)
inner_outer_gap = inner_mean - outer_mean
gap_ci = bootstrap_difference(inner_values, outer_values)

boundary_geo = json.loads((ROOT / "data" / "london_boroughs.geojson").read_text(encoding="utf-8-sig"))
official_links = polygon_links(boundary_geo["features"])
moran, moran_p, permutation_null = permutation_test(values, official_links)
moran_row_standardised = row_standardised_morans_i(values, official_links)

trimmed_values = {name: value for name, value in values.items() if name != "Kensington and Chelsea"}
moran_without_kc, moran_without_kc_p, _ = permutation_test(trimmed_values, official_links)
inner_without_kc = mean(
    row["ratio"] for row in rows if row["sector"] == "Inner" and row["borough"] != "Kensington and Chelsea"
)

ranked = sorted(rows, key=lambda row: row["ratio"], reverse=True)
top_five = [{"borough": row["borough"], "ratio": row["ratio"]} for row in ranked[:5]]
bottom_five = [{"borough": row["borough"], "ratio": row["ratio"]} for row in ranked[-5:]]

# Figure 1 — choropleth on official borough polygons.
W, H = 1500, 950
image = Image.new("RGB", (W, H), "#080b0a")
draw = ImageDraw.Draw(image)
draw.text((70, 55), "LONDON / HOUSING AFFORDABILITY", font=font(20, True), fill="#d7ff50")
draw.text((70, 92), "Where income meets the housing market", font=font(54, True), fill="#f2f0e9")
draw.text((70, 160), "32 boroughs + City of London · official GLA boundary topology", font=font(18), fill="#8e9993")
draw.line((70, 205, 1430, 205), fill="#2c3531", width=2)

def map_colour(value):
    t = max(0, min(1, (value - 9.8) / (24.4 - 9.8)))
    low, high = (44, 86, 71), (215, 255, 80)
    return tuple(int(low[i] + t * (high[i] - low[i])) for i in range(3))

def geometry_polygons(geometry):
    return geometry["coordinates"] if geometry["type"] == "MultiPolygon" else [geometry["coordinates"]]


all_points = [
    point
    for feature in boundary_geo["features"]
    for polygon in geometry_polygons(feature["geometry"])
    for ring in polygon
    for point in ring
]
xmin, xmax = min(point[0] for point in all_points), max(point[0] for point in all_points)
ymin, ymax = min(point[1] for point in all_points), max(point[1] for point in all_points)
left, top, right, bottom = 70, 240, 1030, 875
scale = min((right - left) / (xmax - xmin), (bottom - top) / (ymax - ymin))
map_width, map_height = (xmax - xmin) * scale, (ymax - ymin) * scale
map_left = left + (right - left - map_width) / 2
map_top = top + (bottom - top - map_height) / 2


def screen(point):
    return map_left + (point[0] - xmin) * scale, map_top + (ymax - point[1]) * scale


for feature in boundary_geo["features"]:
    name = feature["properties"]["name"]
    fill = CORAL if name == "Kensington and Chelsea" else map_colour(values[name])
    polygons = geometry_polygons(feature["geometry"])
    for polygon in polygons:
        draw.polygon([screen(point) for point in polygon[0]], fill=fill, outline="#080b0a")
        for hole in polygon[1:]:
            draw.polygon([screen(point) for point in hole], fill="#080b0a")
    # Abbreviate the five most extreme areas to keep the map legible.
    if name in {item["borough"] for item in top_five}:
        largest = max((polygon[0] for polygon in polygons), key=len)
        cx = statistics.mean(point[0] for point in largest)
        cy = statistics.mean(point[1] for point in largest)
        sx, sy = screen((cx, cy))
        label = ABBR[name]
        box = draw.textbbox((0, 0), label, font=font(11, True))
        draw.rectangle((sx - 4, sy - 3, sx + box[2] - box[0] + 5, sy + 14), fill="#080b0a")
        draw.text((sx, sy - 2), label, font=font(11, True), fill="#f2f0e9")

draw.rounded_rectangle((1080, 260, 1425, 700), radius=18, fill="#111714", outline="#2c3531", width=2)
draw.text((1120, 300), "DIAGNOSTIC READOUT", font=font(12, True), fill="#d7ff50")
draw.text((1120, 350), f"{median_ratio:.1f}×", font=font(62, True), fill="#f2f0e9")
draw.text((1120, 420), "borough median", font=font(16), fill="#8e9993")
draw.text((1120, 475), f"Inner mean     {inner_mean:.1f}×", font=font(18, True), fill="#f2f0e9")
draw.text((1120, 515), f"Outer mean     {outer_mean:.1f}×", font=font(18, True), fill="#f2f0e9")
draw.text((1120, 565), f"Moran's I      {moran:.2f}", font=font(18, True), fill="#d7ff50")
draw.text((1120, 605), f"Permutation p  {moran_p:.5f}", font=font(18, True), fill="#d7ff50")
draw.text((1120, 655), "Exploratory · not causal", font=font(14), fill="#8e9993")
draw.text((70, 885), "METHOD  ·  82 shared-boundary links  /  99,999 permutations  /  outlier sensitivity", font=font(13, True), fill="#8e9993")
image.save(OUT / "london_affordability.png", optimize=True)

# Figure 2 — full ranked distribution.
W, H = 1600, 1320
image = Image.new("RGB", (W, H), PAPER)
draw = ImageDraw.Draw(image)
draw.text((78, 58), "FIGURE 02 · DISTRIBUTION", font=font(16, True), fill=COBALT)
draw.text((78, 98), "Affordability pressure is not a single London average", font=font(44, serif=True), fill=INK)
draw.text((80, 166), "Published price-to-workplace-earnings ratios for all 33 local areas. The vertical marker is the borough median.", font=font(18), fill=SOFT)
draw.line((80, 218, 1520, 218), fill=INK, width=2)

chart_left, chart_right = 425, 1490
chart_top, row_height = 252, 29
max_value = 46
for tick in (0, 10, 20, 30, 40):
    x = chart_left + tick / max_value * (chart_right - chart_left)
    draw.line((x, chart_top - 8, x, chart_top + row_height * len(ranked)), fill=HAIR, width=1)
    draw.text((x - 8, chart_top - 34), str(tick), font=font(13), fill=SOFT)

median_x = chart_left + median_ratio / max_value * (chart_right - chart_left)
draw.line((median_x, chart_top - 8, median_x, chart_top + row_height * len(ranked)), fill=COBALT, width=3)
draw.text((median_x + 8, chart_top - 34), f"median {median_ratio:.1f}×", font=font(13, True), fill=COBALT)

for index, row in enumerate(ranked):
    y = chart_top + index * row_height
    draw.text((80, y - 2), f"{index + 1:02d}", font=font(12), fill=SOFT)
    draw.text((122, y - 2), row["borough"], font=font(14), fill=INK)
    x = chart_left + row["ratio"] / max_value * (chart_right - chart_left)
    colour = CORAL if row["borough"] == "Kensington and Chelsea" else (COBALT if row["sector"] == "Inner" else INK)
    draw.line((chart_left, y + 8, x, y + 8), fill=colour, width=5)
    draw.ellipse((x - 6, y + 2, x + 6, y + 14), fill=colour)
    draw.text((x + 13, y - 2), f"{row['ratio']:.1f}×", font=font(13, True), fill=colour)

legend_y = 1240
draw.ellipse((80, legend_y, 92, legend_y + 12), fill=COBALT)
draw.text((102, legend_y - 4), "Inner London", font=font(14), fill=SOFT)
draw.ellipse((255, legend_y, 267, legend_y + 12), fill=INK)
draw.text((277, legend_y - 4), "Outer London", font=font(14), fill=SOFT)
draw.ellipse((445, legend_y, 457, legend_y + 12), fill=CORAL)
draw.text((467, legend_y - 4), "Kensington & Chelsea outlier", font=font(14), fill=SOFT)
draw.text((80, 1282), "Source: GLA, Housing in London 2019, Table 3. Ratios are reported values, not model predictions.", font=font(13), fill=SOFT)
image.save(OUT / "london_ranked_boroughs.png", optimize=True)

# Figure 3 — robustness and inference diagnostics.
W, H = 1600, 920
image = Image.new("RGB", (W, H), PAPER)
draw = ImageDraw.Draw(image)
draw.text((78, 58), "FIGURE 03 · ROBUSTNESS", font=font(16, True), fill=COBALT)
draw.text((78, 98), "What survives alternative analytical choices?", font=font(44, serif=True), fill=INK)
draw.text((80, 166), "Permutation inference, weight standardisation and an explicit outlier check.", font=font(18), fill=SOFT)
draw.line((80, 218, 1520, 218), fill=INK, width=2)

draw.text((80, 252), "A · SPATIAL RANDOMISATION", font=font(14, True), fill=COBALT)
draw.text((80, 286), "Moran's I under 99,999 random reallocations", font=font(27, serif=True), fill=INK)
panel_left, panel_top, panel_right, panel_bottom = 80, 350, 940, 690
lo, hi, bins = min(permutation_null), max(permutation_null + [moran]), 36
counts = [0] * bins
for value in permutation_null:
    bucket = min(bins - 1, int((value - lo) / (hi - lo) * bins))
    counts[bucket] += 1
max_count = max(counts)
for i, count in enumerate(counts):
    x1 = panel_left + i / bins * (panel_right - panel_left)
    x2 = panel_left + (i + 1) / bins * (panel_right - panel_left) - 2
    y1 = panel_bottom - count / max_count * (panel_bottom - panel_top)
    draw.rectangle((x1, y1, x2, panel_bottom), fill=PALE_BLUE)
observed_x = panel_left + (moran - lo) / (hi - lo) * (panel_right - panel_left)
draw.line((observed_x, panel_top - 10, observed_x, panel_bottom + 4), fill=CORAL, width=5)
draw.text((observed_x - 92, panel_top - 42), f"observed I = {moran:.3f}", font=font(14, True), fill=CORAL)
draw.line((panel_left, panel_bottom, panel_right, panel_bottom), fill=INK, width=2)
draw.text((panel_left, panel_bottom + 16), f"null min {lo:.2f}", font=font(12), fill=SOFT)
draw.text((panel_right - 95, panel_bottom + 16), f"max {hi:.2f}", font=font(12), fill=SOFT)
draw.text((80, 752), f"One-sided permutation p = {moran_p:.5f}. Expected I under randomisation = {-1/(len(rows)-1):.3f}.", font=font(16), fill=INK)

x0 = 1030
draw.text((x0, 252), "B · SPECIFICATION CURVE", font=font(14, True), fill=COBALT)
draw.text((x0, 286), "Moran's I by weight rule", font=font(27, serif=True), fill=INK)
trimmed_link_count = sum(a in trimmed_values and b in trimmed_values for a, b in official_links)
items = [
    ("Binary contiguity", len(official_links), moran),
    ("Row-standardised", len(official_links), moran_row_standardised),
    ("Without K&C", trimmed_link_count, moran_without_kc),
]
for i, (label, link_count, statistic) in enumerate(items):
    y = 370 + i * 105
    draw.text((x0, y), label, font=font(15), fill=INK)
    draw.text((x0, y + 30), f"{link_count} usable links", font=font(12), fill=SOFT)
    x1, x2 = x0 + 230, 1490
    draw.line((x1, y + 14, x2, y + 14), fill=HAIR, width=3)
    value_x = x1 + statistic / 0.4 * (x2 - x1)
    draw.ellipse((value_x - 8, y + 6, value_x + 8, y + 22), fill=COBALT)
    draw.text((value_x - 20, y + 31), f"{statistic:.3f}", font=font(13, True), fill=COBALT)

draw.line((x0, 705, 1510, 705), fill=HAIR, width=1)
draw.text((x0, 730), "OUTLIER CHECK", font=font(12, True), fill=COBALT)
draw.text((x0, 762), f"Inner mean: {inner_mean:.1f}× → {inner_without_kc:.1f}× without K&C", font=font(15), fill=INK)
draw.text((x0, 792), f"Moran's I: {moran:.3f} → {moran_without_kc:.3f} without K&C", font=font(15), fill=INK)
draw.text((80, 864), "Interpretation: high-ratio areas are more adjacent than under random reallocation on the official 82-link borough graph.", font=font(14), fill=SOFT)
image.save(OUT / "london_spatial_diagnostics.png", optimize=True)

summary = {
    "status": "exploratory spatial analysis; not causal",
    "sources": [
        "GLA Housing in London 2019, Table 3",
        "GLA London Boroughs dataset e55pw / MapServer layer 3",
    ],
    "local_areas": len(rows),
    "official_london_aggregate_ratio": 12.3,
    "median_ratio": round(median_ratio, 2),
    "iqr": [round(percentile(ratios, 0.25), 2), round(percentile(ratios, 0.75), 2)],
    "inner_mean": round(inner_mean, 2),
    "outer_mean": round(outer_mean, 2),
    "inner_outer_gap": round(inner_outer_gap, 2),
    "inner_outer_mean_gap_iid_resampling_interval_95pct": [round(gap_ci[0], 2), round(gap_ci[1], 2)],
    "inner_median": round(statistics.median(inner_values), 2),
    "outer_median": round(statistics.median(outer_values), 2),
    "inner_mean_without_kensington_chelsea": round(inner_without_kc, 2),
    "morans_i": round(moran, 4),
    "official_shared_boundary_links": len(official_links),
    "morans_i_row_standardised": round(moran_row_standardised, 5),
    "moran_permutation_p_one_sided": round(moran_p, 6),
    "permutations": 99999,
    "expected_i": round(-1 / (len(rows) - 1), 4),
    "moran_without_kensington_chelsea": round(moran_without_kc, 4),
    "moran_without_kensington_chelsea_p_one_sided": round(moran_without_kc_p, 6),
    "top_five": top_five,
    "bottom_five": bottom_five,
    "limitations": [
        "Binary shared-boundary weights ignore distance, flows and cross-boundary intensity.",
        "The 2018 ratio cross-section cannot identify causal effects or dynamics.",
        "Borough averages conceal within-borough heterogeneity and tenure differences.",
        "Workplace-based earnings need not represent the earnings of resident homebuyers.",
    ],
}
(ROOT / "london_summary.json").write_text(json.dumps(summary, indent=2), encoding="utf-8")
print(json.dumps(summary, indent=2))