#!/usr/bin/env python3
"""
HCD APR Table A2 — Ownership x Affordability analysis
======================================================
Question: Of all housing units permitted in California (statewide and in
high-growth counties), what share are BOTH owner-tenure AND affordable
(very-low or low income)? Expressed as "1 in N".

Source (state's own data, self-reported by every jurisdiction, updated weekly):
  https://data.ca.gov/dataset/housing-element-annual-progress-report-apr-data-by-jurisdiction-and-year
  Table A2 CSV (~189 MB):
  https://data.ca.gov/dataset/81b0841f-2802-403e-b48e-2ef4b751f77c/resource/fe505d9b-8c36-42ba-ba30-08bc4f34e022/download/tablea2.csv

Run:  pip install pandas requests && python3 apr_ownership_affordability.py

Design notes:
- Stage = BUILDING PERMITS (the RHNA benchmark HCD itself uses), not
  entitlements or certificates of occupancy.
- Two variants computed, per the methodology discussion:
    STRICT  = deed-restricted (DR) VLI+LI units only
    BROAD   = DR + non-deed-restricted (NDR) VLI+LI units
  NDR includes "assumed affordability" reporting (common for ADUs). Quote
  STRICT as the headline; if BROAD diverges a lot, that gap is itself a
  finding (paper affordability vs. enforceable affordability).
- ADUs broken out separately: jurisdictions frequently count ADUs as
  affordable via assumed-affordability memos, and ADUs are not
  family-purchasable ownership product. Check the with/without-ADU numbers
  before publishing.
- Column names in the raw CSV drift across HCD refreshes, so this script
  DISCOVERS columns by pattern and prints its mapping for you to verify
  before trusting output. If detection fails, it dumps all columns and
  exits — adjust the patterns and rerun.
"""

import re
import sys
from pathlib import Path

import pandas as pd

CSV_URL = ("https://data.ca.gov/dataset/81b0841f-2802-403e-b48e-2ef4b751f77c/"
           "resource/fe505d9b-8c36-42ba-ba30-08bc4f34e022/download/tablea2.csv")
LOCAL = Path(__file__).parent / "tablea2.csv"

HIGH_GROWTH_COUNTIES = {
    # SoCal metro
    "LOS ANGELES", "ORANGE", "SAN DIEGO",
    # Bay Area 9
    "ALAMEDA", "CONTRA COSTA", "MARIN", "NAPA", "SAN FRANCISCO",
    "SAN MATEO", "SANTA CLARA", "SOLANO", "SONOMA",
}


def download():
    if LOCAL.exists() and LOCAL.stat().st_size > 100_000_000:
        print(f"Using cached {LOCAL} ({LOCAL.stat().st_size/1e6:.0f} MB)")
        return
    import requests
    print("Downloading Table A2 (~189 MB)...")
    with requests.get(CSV_URL, stream=True, timeout=600) as r:
        r.raise_for_status()
        with open(LOCAL, "wb") as f:
            for chunk in r.iter_content(chunk_size=1 << 20):
                f.write(chunk)
    print("Done.")


def find_col(cols, *patterns, required=True):
    """Return the first column matching all regex patterns (case-insensitive)."""
    for c in cols:
        if all(re.search(p, c, re.I) for p in patterns):
            return c
    if required:
        print("\nFAILED to find a column matching:", patterns)
        print("All columns:\n", "\n ".join(cols))
        sys.exit(1)
    return None


def main():
    download()
    df = pd.read_csv(LOCAL, low_memory=False)
    cols = list(df.columns)
    print(f"\nLoaded {len(df):,} rows, {len(cols)} columns.")

    # ---- Column discovery (verify this mapping in the printout) ----------
    county = find_col(cols, r"cnty|county")
    year = find_col(cols, r"year")
    tenure = find_col(cols, r"tenure")
    unit_cat = find_col(cols, r"unit.?cat", required=False)  # ADU/SFD/2-4/5+/MH

    # Building-permit-stage income columns. Per the HCD Table A2 data
    # dictionary, HCD reports SIX non-overlapping income bands at each
    # stage (confirmed against the actual dictionary text, not guessed):
    #   BP_ACUTELY_LOW_INCOME_{DR,NDR}    0-15% AMI
    #   BP_EXTREMELY_LOW_INCOME_{DR,NDR}  15-30% AMI
    #   BP_VLOW_INCOME_{DR,NDR}           30-50% AMI
    #   BP_LOW_INCOME_{DR,NDR}            50-80% AMI
    #   BP_MOD_INCOME_{DR,NDR}            moderate
    #   BP_ABOVE_MOD_INCOME               above moderate
    # Acutely-low and extremely-low are ADDITIVE bands, not subsets of
    # very-low — an earlier draft of this script missed them from the
    # total entirely and also mismatched "acutely low" into the "low
    # income" slot via a loose regex. "Very Low Income" (VLI) as commonly
    # used in RHNA reporting = acutely + extremely + vlow (0-50% AMI).
    # "Low Income" (LI) = the 50-80% AMI band alone.
    acutely_dr = find_col(cols, r"acutely", r"bp|permit", r"(^|_)dr$|deed")
    acutely_ndr = find_col(cols, r"acutely", r"bp|permit", r"ndr|non")
    extremely_dr = find_col(cols, r"extremely", r"bp|permit", r"(^|_)dr$|deed")
    extremely_ndr = find_col(cols, r"extremely", r"bp|permit", r"ndr|non")
    vlow_dr = find_col(cols, r"v.?low", r"bp|permit", r"(^|_)dr$|deed")
    vlow_ndr = find_col(cols, r"v.?low", r"bp|permit", r"ndr|non")
    li_dr = find_col(cols, r"(?<!v)(?<!extremely.)(?<!acutely.)low", r"bp|permit",
                     r"(^|_)dr$|deed")
    li_ndr = find_col(cols, r"(?<!v)(?<!extremely.)(?<!acutely.)low", r"bp|permit",
                      r"ndr|non")
    mod_dr = find_col(cols, r"(?<!above.)mod", r"bp|permit", r"(^|_)dr$|deed",
                      required=False)
    mod_ndr = find_col(cols, r"(?<!above.)mod", r"bp|permit", r"ndr|non",
                       required=False)
    amod = find_col(cols, r"above", r"bp|permit")

    vli_dr = [acutely_dr, extremely_dr, vlow_dr]
    vli_ndr = [acutely_ndr, extremely_ndr, vlow_ndr]

    income_cols = [c for c in (acutely_dr, acutely_ndr, extremely_dr,
                               extremely_ndr, vlow_dr, vlow_ndr, li_dr, li_ndr,
                               mod_dr, mod_ndr, amod) if c]
    print("\nColumn mapping (VERIFY THIS):")
    for name, c in [("county", county), ("year", year), ("tenure", tenure),
                    ("unit_category", unit_cat),
                    ("Acutely-low DR (0-15% AMI)", acutely_dr),
                    ("Acutely-low non-DR", acutely_ndr),
                    ("Extremely-low DR (15-30% AMI)", extremely_dr),
                    ("Extremely-low non-DR", extremely_ndr),
                    ("Very-low DR (30-50% AMI)", vlow_dr),
                    ("Very-low non-DR", vlow_ndr),
                    ("Low-income DR (50-80% AMI)", li_dr),
                    ("Low-income non-DR", li_ndr),
                    ("Mod DR", mod_dr), ("Mod non-DR", mod_ndr),
                    ("Above-moderate", amod)]:
        print(f"  {name:32s} -> {c}")

    for c in income_cols:
        df[c] = pd.to_numeric(df[c], errors="coerce").fillna(0)

    df["bp_units"] = df[income_cols].sum(axis=1)
    df = df[df["bp_units"] > 0].copy()          # permit-stage rows only

    # ---- Dedup: jurisdictions re-list still-incomplete permits in every
    # subsequent APR until certificate of occupancy, so the same physical
    # permit (same tracking ID/APN/BP issue date/unit category) shows up
    # under multiple different APR YEAR values. Verified against raw rows
    # (e.g. a Big Bear Lake SFD permit issued 2016-02-17 appears, unit for
    # unit, in the 2018/2019/2020/2021 APRs). Collapse to one row per
    # permit and attribute it to the year it was actually issued
    # (BP_ISSUE_DT1), not the APR reporting year, which is unreliable for
    # repeat-listed rows.
    bp_issue_col = find_col(cols, r"bp.*issue.*dt|bp.*dt.*issue", required=False)
    tracking_col = find_col(cols, r"track.*id|jurs.*track", required=False)
    apn_col = find_col(cols, r"^apn$", required=False) or find_col(cols, r"apn", required=False)
    if bp_issue_col and tracking_col and apn_col:
        bp_issue_year = pd.to_datetime(df[bp_issue_col], errors="coerce").dt.year
        dedup_key = (df[tracking_col].astype(str) + "|" + df[apn_col].astype(str) + "|"
                     + df[bp_issue_col].astype(str) + "|" + df[unit_cat].astype(str))
        before_rows, before_units = len(df), df["bp_units"].sum()
        df = df.assign(_dedup_key=dedup_key, _bp_issue_year=bp_issue_year)
        df = df.sort_values(bp_issue_col, na_position="last")
        df = df.drop_duplicates(subset="_dedup_key", keep="first").copy()
        after_rows, after_units = len(df), df["bp_units"].sum()
        print(f"\n--- Dedup by permit identity (tracking ID + APN + BP issue date + unit category) ---")
        print(f"Rows: {before_rows:,} -> {after_rows:,}  ({before_rows - after_rows:,} dropped)")
        print(f"Units: {before_units:,.0f} -> {after_units:,.0f}  "
              f"({before_units - after_units:,.0f} dropped, "
              f"{(before_units - after_units) / before_units:.1%})")
        # Reassign the analysis year to the actual BP issue year where parseable;
        # fall back to the original APR reporting YEAR otherwise.
        df[year] = df["_bp_issue_year"].fillna(df[year]).astype(int)
    else:
        print("\nWARNING: could not find tracking ID / APN / BP issue date columns "
              "for dedup — proceeding WITHOUT deduplication. Totals will be inflated "
              "by repeat-year reporting of the same permits.")

    # ---- Tenure: report missing/blank before coercing to owner/renter ----
    tenure_raw = df[tenure].astype(str).str.strip().str.upper()
    is_blank_tenure = tenure_raw.isin(["", "NAN", "NONE"]) | df[tenure].isna()
    df["_tenure_missing"] = is_blank_tenure
    df["is_owner"] = tenure_raw.str.startswith("O") & ~is_blank_tenure

    print("\n--- Missing/blank tenure by year (units) ---")
    miss_by_year = (df.assign(missing_units=df["bp_units"].where(df["_tenure_missing"], 0))
                       .groupby(df[year])
                       .agg(total_units=("bp_units", "sum"),
                            missing_units=("missing_units", "sum")))
    miss_by_year["pct_missing"] = miss_by_year.missing_units / miss_by_year.total_units
    print(miss_by_year.to_string(float_format=lambda x: f"{x:,.3f}"))
    overall_missing_pct = df.loc[df["_tenure_missing"], "bp_units"].sum() / df["bp_units"].sum()
    print(f"\nOverall missing-tenure share: {overall_missing_pct:.2%}")
    if overall_missing_pct > 0.05:
        print("MATERIAL (>5%) — reporting ratio both ways (excluded vs included "
              "in denominator) below.")

    df["afford_strict"] = df[vli_dr].sum(axis=1) + df[li_dr]
    df["afford_broad"] = df["afford_strict"] + df[vli_ndr].sum(axis=1) + df[li_ndr]
    df["_county"] = df[county].astype(str).str.strip().str.upper()
    if unit_cat:
        df["is_adu"] = df[unit_cat].astype(str).str.contains("ADU", case=False)
    else:
        df["is_adu"] = False

    def report(frame, label):
        total = frame["bp_units"].sum()
        total_known_tenure = frame.loc[~frame["_tenure_missing"], "bp_units"].sum()
        owner = frame.loc[frame.is_owner, "bp_units"].sum()
        strict = frame.loc[frame.is_owner, "afford_strict"].sum()
        broad = frame.loc[frame.is_owner, "afford_broad"].sum()
        strict_no_adu = frame.loc[frame.is_owner & ~frame.is_adu,
                                  "afford_strict"].sum()
        broad_no_adu = frame.loc[frame.is_owner & ~frame.is_adu,
                                 "afford_broad"].sum()
        print(f"\n=== {label} ===")
        print(f"Total units permitted:                 {total:>12,.0f}")
        print(f"  of which known tenure:                {total_known_tenure:>12,.0f}"
              f"  ({total_known_tenure/total:6.1%})")
        print(f"Owner-tenure units:                    {owner:>12,.0f}"
              f"  ({owner/total:6.1%} of all; {owner/total_known_tenure:6.1%} of known-tenure)")
        for tag, v in [("STRICT (deed-restricted VLI+LI)", strict),
                       ("STRICT excl. ADUs", strict_no_adu),
                       ("BROAD (DR + non-DR VLI+LI)", broad),
                       ("BROAD excl. ADUs", broad_no_adu)]:
            ratio_excl = f"1 in {total_known_tenure/v:,.0f}" if v else "1 in inf"
            ratio_incl = f"1 in {total/v:,.0f}" if v else "1 in inf"
            print(f"Owner + affordable, {tag:32s}{v:>10,.0f}"
                  f"   denom=known-tenure: {ratio_excl:<14s}"
                  f"   denom=all-units: {ratio_incl}")

    report(df, "CALIFORNIA STATEWIDE (all years)")
    hg = df[df["_county"].isin(HIGH_GROWTH_COUNTIES)]
    report(hg, "HIGH-GROWTH COUNTIES (LA, Orange, SD, Bay Area 9)")

    def year_table(frame, label):
        f = frame[(frame[year] >= 2018) & (frame[year] <= 2025)].copy()
        f["oa_strict"] = f["afford_strict"].where(f.is_owner, 0)
        f["oa_strict_noadu"] = f["afford_strict"].where(f.is_owner & ~f.is_adu, 0)
        f["oa_broad"] = f["afford_broad"].where(f.is_owner, 0)
        f["oa_broad_noadu"] = f["afford_broad"].where(f.is_owner & ~f.is_adu, 0)
        g = f.groupby(f[year])[["bp_units", "oa_strict", "oa_strict_noadu",
                                 "oa_broad", "oa_broad_noadu"]].sum()
        for c in ["oa_strict", "oa_strict_noadu", "oa_broad", "oa_broad_noadu"]:
            g[c.replace("oa_", "1_in_")] = (g.bp_units / g[c].replace(0, float("nan"))).round(0)
        print(f"\n--- {label}: by year, all four variants (2018-2025 only) ---")
        print(g.to_string(float_format=lambda x: f"{x:,.3f}"))
        return g

    yr = year_table(df, "STATEWIDE")
    yr_hg = year_table(hg, "HIGH-GROWTH COUNTIES")

    # ---- Sanity check against published statewide annual permit totals ---
    print("\n--- Sanity check: statewide total permitted units by year ---")
    totals_by_year = df.groupby(df[year])["bp_units"].sum()
    print(totals_by_year.to_string(float_format=lambda x: f"{x:,.0f}"))
    print("Compare against published CA annual building-permit counts "
          "(~100-110K units/yr in recent years). Flag if any year is off "
          "by more than ~15%.")

    out = Path(__file__).parent / "apr_ownership_affordability_summary.csv"
    yr.to_csv(out)
    print(f"\nWrote {out}")

    # ---- Spot check: Encinitas La Costa 48 ----
    print("\n--- Spot check: Encinitas rows (looking for La Costa 48-style project) ---")
    print("NOTE: verified against the raw (pre-filter) CSV that the 'La Costa 48' "
          "master row (48 units, 4 deed-restricted VLI, Owner tenure, 55-yr term) "
          "matches known project facts exactly, but is ENTITLEMENT-STAGE ONLY — all "
          "its BP_* columns are zero and BP_ISSUE_DT1 is blank. Per permit-stage-only "
          "scope, it is correctly excluded from this analysis (contributes to neither "
          "numerator nor denominator) because it has not yet pulled a building permit. "
          "A separate cluster of 'La Costa 48 Design Review' single-lot rows (2025) "
          "DOES appear at BP stage below — confirm whether that's a related phase.")
    jur_col = find_col(cols, r"jur|juris", required=False)
    proj_col = find_col(cols, r"proj.*name", required=False)
    if jur_col:
        enc_all = df[df[jur_col].astype(str).str.contains("Encinitas", case=False, na=False)]
        print(f"Found {len(enc_all)} Encinitas rows total (bp_units > 0).")
        show_cols = [c for c in [jur_col, proj_col, year, tenure, unit_cat,
                                  acutely_dr, extremely_dr, vlow_dr, vlow_ndr,
                                  li_dr, li_ndr, "bp_units"] if c]
        enc = enc_all
        if proj_col:
            enc = enc_all[enc_all[proj_col].astype(str).str.contains(
                "la costa", case=False, na=False)]
            print(f"Filtered to {len(enc)} rows matching project name 'La Costa'.")
        if len(enc):
            with pd.option_context("display.max_rows", 200, "display.width", 200):
                print(enc[show_cols].to_string())
        else:
            print("No 'La Costa' project rows found in Encinitas — check project "
                  "name text or inspect enc_all manually.")
    else:
        print("Could not find a jurisdiction column to spot-check Encinitas.")


if __name__ == "__main__":
    main()
