"""Who can vote in New York City: the analysis behind a SparkyData article.

    py src/build_article.py who-can-vote

The question: a week before the 3 November 2026 election, who makes up the people allowed to vote in
the city, how has that changed since 1990, and how many adults cannot vote because they are not
citizens?

Datasets (cards in datasets/):
- ``census_history``: the 1990 and 2000 censuses and the ACS one-year files 2006-2014 (2005's has no
  group quarters and is left out): age, birthplace (born abroad or not) and citizenship of every
  resident, with borough.
- ``acs1_pums``: the ACS one-year files 2015-2024: the same, with country of birth, year of
  naturalization and year of entry, and replicate margins; for 2024 the raw state file's allocation
  flags (``acs1.raw_persons``), which mark the answers the Census Bureau filled in.
- ``cps_basic``: the CPS monthly files to September 2026: the current year, January to September 2026
  against the same months of 2025.
The census records citizenship, never legal status, and it cannot say who is entitled to naturalize.
Its adults include people in prisons and other group quarters, so "adult citizens" is wider than
the people entitled to register (the Board's rules exclude people in prison for a felony), by an
amount these files cannot measure.

Deterministic: no sampling.
"""

from __future__ import annotations

import json

import numpy as np
import pandas as pd

from _paths import REFERENCE
from article_kit import Article, lines, r, strip
from nycdata import acs1, cps, elections, geo, history, pums
from sparkylib.stats import moe_difference, moe_sum

SLUG = "who-can-vote"
SEED = 20261028
FIRST, LATEST = 1990, 2024
CPS_NOW = ((2026, 1), (2026, 9))       # January to September 2026, the current year
CPS_THEN = ((2025, 1), (2025, 9))      # the same months of 2025
BOROUGHS = ["Bronx", "Brooklyn", "Manhattan", "Queens", "Staten Island"]
AGES = [("18_29", "18 to 29", 18, 29), ("30_44", "30 to 44", 30, 44), ("45_64", "45 to 64", 45, 64), ("65_plus", "65 and over", 65, 200)]
REGIONS = [("caribbean_central", "Caribbean, Mexico and Central America"), ("asia", "Asia"), ("europe", "Europe"),
           ("south_america", "South America"), ("africa", "Africa"), ("other_world", "Canada, Oceania and elsewhere")]
MONTHS = ["January", "February", "March", "April", "May", "June", "July", "August", "September", "October",
          "November", "December"]
COURT = ("nyctapp_fossella_v_adams_2025", "nyctapp_fossella_v_adams_2025-03-20.pdf")
QUALIFY = ("nyc_boe_voter_qualifications", "nyc_boe/voter_qualifications_2026.html")
DEADLINES = ("nyc_boe_registration_deadlines_2026", "nyc_boe/registration_deadlines_2026.html")
BOE25 = ("nyc_boe_annual_report_2025", "nyc_boe/BOE_Annual_Report_2025.pdf")
B05003 = REFERENCE / "tables" / "acs_b05003_nyc_2024_1yr.json"      # source acs_b05003_nyc_2024_1yr
CITY = "16000US3651000"
#: B05003 cells for adults (18 and over), men and women: native, naturalized, not a citizen
B05003_ADULTS = {"native": ("B05003009", "B05003020"), "naturalized": ("B05003011", "B05003022"),
                 "noncitizen": ("B05003012", "B05003023")}


# ------------------------------------------------------------------ the census years


def sample(year: int):
    """One survey year's adults in the city: the frame, the estimator, and the masks every estimate uses
    (native-born citizens, naturalized citizens, non-citizens), with borough and age."""
    if year < 2015:
        s = history.load(year, ["age", "foreign_born", "noncitizen"])
        d = s.data
        nyc = np.ones(len(d), bool)
        age = d["age"].to_numpy()
        fb, nc = (d["foreign_born"] == 1).to_numpy(), (d["noncitizen"] == 1).to_numpy()
        boro = d["borough"].to_numpy()
        extra = None
    else:
        s = acs1.load(year, "person", ["AGEP", "NATIVITY", "CIT"] + (["POBP", "CITWP", "YOEP"] if year == LATEST else []))
        d = s.data
        assert d["CIT"].between(1, 5).all() and d["NATIVITY"].isin([1, 2]).all(), f"{year}: invalid citizenship or nativity"
        nyc = s.in_nyc()
        age = d["AGEP"].to_numpy()
        fb, nc = (d["NATIVITY"] == 2).to_numpy(), (d["CIT"] == 5).to_numpy()
        # every foreign-born resident is a naturalized citizen (CIT 4) or not a citizen (5)
        assert ((d["CIT"] == 4) == (fb & ~nc)).all(), f"{year}: CIT 4 should be exactly the foreign-born citizens"
        boro = s.borough().to_numpy()
        extra = d
    adult = nyc & (age >= 18)
    m = {"adult": adult, "citizen": adult & ~nc, "native": adult & ~fb, "naturalized": adult & fb & ~nc,
         "noncitizen": adult & nc}
    assert not (m["native"] & nc).any(), f"{year}: a native-born non-citizen"
    return s, m, age, boro, extra


def electorate_count(R, mask) -> dict:
    e = R.count(mask)
    return {"count": e.value, "moe": e.moe}


def electorate_share(R, base, num) -> dict:
    e = R.share(base, num)
    return {"share": 100 * e.value, "moe": 100 * e.moe}


def electorate_year(year: int) -> dict:
    s, m, age, boro, _ = sample(year)
    R = s.estimator(min_records=0)
    out = {k: electorate_count(R, v) for k, v in m.items()}
    gap = R.count(m["native"]) - R.count(m["naturalized"])          # replicate by replicate, one file
    out["native_minus_naturalized"] = {"count": gap.value, "moe": gap.moe}
    out["noncitizen_share"] = electorate_share(R, m["adult"], m["noncitizen"])
    out["naturalized_share"] = electorate_share(R, m["citizen"], m["naturalized"])
    out["records"] = int(m["adult"].sum())
    out["borough"] = {}
    for b in BOROUGHS:
        ib = boro == b
        out["borough"][b] = {"naturalized_share": electorate_share(R, m["citizen"] & ib, m["naturalized"]),
                             "noncitizen_share": electorate_share(R, m["adult"] & ib, m["noncitizen"]),
                             "citizens": electorate_count(R, m["citizen"] & ib)}
    out["age"] = {}
    for k, _, lo, hi in AGES:
        ia = (age >= lo) & (age <= hi)
        out["age"][k] = {"naturalized_share": electorate_share(R, m["citizen"] & ia, m["naturalized"]),
                         "noncitizen_share": electorate_share(R, m["adult"] & ia, m["noncitizen"])}
    return out


def latest_detail() -> dict:
    """The latest file in detail: where the naturalized were born, when they naturalized, and when the
    non-citizens arrived."""
    s, m, age, boro, d = sample(LATEST)
    R = s.estimator(min_records=0)
    out = {}
    region = geo.world_region(d["POBP"])
    out["region"] = {k: electorate_share(R, m["naturalized"], region == k) | {"n": round(R.count(m["naturalized"] & (region == k)).value, -3)}
                     for k, _ in REGIONS}
    labels = pums.code_labels("POBP")
    pobp = d["POBP"].astype(int).to_numpy()
    w = d[s.weight].to_numpy()
    tot = pd.Series(w[m["naturalized"]], index=pobp[m["naturalized"]]).groupby(level=0).sum().sort_values(ascending=False)
    top = []
    for code in tot.index[:8]:
        c = int(code)
        e = R.count(m["naturalized"] & (pobp == c))
        top.append({"code": c, "country": labels.get(f"{c:03d}", str(c)), "count": e.value, "moe": e.moe,
                    "share_of_citizens": 100 * R.share(m["citizen"], m["naturalized"] & (pobp == c)).value,
                    "records": int((m["naturalized"] & (pobp == c)).sum())})
    assert not (m["naturalized"] & (pobp < 100)).any(), "a naturalized adult with a domestic birthplace code (the recode's 'us')"
    out["countries"] = top
    c0, c1 = top[0]["code"], top[1]["code"]
    gap = R.count(m["naturalized"] & (pobp == c0)) - R.count(m["naturalized"] & (pobp == c1))
    out["top_two_contrast"] = {"diff": gap.value, "moe": gap.moe, "lo": gap.value - gap.moe, "hi": gap.value + gap.moe}
    citwp = pd.to_numeric(d["CITWP"], errors="coerce").to_numpy()
    out["naturalized_2015_on"] = electorate_share(R, m["citizen"], m["naturalized"] & (citwp >= 2015))
    out["naturalized_2015_on_count"] = electorate_count(R, m["naturalized"] & (citwp >= 2015))
    out["naturalized_2015_on_of_naturalized"] = electorate_share(R, m["naturalized"], citwp >= 2015)
    out["naturalized_before_2000_of_naturalized"] = electorate_share(R, m["naturalized"], citwp < 2000)
    yoep = pd.to_numeric(d["YOEP"], errors="coerce").to_numpy()
    out["noncitizen_arrived_2015_on"] = electorate_share(R, m["noncitizen"], yoep >= 2015)
    out["noncitizen_arrived_before_2000"] = electorate_share(R, m["noncitizen"], yoep < 2000)
    out["noncitizen_arrived_before_2010"] = electorate_share(R, m["noncitizen"], yoep < 2010)
    return out


def filled_in_answers() -> dict:
    """How many of the latest file's answers the Census Bureau filled in (the codebook's allocation flags,
    read from the raw state file: the extract does not keep them), and the shares the page prints with
    those answers left out. Point estimates, full-sample weight."""
    d = acs1.raw_persons(LATEST, ["AGEP", "CIT", "CITWP", "YOEP", "FCITP", "FCITWP", "FYOEP"])
    w = d["PWGTP"].to_numpy(float)
    adult = (d["in_nyc"] & (d["AGEP"] >= 18)).to_numpy()
    cit_code = d["CIT"].to_numpy()
    cit, nat, nc = adult & (cit_code <= 4), adult & (cit_code == 4), adult & (cit_code == 5)
    flag = {f: (d[f] == 1).to_numpy() for f in ("FCITP", "FCITWP", "FYOEP")}
    citwp, yoep = pd.to_numeric(d["CITWP"], errors="coerce").to_numpy(), pd.to_numeric(d["YOEP"], errors="coerce").to_numpy()

    def pct(base, num):
        return 100 * w[base & num].sum() / w[base].sum()
    return {"adults": float(w[adult].sum()),
            "citizenship_filled_in": pct(adult, flag["FCITP"]),
            "year_naturalized_filled_in": pct(nat, flag["FCITWP"]),
            "year_of_entry_filled_in": pct(nc, flag["FYOEP"]),
            "naturalized_share": pct(cit, cit_code == 4),
            "naturalized_share_as_answered": pct(cit & ~flag["FCITP"], cit_code == 4),
            "noncitizen_arrived_2015_on": pct(nc, yoep >= 2015),
            "noncitizen_arrived_2015_on_as_answered": pct(nc & ~flag["FYOEP"], yoep >= 2015),
            "naturalized_2015_on_of_naturalized": pct(nat, citwp >= 2015),
            "naturalized_2015_on_of_naturalized_as_answered": pct(nat & ~flag["FCITWP"], citwp >= 2015)}


# ------------------------------------------------------------------ the current year


def current_year() -> dict:
    """January to September 2026 against the same months of 2025, from the CPS: adults (18 and over) who
    are not citizens, and naturalized citizens among adult citizens (final weight, household-clustered
    margins)."""
    months = cps.months()
    now = [mm for mm in months if CPS_NOW[0] <= mm <= CPS_NOW[1]]
    then = [mm for mm in months if CPS_THEN[0] <= mm <= CPS_THEN[1]]
    assert len(now) == 9 and len(then) == 9, "January to September of both years"
    d = cps.load(["HRHHID", "HRHHID2", "GESTFIPS", "GTCO", "PRTAGE", "PRCITSHP", "PWSSWGT", "HRYEAR4", "HRMONTH"], months=then + now)
    # members of the armed forces have person records but no final weight; the survey's civilian
    # estimates leave them out, and so does this one
    d = d[d["PWSSWGT"] > 0].reset_index(drop=True)
    assert d["PRCITSHP"].between(1, 5).all(), "a citizenship code outside 1-5"
    des = cps.Design(d)
    ym = cps.yyyymm(d)
    nyc = cps.in_nyc(d) & (d["PRTAGE"] >= 18).to_numpy()
    nc = (d["PRCITSHP"] == 5).to_numpy().astype(float)
    nat = (d["PRCITSHP"] == 4).to_numpy().astype(float)
    cit = (d["PRCITSHP"] <= 4).to_numpy()
    one = np.ones(len(d))
    out = {}
    for k, mm in (("then", then), ("now", now)):
        inm = np.isin(ym, [y * 100 + m for y, m in mm])
        a, ma = cps.ratio_moe(d, y=nc, x=one, weight="PWSSWGT", mask=inm & nyc, design=des)
        b, mb = cps.ratio_moe(d, y=nat, x=one, weight="PWSSWGT", mask=inm & nyc & cit, design=des)
        out[k] = {"noncitizen_share": r(100 * a, 2), "noncitizen_moe": r(100 * ma, 2), "naturalized_share": r(100 * b, 2),
                  "naturalized_moe": r(100 * mb, 1), "records": int((inm & nyc).sum()),
                  "first": f"{MONTHS[mm[0][1] - 1]} {mm[0][0]}", "last": f"{MONTHS[mm[-1][1] - 1]} {mm[-1][0]}"}
    for k in ("noncitizen", "naturalized"):
        dv = out["now"][f"{k}_share"] - out["then"][f"{k}_share"]
        mv = moe_difference(out["now"][f"{k}_moe"], out["then"][f"{k}_moe"])
        out[f"{k}_change"] = {"pts": r(dv, 1), "lo": r(dv - mv, 1), "hi": r(dv + mv, 1)}
    return out


# ------------------------------------------------------------------ main


def main(argv: list[str] | None = None) -> int:
    """Run the analysis and write every file (no options)."""
    A = Article(SLUG, __file__, seed=SEED)
    P = A.printed
    Q = A.quotes()
    # 2005 is left out: that year's survey did not count group quarters, which the censuses and every later
    # year do (elections.CITIZEN_YEARS_LEFT_OUT; card census_history)
    years = sorted(y for y in history.years() + acs1.years() if FIRST <= y <= LATEST and y not in elections.CITIZEN_YEARS_LEFT_OUT)
    est = {y: electorate_year(y) for y in years}

    def rr(blk, nd=2):
        """Round an estimate block for the results file."""
        if isinstance(blk, dict) and ("count" in blk or "share" in blk):
            k = "count" if "count" in blk else "share"
            nd_ = -3 if k == "count" else nd
            return {k: r(blk[k], nd_), "moe": None if blk["moe"] is None or np.isnan(blk["moe"]) else r(blk["moe"], nd_)}
        if isinstance(blk, dict):
            return {kk: rr(v, nd) for kk, v in blk.items()}
        return blk

    A.results["long"] = {str(y): rr(v) for y, v in est.items()}
    # the Bureau's published adults by nativity and citizenship (B05003), reproduced from the 2024 file
    pub = json.loads(B05003.read_text(encoding="utf-8"))["data"][CITY]["B05003"]
    for k, cells in B05003_ADULTS.items():
        val = sum(pub["estimate"][c] for c in cells)
        moe_ = float(np.sqrt(sum(pub["error"][c] ** 2 for c in cells)))
        A.reproduce(f"Adults 18 and over, {k.replace('noncitizen', 'not a citizen')}, New York City, 2024 (B05003)",
                    published=val, published_moe90=r(moe_, 0), reproduced=r(est[LATEST][k]["count"], 0),
                    reproduced_moe90=r(est[LATEST][k]["moe"], 0), kind="count")
    a, b = est[FIRST], est[LATEST]
    g_cit = b["citizen"]["count"] - a["citizen"]["count"]
    g_nat = b["naturalized"]["count"] - a["naturalized"]["count"]
    g_nat_ = b["native"]["count"] - a["native"]["count"]
    A.results["growth"] = {"citizens": r(g_cit, -3), "naturalized": r(g_nat, -3), "native": r(g_nat_, -3),
                           "naturalized_share_of_growth": r(100 * g_nat / g_cit, 0),
                           "citizens_1990": r(a["citizen"]["count"], -3), "citizens_2024": r(b["citizen"]["count"], -3),
                           "naturalized_1990": r(a["naturalized"]["count"], -3), "naturalized_2024": r(b["naturalized"]["count"], -3),
                           "native_1990": r(a["native"]["count"], -3), "native_2024": r(b["native"]["count"], -3),
                           "native_2024_moe": r(b["native"]["moe"], -3), "naturalized_2024_moe": r(b["naturalized"]["moe"], -3)}
    A.results["growth"]["native_change_pct"] = r(100 * g_nat_ / a["native"]["count"], 1)
    # the net change depends on the base year: the native-born fell in the 1990s and grew after 2006
    for base in (2000, 2012):
        e0 = est[base]
        gc = b["citizen"]["count"] - e0["citizen"]["count"]
        gn = b["naturalized"]["count"] - e0["naturalized"]["count"]
        gv = b["native"]["count"] - e0["native"]["count"]
        A.results["growth"][f"since_{base}"] = {"citizens": r(gc, -3), "naturalized": r(gn, -3), "native": r(gv, -3),
                                                "naturalized_share_of_growth": r(100 * gn / gc, 0)}
    nat_counts = {y: v["native"]["count"] for y, v in est.items()}
    A.results["growth"]["native_2000"] = r(nat_counts[2000], -3)
    A.results["growth"]["native_fall_1990s"] = r(nat_counts[1990] - nat_counts[2000], -3)
    low = min((y for y in nat_counts if y >= 2000), key=nat_counts.get)
    A.results["growth"]["native_low_year"] = low
    assert low == 2006 and nat_counts[2024] > nat_counts[1990], "the page says the native-born fell to 2006 and passed 1990 by 2024"
    assert A.results["growth"]["since_2012"]["native"] > A.results["growth"]["since_2012"]["naturalized"], \
        "the page says the native-born added more than the naturalized since 2012"
    # the changes the page states between survey years carry their ranges (two files, independent samples);
    # changes from the 1990 and 2000 censuses have no margin calculated and are compared as estimates only
    def change(y0, y1, k):
        dv = est[y1][k]["count"] - est[y0][k]["count"]
        mv = moe_sum(est[y1][k]["moe"], est[y0][k]["moe"])
        return {"diff": r(dv, -3), "moe": r(mv, -3), "lo": r(dv - mv, -3), "hi": r(dv + mv, -3)}
    A.results["changes"] = {"native_2006_2024": change(2006, LATEST, "native"),
                            "native_since_2012": change(2012, LATEST, "native"),
                            "naturalized_since_2012": change(2012, LATEST, "naturalized"),
                            "native_gain_over_naturalized_since_2012": change(2012, LATEST, "native_minus_naturalized")}
    for k in A.results["changes"]:
        A.claim(f"changes.{k}", "up")
    # the methodologist's review of 24 October: the native-born lead since 2012 clears zero on the 2024 file (and just
    # on 2019's); the same contrast ending in 2021, 2022 or 2023 includes zero, so the page says so
    by_end = {str(y): change(2012, y, "native_minus_naturalized") for y in (2019, 2021, 2022, 2023, LATEST) if y in est}
    A.results["changes"]["native_gain_over_naturalized_since_2012_by_end_year"] = by_end
    clear = [y for y, v in by_end.items() if v["lo"] > 0]
    assert clear == ["2019", str(LATEST)], by_end
    # the age gradient: the oldest group's share less the youngest's, 2000 and 2024, and the group with the highest
    # share each year (in 2000 it was 45 to 64, not 65 and over)
    ag = {y: {k: est[y]["age"][k]["naturalized_share"]["share"] for k, *_ in AGES} for y in (2000, LATEST)}
    A.results["age_gradient"] = {str(y): {"oldest_minus_youngest": r(ag[y]["65_plus"] - ag[y]["18_29"], 1),
                                          "spread": r(max(ag[y].values()) - min(ag[y].values()), 1),
                                          "highest_group": max(ag[y], key=ag[y].get)} for y in ag}
    assert A.results["age_gradient"]["2000"]["highest_group"] == "45_64" and A.results["age_gradient"][str(LATEST)]["highest_group"] == "65_plus"
    assert A.results["age_gradient"][str(LATEST)]["oldest_minus_youngest"] > 2 * A.results["age_gradient"]["2000"]["oldest_minus_youngest"]
    # the counts on one population basis. Each survey file is weighted to the population estimates of its own vintage,
    # and the 2020 census revised the city's upward (about 5% for the 2018 and 2019 files); a share within one file is
    # unaffected, a count compared across years is not. elections.population_controls gives each file's count of the
    # city against the Bureau's revised estimate; multiplying by the ratio puts the counts on one series, assuming the
    # file's shares by citizenship are unaffected. A check on every count comparison the page makes, not a replacement
    # for the files' own counts (style guide: disclose a methods break and test the conclusion on one basis).
    pc = elections.population_controls().set_index("year")["ratio"]
    reb = {y: {k: est[y][k]["count"] * pc.loc[y] for k in ("citizen", "native", "naturalized", "noncitizen")} for y in est if y in pc.index}
    assert set(reb) == set(est), "every measured year has a control ratio"
    low_r = min((y for y in reb if y >= 2000), key=lambda y: reb[y]["native"])
    gc_r, gn_r = reb[LATEST]["citizen"] - reb[FIRST]["citizen"], reb[LATEST]["naturalized"] - reb[FIRST]["naturalized"]
    A.results["rebased"] = {
        "ratio": {str(y): r(pc.loc[y], 4) for y in reb},
        "native_1990_2000": r(reb[2000]["native"] - reb[1990]["native"], -3),
        "native_low_year": low_r,
        "native_2006_2024": r(reb[LATEST]["native"] - reb[2006]["native"], -3),
        "growth_1990_2024": {"citizens": r(gc_r, -3), "naturalized": r(gn_r, -3), "naturalized_share_of_growth": r(100 * gn_r / gc_r, 0)},
        "since_2012": {"native": r(reb[LATEST]["native"] - reb[2012]["native"], -3),
                       "naturalized": r(reb[LATEST]["naturalized"] - reb[2012]["naturalized"], -3)},
        "noncitizen_2024_between_1990_and_2000": bool(reb[1990]["noncitizen"] < reb[LATEST]["noncitizen"] < reb[2000]["noncitizen"]),
    }
    RB = A.results["rebased"]
    assert RB["native_1990_2000"] < 0 and low_r == 2006 and RB["native_2006_2024"] > 0, RB
    assert 85 <= RB["growth_1990_2024"]["naturalized_share_of_growth"] <= 95, RB
    assert RB["since_2012"]["native"] > RB["since_2012"]["naturalized"] and RB["noncitizen_2024_between_1990_and_2000"], RB
    assert 1.05 < pc.loc[2018] < 1.06 and 1.05 < pc.loc[2019] < 1.06, "the page says about 5%"
    nc_series = {y: v["noncitizen_share"]["share"] for y, v in est.items()}
    peak = max(nc_series, key=nc_series.get)
    A.results["noncitizen"] = {"1990": r(nc_series[1990], 2), "2000": r(nc_series[2000], 2), "2024": r(nc_series[LATEST], 2),
                               "peak_year": peak, "peak": r(nc_series[peak], 1),
                               "adults_2024": r(b["noncitizen"]["count"], -3), "adults_2024_moe": r(b["noncitizen"]["moe"], -3),
                               "adults_2000": r(est[2000]["noncitizen"]["count"], -3)}
    nat_series = {y: v["naturalized_share"]["share"] for y, v in est.items()}
    # the page: the share has held near a third since 2012 (every year from 2012 within two points of a third)
    assert all(abs(v - 100 / 3) < 2.5 for y, v in nat_series.items() if y >= 2012), {y: v for y, v in nat_series.items() if y >= 2012}
    A.results["naturalized"] = {"1990": r(nat_series[1990], 2), "2000": r(nat_series[2000], 2), "2024": r(nat_series[LATEST], 2),
                                "2024_moe": r(est[LATEST]["naturalized_share"]["moe"], 1),
                                "max_year": max(nat_series, key=nat_series.get), "max": r(max(nat_series.values()), 1),
                                "since_2012_range": [r(min(v for y, v in nat_series.items() if y >= 2012), 1),
                                                     r(max(v for y, v in nat_series.items() if y >= 2012), 1)],
                                "since_2012_min": r(min(v for y, v in nat_series.items() if y >= 2012), 2),
                                "since_2012_max": r(max(v for y, v in nat_series.items() if y >= 2012), 2)}
    A.results["detail"] = rr(latest_detail())
    F = filled_in_answers()
    # the raw file's city adults are the extract's
    assert abs(F["adults"] - est[LATEST]["adult"]["count"]) < 1, (F["adults"], est[LATEST]["adult"]["count"])
    A.results["filled_in"] = {k: r(v, 2) for k, v in F.items()}
    dict_line = pums.dictionary_entries(["FCITP"])["FCITP"][0]
    assert dict_line.endswith(",Citizenship allocation flag"), dict_line
    A.results["filled_in"]["flag_label"] = "Citizenship allocation flag"
    # the page: leaving the filled-in answers out moves the three shares checked by at most the largest move printed
    # (the arrivals since 2015, about 2.2 points) and the naturalized share by less than a point
    moves = {a_: abs(F[a_] - F[b_]) for a_, b_ in (("naturalized_share", "naturalized_share_as_answered"),
                                                   ("noncitizen_arrived_2015_on", "noncitizen_arrived_2015_on_as_answered"),
                                                   ("naturalized_2015_on_of_naturalized", "naturalized_2015_on_of_naturalized_as_answered"))}
    A.results["filled_in_largest_move"] = {"share": max(moves, key=moves.get), "points": r(max(moves.values()), 2)}
    assert max(moves, key=moves.get) == "noncitizen_arrived_2015_on" and max(moves.values()) < 2.5, moves
    # and: "the other shares checked move by less than a point"
    assert abs(F["naturalized_share"] - F["naturalized_share_as_answered"]) < 1
    assert abs(F["naturalized_2015_on_of_naturalized"] - F["naturalized_2015_on_of_naturalized_as_answered"]) < 1
    A.results["cps_2026"] = current_year()
    C = A.results["cps_2026"]
    A.latest("2026-09", "CPS basic monthly files, January to September 2026 against the same months of 2025 (citizenship of adults)")
    A.long_view(FIRST, LATEST, "Census 1990 and 2000 PUMS, ACS one-year PUMS 2006-2024 (2005 left out: no group "
                "quarters): adults by citizenship and birthplace")

    # ---- the borough and age splits, 2024 (and 2000 beside it)
    BB = A.results["long"][str(LATEST)]["borough"]
    nb = {k: v["naturalized_share"]["share"] for k, v in BB.items()}
    A.results["borough_order"] = sorted(nb, key=nb.get, reverse=True)
    AG = A.results["long"][str(LATEST)]["age"]

    # ---- the documents
    hist = {}
    hist["court_holding"] = Q.quote(*COURT, "Article II, Section 1 of our Constitution limits voting to citizens")
    hist["qualifications"] = Q.quote(*QUALIFY, "Be a citizen of the United States")
    hist["deadline_2026"] = Q.quote(*DEADLINES, "Last day to register to vote in person is October 24, 2026")
    hist["law_definition"] = Q.quote(*COURT, "is either a lawful permanent resident or authorized to work in the United States")
    hist["enjoined"] = Q.quote(*COURT, "permanently enjoined the City from implementing the law")
    hist["decided"] = Q.quote(*COURT, "Decided March 20, 2025")
    hist["effective"] = Q.quote(*COURT, "became effective on January 9, 2022")
    hist["residence"] = Q.quote(*COURT, "will have been such a resident for 30 consecutive days or longer")
    hist["law_other_qualifications"] = Q.quote(*COURT, "meets all qualifications for registering or preregistering to vote "
                                                       "under the election law, except for possessing United States citizenship")
    hist["not_a_citizen_on_election_day"] = Q.quote(*COURT, "a person who is not a United States citizen on the date of the "
                                                            "election on which he or she is voting")
    hist["affirmed"] = Q.quote(*COURT, "limits voting to citizens, and affirm")
    hist["law_registered"] = Q.quote(*COURT, "has registered or preregistered to vote with the board of elections in the city of New York")
    hist["rules"] = Q.quote(*QUALIFY, "Be a New York City resident for at least 30 days. Be 18 years old on or before Election Day. "
                                      "Not be in prison for a felony conviction. Not be adjudged mentally incompetent by a court.")
    hist["territories"] = Q.quote(*QUALIFY, "those persons born in Puerto Rico, Guam and the U.S. Virgin Islands")
    hist["offices_2026"] = Q.quote(*BOE25, "voters will elect members of the New York State Legislature, the Governor")
    hist["acs_2005_gq"] = Q.quote("census_acs_comparing_2006", "census_acs_comparing_2006.html",
                                  "is included in the 2006 ACS and not included in the 2005 ACS")
    hist["qualifications_territories"] = Q.quote(*QUALIFY, "includes those persons born in Puerto Rico, Guam and the U.S. Virgin Islands")
    # the Bureau's update of 22 October 2026 (read 26 October): the 2025 one-year estimates scheduled for 29 October, the
    # day after this page is released; and its microdata folder as listed on 26 October 2026, 2024 still the last year
    hist["acs_2025_scheduled"] = Q.quote("census_acs_2025_release_scheduled", "census_acs_updates_2026-10-26.html",
                                         "scheduled for Thursday, October 29, 2026")
    hist["pums_folder_latest"] = Q.quote("census_acs_pums_directory_2026_10_26", "census_acs_pums_directory_2026-10-26.html",
                                         "2023/ 2025-01-23 07:24 - 2024/ 2026-01-29 10:06 -")
    assert "2025/" not in (REFERENCE / "reports" / "census_acs_pums_directory_2026-10-26.html").read_text(encoding="utf-8", errors="replace")
    A.results["quotes"] = {k: {"quote": v["quote"], "page": v.get("page"), "source": v.get("source")} for k, v in hist.items()}
    # the codebook's own words for the five citizenship groups (data/reference/PUMS_Data_Dictionary_2020-2024.csv)
    labels = pums.code_labels("CIT")
    A.results["cit_labels"] = labels
    assert labels["2"] == "Born in Puerto Rico, Guam, the U.S. Virgin Islands, or the Northern Marianas", labels

    # ---- printed figures
    G, N, NT, D = A.results["growth"], A.results["noncitizen"], A.results["naturalized"], A.results["detail"]
    P("growth.citizens_2024", f"{G['citizens_2024'] / 1e6:.1f} million")
    P("growth.citizens", f"{G['citizens'] / 1e6:.2f} million")
    P("growth.naturalized_1990", f"{G['naturalized_1990']:,.0f}")
    P("growth.naturalized_2024", f"{G['naturalized_2024']:,.0f}")
    P("growth.native_1990", f"{G['native_1990'] / 1e6:.2f} million")
    P("growth.native_2024", f"{G['native_2024'] / 1e6:.2f} million")
    A.results["growth"]["native_2006"] = r(nat_counts[2006], -3)
    P("growth.native_2006", f"{A.results['growth']['native_2006'] / 1e6:.2f} million")
    A.results["noncitizen"]["adults_1990"] = r(est[1990]["noncitizen"]["count"], -3)
    P("noncitizen.adults_1990", f"{A.results['noncitizen']['adults_1990'] / 1e6:.2f} million")
    P("noncitizen.adults_2000", f"{N['adults_2000'] / 1e6:.2f} million")
    P("naturalized.since_2012_min", f"{NT['since_2012_min']:.0f}%")
    P("naturalized.since_2012_max", f"{NT['since_2012_max']:.0f}%")
    for k in ("18_29", "65_plus"):
        P(f"long.2000.age.{k}.naturalized_share.share", f"{A.results['long']['2000']['age'][k]['naturalized_share']['share']:.0f}%")
    P("detail.naturalized_2015_on_count.count", f"{D['naturalized_2015_on_count']['count']:,.0f}")
    P("growth.since_2012.native", f"{G['since_2012']['native']:,.0f}")
    P("growth.since_2012.naturalized", f"{G['since_2012']['naturalized']:,.0f}")
    P("naturalized.1990", f"{NT['1990']:.0f}%")
    P("naturalized.2000", f"{NT['2000']:.0f}%")
    P("naturalized.2024", f"{NT['2024']:.0f}%")
    P("noncitizen.1990", f"{N['1990']:.0f}%")
    P("noncitizen.2000", f"{N['2000']:.0f}%")
    P("noncitizen.2024", f"{N['2024']:.0f}%")
    P("noncitizen.adults_2024", f"{N['adults_2024'] / 1e6:.1f} million")
    assert est[1990]["noncitizen"]["count"] < est[LATEST]["noncitizen"]["count"] < est[2000]["noncitizen"]["count"], \
        "the page says the count is above 1990's and below 2000's"
    P("cps_2026.now.noncitizen_share", f"{C['now']['noncitizen_share']:.0f}%")
    P("cps_2026.then.noncitizen_share", f"{C['then']['noncitizen_share']:.0f}%")
    P("cps_2026.now.naturalized_share", f"{C['now']['naturalized_share']:.0f}%")
    P("cps_2026.then.naturalized_share", f"{C['then']['naturalized_share']:.0f}%")
    for bname in ("Queens", "Manhattan", "Staten Island"):
        P(f"long.{LATEST}.borough.{bname}.naturalized_share.share", f"{BB[bname]['naturalized_share']['share']:.0f}%")
    P("long.2000.borough.Staten Island.naturalized_share.share",
      f"{A.results['long']['2000']['borough']['Staten Island']['naturalized_share']['share']:.0f}%")
    for k in ("18_29", "65_plus"):
        P(f"long.{LATEST}.age.{k}.naturalized_share.share", f"{AG[k]['naturalized_share']['share']:.0f}%")
    P("detail.naturalized_2015_on.share", f"{D['naturalized_2015_on']['share']:.0f}%")
    P("detail.noncitizen_arrived_2015_on.share", f"{D['noncitizen_arrived_2015_on']['share']:.0f}%")
    # the two largest countries of birth are within each other's margins: printed as "about 200,000 each"
    c0, c1 = D["countries"][0], D["countries"][1]
    TT = D["top_two_contrast"]
    tie = TT["lo"] < 0 < TT["hi"]
    assert tie and {c0["country"], c1["country"]} == {"Dominican Republic", "China"}, (c0, c1, TT)
    D["top_two"] = {"countries": sorted([c0["country"], c1["country"]]), "about": r((c0["count"] + c1["count"]) / 2, -5)}
    P("detail.top_two.about", f"{D['top_two']['about']:,.0f}")
    # the Caribbean (with Mexico and Central America) and Asia: about a third each, not ranked
    rc, ra = D["region"]["caribbean_central"], D["region"]["asia"]
    assert abs(rc["share"] - ra["share"]) < 2 and 28 < min(rc["share"], ra["share"]), (rc, ra)
    # the page's growth claim is the net change from 1990: printed with its base years beside it
    P("growth.naturalized_share_of_growth", f"{G['naturalized_share_of_growth']:.0f}%")
    CG = A.results["changes"]
    P("changes.native_2006_2024.diff", f"{CG['native_2006_2024']['diff']:,.0f}")
    P("changes.native_2006_2024.moe", f"{CG['native_2006_2024']['moe']:,.0f}")
    P("changes.native_gain_over_naturalized_since_2012.diff", f"{CG['native_gain_over_naturalized_since_2012']['diff']:,.0f}")
    P("changes.native_gain_over_naturalized_since_2012.moe", f"{CG['native_gain_over_naturalized_since_2012']['moe']:,.0f}")
    FI = A.results["filled_in"]
    P("filled_in.citizenship_filled_in", f"{FI['citizenship_filled_in']:.0f}%")
    P("filled_in.year_naturalized_filled_in", f"{FI['year_naturalized_filled_in']:.0f}%")
    P("filled_in.year_of_entry_filled_in", f"{FI['year_of_entry_filled_in']:.0f}%")
    P("filled_in.noncitizen_arrived_2015_on_as_answered", f"{FI['noncitizen_arrived_2015_on_as_answered']:.0f}%")
    P("naturalized.2024_moe", f"{NT['2024_moe']:.1f}")
    P("rebased.native_2006_2024", f"{r(RB['native_2006_2024'], -4):,.0f}")
    P("rebased.growth_1990_2024.citizens", f"{RB['growth_1990_2024']['citizens'] / 1e6:.2f} million")
    P("rebased.growth_1990_2024.naturalized_share_of_growth", f"{RB['growth_1990_2024']['naturalized_share_of_growth']:.0f}%")
    P("rebased.ratio.2018", f"{100 * (RB['ratio']['2018'] - 1):.0f}%")
    for y in ("2000", str(LATEST)):
        P(f"age_gradient.{y}.spread", f"{A.results['age_gradient'][y]['spread']:.0f}")
    P("filled_in_largest_move.points", f"{A.results['filled_in_largest_move']['points']:.1f}")
    BE = A.results["changes"]["native_gain_over_naturalized_since_2012_by_end_year"]
    P("changes.native_gain_over_naturalized_since_2012_by_end_year.2023.diff", f"{abs(BE['2023']['diff']):,.0f}")
    P("changes.native_gain_over_naturalized_since_2012_by_end_year.2023.moe", f"{BE['2023']['moe']:,.0f}")
    for k in ("noncitizen", "naturalized"):          # a minus sign, not a hyphen, on the page
        P(f"cps_2026.{k}_change.lo", f"{C[f'{k}_change']['lo']:.1f}".replace("-", "−"))
        P(f"cps_2026.{k}_change.hi", f"{C[f'{k}_change']['hi']:.1f}".replace("-", "−"))

    # ---- charts
    xs = [[str(y), str(y) + (" (census)" if y < 2005 else "")] for y in years]
    ticks = [str(y) for y in years if y in (1990, 2000, 2010, 2020, LATEST)]
    L = A.results["long"]
    A.chart("electorate", lines(
        title="Naturalized citizens rose from about a sixth of the city's adult citizens to a third; the native-born fell, then grew again",
        subtitle=f"Adult citizens (18 and over) in New York City by birthplace, {FIRST} to {LATEST}, millions",
        x=xs, series=[["native", "Born in the United States, its territories or abroad to American parents", "Native-born"],
                      ["naturalized", "Born abroad, naturalized", "Naturalized"]],
        cells={"": {k: {**{str(y): r(L[str(y)][k]["count"] / 1e6, 2) for y in years},
                        **{f"{y}_moe": r(L[str(y)][k]["moe"] / 1e6, 2) for y in years if L[str(y)][k]["moe"] is not None}}
                    for k in ("native", "naturalized")}},
        format="m", axis="Adult citizens (millions)", labels="legend", x_ticks=ticks, x_tick_labels={t: t for t in ticks},
        note="The 1990 and 2000 censuses (5% samples, no margin calculated), then the American Community Survey one-year files "
             "from 2006 (2001-2005 have no comparable file: 2005's survey left out people in group quarters, such as prisons "
             "and dormitories). Native-born includes people born in Puerto Rico and the other territories and abroad to "
             "American parents. Each survey year is weighted to that year's population estimates, which the 2020 census "
             "revised upward for the city, so the counts before 2020 sit on a lower base than those after (the questions "
             "below). Ticks are 90% margins where the file carries them; the 2020 file, an experimental release, has none.",
        provenance="U.S. Census Bureau, decennial census and American Community Survey microdata"))
    A.chart("noncitizen_share", lines(
        title="One adult New Yorker in six is not a citizen, fewer than in 2000",
        subtitle=f"Adults (18 and over) who are not U.S. citizens, share of the city's adults, {FIRST} to {LATEST}",
        x=xs, series=[["noncitizen", "Not a citizen", "Not a citizen"]],
        cells={"": {"noncitizen": {**{str(y): r(L[str(y)]["noncitizen_share"]["share"], 1) for y in years},
                                   **{f"{y}_moe": L[str(y)]["noncitizen_share"]["moe"] for y in years
                                      if L[str(y)]["noncitizen_share"]["moe"] is not None}}}},
        format="pct", axis="Share of adults", x_ticks=ticks, x_tick_labels={t: t for t in ticks},
        note="A share within one file, which the scalar rebasing in the questions leaves unchanged. Ticks are 90% margins "
             "from the replicate weights, where the file carries them; the 2020 file, an experimental release, has none.",
        provenance="U.S. Census Bureau, decennial census and American Community Survey microdata"))
    A.chart("boroughs", strip(
        title="In Queens, more than four in ten adult citizens are naturalized",
        subtitle=f"Naturalized citizens, share of each borough's adult citizens, 2000 and {LATEST}",
        x=[["2000", "2000 (census)", 1], [str(LATEST), str(LATEST), 0]],
        series=[[bname, bname, bname] for bname in BOROUGHS],
        cells={"": {bname: {"2000": r(L["2000"]["borough"][bname]["naturalized_share"]["share"], 1),
                            str(LATEST): BB[bname]["naturalized_share"]["share"],
                            f"{LATEST}_moe": BB[bname]["naturalized_share"]["moe"]} for bname in BOROUGHS}},
        format="pct", axis="Share of adult citizens", value_x=str(LATEST), value_label=str(LATEST), moe_label="90% margin",
        note="The 2000 census (no margin calculated) and the 2024 American Community Survey one-year file, by the borough of residence.",
        provenance="U.S. Census Bureau, decennial census and American Community Survey microdata"))
    A.chart("by_age", strip(
        title="The youngest adult citizens are the least likely to be naturalized",
        subtitle=f"Naturalized citizens, share of adult citizens in each age group, New York City, 2000 and {LATEST}",
        x=[["2000", "2000 (census)", 1], [str(LATEST), str(LATEST), 0]],
        series=[[k, lab, lab] for k, lab, *_ in AGES],
        cells={"": {k: {"2000": r(L["2000"]["age"][k]["naturalized_share"]["share"], 1),
                        str(LATEST): AG[k]["naturalized_share"]["share"],
                        f"{LATEST}_moe": AG[k]["naturalized_share"]["moe"]} for k, *_ in AGES}},
        format="pct", axis="Share of adult citizens", value_x=str(LATEST), value_label=str(LATEST), moe_label="90% margin",
        order="series",
        note="The 2000 census and the 2024 American Community Survey one-year file.",
        provenance="U.S. Census Bureau, decennial census and American Community Survey microdata"))
    A.chart("naturalized_share", lines(
        title="Naturalized citizens have been about a third of the city's adult citizens since 2012, up from a sixth in 1990",
        subtitle=f"Naturalized citizens, share of adult citizens (18 and over), New York City, {FIRST} to {LATEST}",
        x=xs, series=[["naturalized", "Naturalized, share of adult citizens", "Naturalized"]],
        cells={"": {"naturalized": {**{str(y): r(L[str(y)]["naturalized_share"]["share"], 1) for y in years},
                                    **{f"{y}_moe": L[str(y)]["naturalized_share"]["moe"] for y in years
                                       if L[str(y)]["naturalized_share"]["moe"] is not None}}}},
        format="pct", axis="Share of adult citizens", x_ticks=ticks, x_tick_labels={t: t for t in ticks},
        note="The 1990 and 2000 censuses (5% samples, no margin calculated), then the American Community Survey one-year files "
             "from 2006 (2001-2005 have no comparable file: 2005's survey left out people in group quarters, such as prisons "
             "and dormitories; 2020's was an experimental release, with no margin). A share within one file, which the "
             "scalar rebasing in the questions leaves unchanged. Ticks are 90% margins from the replicate weights.",
        provenance="U.S. Census Bureau, decennial census and American Community Survey microdata"))

    rows = []
    for y in years:
        e = L[str(y)]
        for k in ("adult", "citizen", "native", "naturalized", "noncitizen"):
            rows.append({"table": "adults_by_citizenship", "year": y, "group": k, "area": "New York City",
                         "value": e[k]["count"], "unit": "people", "moe90": e[k]["moe"], "moe90_unit": "people"})
        for bname in BOROUGHS:
            rows.append({"table": "naturalized_share_of_adult_citizens", "year": y, "group": "naturalized", "area": bname,
                         "value": e["borough"][bname]["naturalized_share"]["share"], "unit": "percent of adult citizens",
                         "moe90": e["borough"][bname]["naturalized_share"]["moe"], "moe90_unit": "percentage points"})
    for c in D["countries"]:
        rows.append({"table": f"naturalized_adult_citizens_by_country_{LATEST}", "year": LATEST, "group": c["country"],
                     "area": "New York City", "value": c["count"], "unit": "people", "moe90": c["moe"], "moe90_unit": "people"})
    A.write(data=pd.DataFrame(rows))
    print(f"{SLUG}: wrote {A.out}")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())
