"""exercise_3b.py

Exercise III-B. Blind-spot mapping and the reliability screen.

Computes the coefficient of variation for every median household income
estimate at three geographies, counts the units failing a thirty percent
coefficient of variation, a reliability cutoff commonly applied to American
Community Survey estimates, finds what the ranking of the poorest
block groups rested on, and tests whether small and large tracts are
represented by evidence of the same quality.

    python exercise_3b.py
"""
from __future__ import annotations

import sys
from pathlib import Path

import numpy as np
import pandas as pd

sys.path.insert(0, str(Path(__file__).resolve().parent))
import common as C                                     # noqa: E402

CV_LIMIT = 30.0
MOE_TO_SD = 1.645        # the Census Bureau publishes a 90 percent margin


def cv_of(d):
    """Coefficient of variation in percent.

    The published margin is a 90 percent interval, so dividing by 1.645
    recovers the standard error. Using 1.96 here would treat it as a 95
    percent interval and understate every coefficient by about sixteen
    percent, with no error message and a map that still looks reasonable.
    """
    return 100 * (d.medinc_moe / MOE_TO_SD) / d.medinc


def main():
    C.ensure_dirs()
    m = C.load_acs()
    C.banner("exercise III-B: blind-spot mapping and the reliability screen",
             "ACS 2023 5-year, Mississippi, estimates and margins of error",
             "acs2023_mississippi_asi.csv")
    out = {"cv_limit_pct": CV_LIMIT}

    print("\n[1] THE SCREEN AT THREE GEOGRAPHIES")
    print(f"    {'geography':<12}{'present':>9}{'suppressed':>12}"
          f"{'CV over 30%':>13}{'median CV':>11}")
    rows = []
    for lev in ("block group", "tract", "county"):
        d = m[m.level == lev]
        dd = d.dropna(subset=["medinc", "medinc_moe"])
        cv = cv_of(dd)
        r = {"geography": lev, "present": int(len(dd)),
             "suppressed": int(d.medinc.isna().sum()),
             "pct_over_limit": round(100 * float((cv > CV_LIMIT).mean()), 1),
             "median_cv": round(float(cv.median()), 1)}
        rows.append(r)
        print(f"    {lev:<12}{r['present']:>9d}{r['suppressed']:>12d}"
              f"{r['pct_over_limit']:>12.1f}%{r['median_cv']:>10.1f}%")
    out["screen"] = rows

    print("\n[2] WHAT THE RANKING RESTED ON")
    bg = m[m.level == "block group"].dropna(subset=["medinc", "medinc_moe"])
    bg = bg.assign(cv=cv_of(bg))
    poorest = bg.nsmallest(10, "medinc")
    fail = int((poorest.cv > CV_LIMIT).sum())
    print(f"    the ten lowest income block groups in the state")
    print(f"    failing the reliability screen: {fail} of 10")
    print(f"    their median coefficient of variation: "
          f"{poorest.cv.median():.1f}%")
    out["poorest10"] = {"failing": fail,
                        "median_cv": round(float(poorest.cv.median()), 1)}

    print("\n[3] INFORMATION PARITY")
    tr = m[m.level == "tract"].dropna(subset=["medinc", "medinc_moe", "pop"])
    tr = tr[tr["pop"] > 0].assign(cv=cv_of(tr[tr["pop"] > 0]))
    tr["pop_quintile"] = pd.qcut(tr["pop"], 5, labels=[1, 2, 3, 4, 5])
    print(f"    {'population quintile':<22}{'tracts':>8}{'median CV':>11}"
          f"{'CV over 30%':>13}")
    par = []
    for q, g in tr.groupby("pop_quintile", observed=True):
        r = {"quintile": int(q), "tracts": int(len(g)),
             "median_cv": round(float(g.cv.median()), 1),
             "pct_over_limit": round(100 * float((g.cv > CV_LIMIT).mean()), 1)}
        par.append(r)
        tag = "  smallest tracts" if q == 1 else ("  largest tracts"
                                                  if q == 5 else "")
        print(f"    {int(q):<22}{r['tracts']:>8d}{r['median_cv']:>10.1f}%"
              f"{r['pct_over_limit']:>12.1f}%{tag}")
    out["parity"] = par

    # The three-panel map is drawn by the shared figure module so that the
    # exercise and the workbook show the same image.
    import make_figures as MF
    MF.fig_reliability()

    tbl = bg[["code", "medinc", "medinc_moe", "cv"]].copy()
    tbl["passes_screen"] = tbl.cv <= CV_LIMIT
    C.write_handoff("exercise_3b", tbl, {
        "reliability_limit_pct": CV_LIMIT,
        "moe_to_standard_error": "divide by 1.645, a 90 percent margin",
        "suppressed_units": "shown in gray, never imputed",
        "classification": "quantile, five classes",
        "color_ramp": "YlGnBu for the estimate, OrRd for reliability, "
                      "unjustified",
    }, [
        "The screen is a threshold on a continuous quantity. A unit at 29 "
        "percent passes and one at 31 fails. Should the map show that "
        "boundary as sharply as it does?",
        "Suppressed units are gray on every panel. Does gray read as missing "
        "or as low to someone who has not read the legend?",
        "Which of the three maps would you publish, and what goes in the "
        "caption of the one you chose?",
        "Small tracts carry weaker evidence. Who lives in them, and does the "
        "screened map make those communities disappear?",
    ])
    C.save_checkvalues("exercise_3b", out)


if __name__ == "__main__":
    main()