perfect-postcode/analysis/cheaper_twins.py

#!/usr/bin/env python3
"""Cheaper-twin / name-premium index over all-England property data.

This is the ROOT growth artifact: every page, OG card, video and outreach number derives from its
output. It reads the local property data, aggregates to POSTCODE SECTOR grain (e.g. "N1 1", the same
grain the homepage TwinProof block uses, which is load-bearing), and finds "cheaper twins": nearby
sectors that match on property type, build era, school provision and station access but differ
materially in estimated price per square metre.

Defensibility rules baked in (see growth/README.md):
  - Sector aggregation only, never address-level output (Royal Mail / OS rights).
  - Minimum sample sizes per sector (--min-props, --min-recorded).
  - Robust statistics: median £/sqm; a sector needs real recorded sales, not just modelled estimates.
  - England only (ctry25cd starts with "E").
  - Every row stamps its N.
  - "Twin" requires genuine like-for-like matching before any price claim is made.

Outputs (analysis/out/):
  - sector_index.parquet / .csv:  per-sector value table (powers "best value", "£100k buys X m²", etc.)
  - cheaper_twins.parquet / .csv:  ranked twin pairs (pricey name -> cheaper twin)
  - national_facts.json:           headline stats for collateral/finding placeholders

Run:  source .venv/bin/activate && python analysis/cheaper_twins.py
"""
from __future__ import annotations

import argparse
import json
import math
from pathlib import Path

import numpy as np
import polars as pl
from scipy.spatial import cKDTree

DATA = Path("property-data")
OUT = Path("analysis/out")
# Postcode sector = outward code + first inward digit, e.g. "N1 1", "M20 2".
SECTOR_RE = r"^([A-Z]{1,2}[0-9][A-Z0-9]? [0-9])"

STATION_DIST_COLS = [
    "Distance to nearest amenity (Rail station) (km)",
    "Distance to nearest amenity (Tube station) (km)",
    "Distance to nearest amenity (DLR station) (km)",
    "Distance to nearest amenity (Tram & Metro stop) (km)",
]
AREA_COLS = [
    "lat",
    "lon",
    "Good+ primary school catchments",
    "Good+ secondary school catchments",
    "Outstanding primary school catchments",
    "Outstanding secondary school catchments",
    "Serious crime (/yr, 7y)",
    "Minor crime (/yr, 7y)",
    "Noise (dB)",
    "Max available download speed (Mbps)",
    "Median age",
    "% Owner occupied",
    "% Degree or higher",
    *STATION_DIST_COLS,
]


def _collect(lf: pl.LazyFrame) -> pl.DataFrame:
    """Collect with streaming if the installed polars supports it, else fall back."""
    try:
        return lf.collect(streaming=True)
    except Exception:
        return lf.collect()


def property_aggregates() -> pl.DataFrame:
    props = (
        pl.scan_parquet(DATA / "properties.parquet")
        .with_columns(pl.col("Postcode").str.extract(SECTOR_RE, 1).alias("Sector"))
        .filter(pl.col("Sector").is_not_null())
    )
    agg = props.group_by("Sector").agg(
        pl.len().alias("n_props"),
        pl.col("Price per sqm").drop_nulls().len().alias("n_recorded"),
        pl.col("Est. price per sqm").median().alias("est_psqm"),
        pl.col("Price per sqm").median().alias("recorded_psqm"),
        pl.col("Last known price").median().alias("median_price"),
        pl.col("Estimated current price").median().alias("est_price"),
        pl.col("Total floor area (sqm)").median().alias("median_floor"),
        pl.col("Number of bedrooms & living rooms").median().alias("median_rooms"),
        pl.col("Construction year")
        .filter(pl.col("Construction year") > 1800)
        .median()
        .alias("median_build_year"),
    )

    # Dominant property type + its share, per sector (robust mode without ordering assumptions).
    tc = props.group_by(["Sector", "Property type"]).agg(pl.len().alias("c"))
    tc = tc.with_columns(
        (pl.col("c") / pl.col("c").sum().over("Sector")).alias("share"),
        pl.col("c").max().over("Sector").alias("maxc"),
    )
    dom = (
        tc.filter(pl.col("c") == pl.col("maxc"))
        .unique(subset="Sector", keep="first")
        .select(
            pl.col("Sector"),
            pl.col("Property type").alias("dominant_type"),
            pl.col("share").alias("dominant_share"),
        )
    )
    return _collect(agg).join(_collect(dom), on="Sector", how="left")


def area_aggregates() -> pl.DataFrame:
    # Property counts per postcode unit -> weights, so area features are weighted by housing stock.
    unit_counts = _collect(
        pl.scan_parquet(DATA / "properties.parquet")
        .group_by("Postcode")
        .agg(pl.len().alias("n"))
    )

    pc = pl.read_parquet(DATA / "postcode.parquet", columns=["Postcode", "ctry25cd", *AREA_COLS])
    pc = pc.filter(pl.col("ctry25cd").str.starts_with("E"))  # England only
    pc = pc.join(unit_counts, on="Postcode", how="inner")  # only units that contain homes
    pc = pc.with_columns(
        pl.col("Postcode").str.extract(SECTOR_RE, 1).alias("Sector"),
        pl.min_horizontal(STATION_DIST_COLS).alias("dist_station_km"),
    ).filter(pl.col("Sector").is_not_null())

    def wmean(col: str, alias: str) -> pl.Expr:
        # weighted mean ignoring nulls: sum(col*n) / sum(n where col not null)
        num = (pl.col(col) * pl.col("n")).sum()
        den = pl.col("n").filter(pl.col(col).is_not_null()).sum()
        return (num / den).alias(alias)

    return pc.group_by("Sector").agg(
        pl.col("n").sum().alias("area_n"),
        wmean("lat", "lat"),
        wmean("lon", "lon"),
        wmean("Good+ primary school catchments", "good_primary"),
        wmean("Good+ secondary school catchments", "good_secondary"),
        wmean("Outstanding primary school catchments", "outstanding_primary"),
        wmean("Outstanding secondary school catchments", "outstanding_secondary"),
        wmean("Serious crime (/yr, 7y)", "serious_crime"),
        wmean("Noise (dB)", "noise_db"),
        wmean("Max available download speed (Mbps)", "broadband_mbps"),
        wmean("Median age", "median_age"),
        wmean("% Owner occupied", "pct_owner"),
        wmean("% Degree or higher", "pct_degree"),
        wmean("dist_station_km", "dist_station_km"),
    )


def build_index(args) -> pl.DataFrame:
    idx = property_aggregates().join(area_aggregates(), on="Sector", how="inner")
    idx = idx.filter(
        (pl.col("n_props") >= args.min_props)
        & (pl.col("n_recorded") >= args.min_recorded)
        & pl.col("est_psqm").is_not_null()
        & (pl.col("est_psqm") > 0)
        & pl.col("lat").is_not_null()
        & pl.col("lon").is_not_null()
        & pl.col("dist_station_km").is_not_null()
    )
    idx = idx.with_columns(
        (100_000 / pl.col("est_psqm")).round(1).alias("sqm_per_100k"),
        pl.col("est_psqm").round().cast(pl.Int64),
        pl.col("recorded_psqm").round().cast(pl.Int64),
    )
    return idx.sort("est_psqm", descending=True)


def find_twins(idx: pl.DataFrame, args) -> pl.DataFrame:
    d = idx.to_dict(as_series=False)
    n = len(d["Sector"])
    lat = np.array(d["lat"], float)
    lon = np.array(d["lon"], float)
    psqm = np.array(d["est_psqm"], float)
    floor = np.array([f if f is not None else np.nan for f in d["median_floor"]], float)
    build = np.array([b if b is not None else np.nan for b in d["median_build_year"]], float)
    dom = d["dominant_type"]
    good_sec = np.array(d["good_secondary"], float)
    good_pri = np.array(d["good_primary"], float)
    crime = np.array(d["serious_crime"], float)
    station = np.array(d["dist_station_km"], float)
    owner = np.array(d["pct_owner"], float)
    degree = np.array(d["pct_degree"], float)
    age = np.array(d["median_age"], float)

    # Planar projection (km) for a local KD-tree radius search.
    lat0 = math.radians(float(np.nanmean(lat)))
    x = lon * 111.320 * math.cos(lat0)
    y = lat * 110.574
    tree = cKDTree(np.column_stack([x, y]))

    rows = []
    for i in range(n):
        neigh = tree.query_ball_point([x[i], y[i]], args.max_km)
        best_j, best_gap = -1, -1.0
        for j in neigh:
            if j == i or psqm[i] <= psqm[j]:
                continue  # i must be the pricier side
            gap = 1.0 - psqm[j] / psqm[i]
            # A genuine "twin" sits in a believable band: below min it's not a story,
            # above max it's a different market tier (city-centre premium / prime), not a twin.
            if gap < args.min_gap or gap > args.max_gap:
                continue
            if dom[i] != dom[j]:
                continue
            if not (np.isfinite(build[i]) and np.isfinite(build[j])) or abs(build[i] - build[j]) > args.build_band:
                continue
            if abs(good_sec[i] - good_sec[j]) > args.school_tol or abs(good_pri[i] - good_pri[j]) > args.school_tol:
                continue
            if station[i] > args.station_max or station[j] > args.station_max or abs(station[i] - station[j]) > args.station_tol:
                continue
            # Similarity gates: the two must be the SAME KIND of neighbourhood, so the gap
            # reads as a name premium, not a tier jump (deprivation/tenure/education/age/safety/size).
            if crime[j] > crime[i] * args.crime_ratio or crime[i] > crime[j] * args.crime_ratio:
                continue
            if abs(owner[i] - owner[j]) > args.owner_tol:
                continue
            if abs(degree[i] - degree[j]) > args.degree_tol:
                continue
            if np.isfinite(age[i]) and np.isfinite(age[j]) and abs(age[i] - age[j]) > args.age_tol:
                continue
            if np.isfinite(floor[i]) and np.isfinite(floor[j]) and floor[j] > 0:
                fr = floor[i] / floor[j]
                if fr < args.floor_ratio or fr > 1.0 / args.floor_ratio:
                    continue
            if gap > best_gap:
                best_gap, best_j = gap, j
        if best_j < 0:
            continue
        j = best_j
        avg_floor = np.nanmean([floor[i], floor[j]])
        if not np.isfinite(avg_floor):
            avg_floor = 90.0
        rows.append(
            {
                "pricey_sector": d["Sector"][i],
                "twin_sector": d["Sector"][j],
                "pricey_psqm": int(psqm[i]),
                "twin_psqm": int(psqm[j]),
                "gap_pct": round(best_gap * 100, 1),
                "gap_per_sqm": int(psqm[i] - psqm[j]),
                "gap_on_avg_home": int((psqm[i] - psqm[j]) * avg_floor),
                "gap_on_90sqm": int((psqm[i] - psqm[j]) * 90),
                "dist_km": round(math.hypot(x[i] - x[j], y[i] - y[j]), 2),
                "dominant_type": dom[i],
                "build_year": None if not np.isfinite(build[i]) else int(build[i]),
                "good_secondary": round(float(good_sec[i]), 1),
                "station_km": round(float(station[i]), 2),
                "pricey_lat": round(float(lat[i]), 5),
                "pricey_lon": round(float(lon[i]), 5),
                "twin_lat": round(float(lat[j]), 5),
                "twin_lon": round(float(lon[j]), 5),
                "pricey_n": int(d["n_props"][i]),
                "twin_n": int(d["n_props"][j]),
            }
        )

    if not rows:
        return pl.DataFrame()
    tw = pl.DataFrame(rows)
    # Dedup unordered pairs, keep the biggest gap.
    tw = tw.with_columns(
        pl.concat_list(
            [
                pl.min_horizontal("pricey_sector", "twin_sector"),
                pl.max_horizontal("pricey_sector", "twin_sector"),
            ]
        )
        .list.join("|")
        .alias("_pair")
    )
    tw = tw.sort("gap_pct", descending=True).unique(subset="_pair", keep="first").drop("_pair")
    return tw.filter(pl.col("gap_on_90sqm") >= args.min_abs_gap).sort("gap_pct", descending=True)


def national_facts(idx: pl.DataFrame, twins: pl.DataFrame, args) -> dict:
    valid = idx.filter(pl.col("n_props") >= args.min_props)
    cheapest = valid.sort("est_psqm").head(1).to_dicts()[0]
    dearest = valid.sort("est_psqm", descending=True).head(1).to_dicts()[0]
    facts = {
        "generated_with": "analysis/cheaper_twins.py",
        "params": vars(args),
        "n_sectors": idx.height,
        "n_twin_pairs": twins.height,
        "attribution": "Contains HM Land Registry data © Crown copyright and database right. Licensed under the Open Government Licence v3.0.",
        "best_value_sector": {"sector": cheapest["Sector"], "est_psqm": cheapest["est_psqm"], "sqm_per_100k": cheapest["sqm_per_100k"], "n": cheapest["n_props"]},
        "dearest_sector": {"sector": dearest["Sector"], "est_psqm": dearest["est_psqm"], "sqm_per_100k": dearest["sqm_per_100k"], "n": dearest["n_props"]},
    }
    if twins.height:
        top = twins.head(10).to_dicts()
        facts["biggest_twin_gap"] = top[0]
        facts["top_twins"] = top
    return facts


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument("--min-props", type=int, default=150, help="min properties per sector")
    p.add_argument("--min-recorded", type=int, default=40, help="min recorded sales (with floor area) per sector")
    p.add_argument("--max-km", type=float, default=3.0, help="max centroid distance for a twin (km)")
    p.add_argument("--min-gap", type=float, default=0.15, help="min fractional £/sqm gap")
    p.add_argument("--max-gap", type=float, default=0.45, help="max fractional gap (above this it's a tier jump, not a twin)")
    p.add_argument("--build-band", type=float, default=30, help="max build-year difference")
    p.add_argument("--school-tol", type=float, default=1.5, help="max difference in good-catchment counts")
    p.add_argument("--station-max", type=float, default=1.5, help="both sectors must be within this many km of a station")
    p.add_argument("--station-tol", type=float, default=0.9, help="max difference in station distance (km)")
    p.add_argument("--crime-ratio", type=float, default=1.5, help="serious crime must be within this ratio either way")
    p.add_argument("--owner-tol", type=float, default=22, help="max difference in %% owner-occupied")
    p.add_argument("--degree-tol", type=float, default=22, help="max difference in %% degree-or-higher")
    p.add_argument("--age-tol", type=float, default=12, help="max difference in median age (years)")
    p.add_argument("--floor-ratio", type=float, default=0.72, help="median floor area must be within this ratio either way")
    p.add_argument("--min-abs-gap", type=int, default=20000, help="min £ gap on a 90 sqm home for the twin list")
    args = p.parse_args()

    OUT.mkdir(parents=True, exist_ok=True)
    print("Aggregating 22.4M properties to sector grain ...")
    idx = build_index(args)
    print(f"  {idx.height} valid England sectors (>= {args.min_props} props, >= {args.min_recorded} recorded sales)")
    idx.write_parquet(OUT / "sector_index.parquet")
    idx.write_csv(OUT / "sector_index.csv")

    print("Matching cheaper twins ...")
    twins = find_twins(idx, args)
    print(f"  {twins.height} twin pairs")
    if twins.height:
        twins.write_parquet(OUT / "cheaper_twins.parquet")
        twins.write_csv(OUT / "cheaper_twins.csv")

    facts = national_facts(idx, twins, args)
    (OUT / "national_facts.json").write_text(json.dumps(facts, indent=2, default=str))
    print(f"Wrote outputs to {OUT}/")
    if twins.height:
        print("\nTop 10 cheaper twins (pricey -> twin, gap%, £ on 90sqm):")
        for r in twins.head(10).to_dicts():
            print(f"  {r['pricey_sector']:>7} -> {r['twin_sector']:<7}  {r['gap_pct']:>5}%  £{r['gap_on_90sqm']:,}  ({r['dominant_type']}, ~{r['build_year']})")


if __name__ == "__main__":
    main()