"""Raising a child in New York: where the city's families live, who minds the children,
what care costs and who provides it, and which families stay — every number and figure
in the article, with the working.

    py src/article_parents.py

This file is the article's notebook. The article keeps to what was found; the how and
the why live here, next to the code that computes them.

======================================================================================
1. The question, and the data that answer it
======================================================================================

Half a million children under six live in the five boroughs. The published census
table that comes closest to describing their parents (B23008, own children under 18
by living arrangement and the employment status of their parents) says how many live
with two working parents, with one parent at home, or with a lone parent. It does not
say how that changes with the child's age, where those families live, what their
households earn, whether the three-year-olds are in a preschool, who the nannies are,
or how many of the families leave. The microdata behind the table does.

Six sources, and the article says which one each number comes from.

* **The published table.** B23008 for New York city from the one-year survey: 2019
  from the Bureau's sequence-format summary file (sequence 101), 2021-2023 from the
  table-based summary files, 2024 from the Census Reporter mirror. The article
  reproduces every row before going anywhere the table does not.

* **The census microdata.** The 2015-2024 one-year PUMS files for the metro study
  area (data/interim/acs1/metro_person_<year>.csv.gz and metro_housing_<year>.csv.gz,
  with their 80 replicate weights), cut to residents of the five boroughs. Children
  carry the employment status of the parents they live with (the ESP item, whose
  universe is the published table's), school enrolment from age three, and the
  household's income, tenure and rent. Adults carry their own employment, hours,
  earnings, occupation and industry, which is how the nannies and the day-care
  workers are found, and whether a woman gave birth in the past year.

* **Leavers**, from the national one-year files (data/interim/migration/
  leavers_<year>.csv.gz, built by extract_migration.py): every person who lived in
  the five boroughs a year earlier and answered the census from elsewhere in the
  United States, joined to the household they live in now.

* **The State's child-care rate schedule** (OCFS market rates effective 1 October
  2024, county group 4, fetched from the Day Care Council of New York's copy because
  ocfs.ny.gov blocks scripted clients) and the City Comptroller's January 2025
  spotlight, which multiplies the same weekly rates by 52. Published parameters, used
  to compute what full-time care costs against the incomes the census measures.

* **The Census Bureau's population estimates** by age (Vintage 2020 for July 2010 to
  July 2020, Vintage 2025 for July 2020 to July 2025): the under-five population of
  each borough, an administrative series that runs a year past the survey.

* **Two City series.** Live births 2010-2023 from the Health Department's Summary of
  Vital Statistics 2023 (transcribed), and public-school enrolment in 3-K, pre-K and
  kindergarten for 2021-22 to 2025-26 from the Department of Education's demographic
  snapshot (read from the workbook). The State Comptroller's January 2026 report on
  child care for children under five supplies the voucher and seat counts quoted in
  the FAQ (transcribed).

======================================================================================
2. Universe and definitions
======================================================================================

* **Children**: residents of the five boroughs under 18. *Own children* are those
  living with at least one parent as the child of the householder or in a subfamily,
  the universe of table B23008 and of the parents'-employment item. Children living
  with grandparents or other relatives and children in group quarters are counted
  and shown as "not living with a parent".
* **Parents' work arrangement** (from the child's record): two parents both in the
  labor force; two parents with one not in the labor force (the "parent at home");
  two parents with neither in the labor force; one parent in the labor force; one
  parent not in the labor force. "All resident parents at work" groups the first and
  fourth; "a parent at home" groups the second, third and fifth.
* **Parents** (adult measures): the householder, spouse or partner in a household
  where a child of the householder under 18 lives; the *youngest child's age* is the
  youngest such child. Employed, in the labor force and working 35 hours or more are
  the person's own items. A *new mother* is a woman aged 15 to 50 who reports giving
  birth in the past twelve months.
* **Preschool**: among three- and four-year-olds, enrolled in school (public or
  private) in the past three months; at those ages that is nursery school, preschool
  or kindergarten.
* **The child-care workforce**: employed people whose occupation is childcare worker,
  split by industry into private households (the census's nannies and babysitters
  employed directly by a family), child day care services (centres and licensed
  homes) and everywhere else (schools, religious organisations, recreation); plus
  preschool and kindergarten teachers working in child day care services. Hourly
  earnings are annual earnings divided by usual hours times weeks worked, for people
  who report all three.
* **Household income**: the household's total income in 2024 dollars (the Bureau's
  within-year factor, then the New York-metro CPI). *Rent burden*: gross rent 30% or
  more of household income, renters with a computed ratio. *Poverty*: the child's
  family income below the federal poverty threshold, as the census computes it.
* **Leavers and out-migration rate**: as in the migration article, a leaver lived in
  the five boroughs a year ago and now lives elsewhere in the United States; the rate
  divides leavers by everyone who lived in the city a year ago and is still in a US
  census file. A child's household income, tenure and destination are those of the
  household they live in now, after the move.
* **Cost of care (computed)**: the State's weekly maximum rate times 52 weeks, as the
  Comptroller states annual costs: $26,000 for an infant in a centre, $20,800 in a
  family day care home; compared with the census-measured *family* income (the income
  of the householder and everyone related to them, the Bureau's FINCP, which is the
  income the 7% benchmark is defined on) of households whose youngest child is under
  three, at 7%, 10%, 20% and 30% of income. The household-income version (everyone in
  the unit, related or not) is kept as a sensitivity. The 7% figure is the benchmark
  the City Comptroller uses, drawn from the federal child-care subsidy rule of March
  2024; the federal rule making it mandatory was rescinded effective 13 July 2026.
* **Other household income** (for the parent-at-home gradient): the household's income
  less the parent's own total personal income, so that a mother's employment is not
  compared with an income measure that contains her own earnings.

======================================================================================
3. Margins
======================================================================================

PUMS ships 80 replicate weights per person and per housing unit. For any statistic θ
computed with the full weight, the same statistic is recomputed with each replicate
weight r, and Var(θ) = (4 / 80) · Σ_r (θ_r − θ)². Margins are 90% (1.645 SE), the
Bureau's convention. Medians are weighted medians without a margin. Rates are shares
computed replicate by replicate. The population estimates, the vital statistics and
the school enrolment carry no sampling margin.

Outputs (all regenerable; nothing is edited by hand):
    output/articles/raising_a_child_in_new_york/results.json      every number quoted in the prose
    output/articles/raising_a_child_in_new_york/pums_direct.csv   census-side estimates with margins
    output/articles/raising_a_child_in_new_york/reproduction.csv  the published-table reproduction
    output/articles/raising_a_child_in_new_york/charts.json       cubes behind the interactive charts
    output/articles/raising_a_child_in_new_york/fig*.svg          static twins of the same figures
    output/articles/raising_a_child_in_new_york/data.csv          companion dataset (re-randomised ids)
    site/public/articles/raising-a-child-in-new-york/             the same, copied for the website
"""

from __future__ import annotations

import argparse
import json
import shutil
from pathlib import Path

import numpy as np
import pandas as pd

from _paths import P, INTERIM, OUTPUT, RAW, REFERENCE
from article_commute import weighted_median
from article_migration import (Census, BOROUGHS, INCOME_BANDS, MIN_RECORDS, NYC_MIGPUMA, STATE_NAMES, SUBURB_REGIONS, SF_2019, Z90, MIG_DIR,
                               _geo, vintage, cpi_deflators, race_label, income_band, state_group, _grouped, _mpl, _save, INK, GRID, BLUE, ORANGE, GREEN, TICK, FIG_W)
import article_parents_charts as CH

SLUG = "raising-a-child-in-new-york"
OUT = OUTPUT / "articles" / "raising_a_child_in_new_york"
SITE_OUT = P.root / "site" / "public" / "articles" / SLUG
YEARS = list(range(2015, 2025))
FOCUS = (2015, 2019, 2024)
TABLE_2024 = REFERENCE / "tables" / "acs_b23008_nyc_us_2024_1yr.json"
DOE_SNAPSHOT = RAW / "doe" / "demographic_snapshot_2021_22_to_2025_26.xlsx"
NYC_COUNTIES = {"005": "Bronx", "047": "Brooklyn", "061": "Manhattan", "081": "Queens", "085": "Staten Island"}

AGE_BANDS = [("a0", "Under 1", 0, 1), ("a1", "Age 1", 1, 2), ("a2", "Age 2", 2, 3), ("a3", "Age 3", 3, 4), ("a4", "Age 4", 4, 5), ("a5", "Age 5", 5, 6),
             ("a6_12", "6 to 12", 6, 13), ("a13_17", "13 to 17", 13, 18)]
YOUNGEST_BANDS = AGE_BANDS
LEAVE_BANDS = [("c1_5", "Children 1 to 5", 1, 6), ("c6_12", "Children 6 to 12", 6, 13), ("c13_17", "Children 13 to 17", 13, 18), ("adults_25_44", "Adults 25 to 44", 25, 45),
               ("adults_45_plus", "Adults 45 and over", 45, 200)]
ESP_CODES = {1: "both_lf", 2: "father_lf", 3: "mother_lf", 4: "neither_lf", 5: "lone_father_lf", 6: "lone_father_nilf", 7: "lone_mother_lf", 8: "lone_mother_nilf"}
ARRANGEMENTS = [("both_lf", "Two parents, both in the labor force"), ("one_home", "Two parents, one not in the labor force"), ("neither_lf", "Two parents, neither in the labor force"),
                ("lone_lf", "One parent, in the labor force"), ("lone_home", "One parent, not in the labor force"), ("no_parent", "Not living with a parent")]
ARR_GROUP = {"both_lf": "both_lf", "father_lf": "one_home", "mother_lf": "one_home", "neither_lf": "neither_lf", "lone_father_lf": "lone_lf", "lone_mother_lf": "lone_lf",
             "lone_father_nilf": "lone_home", "lone_mother_nilf": "lone_home"}
ALL_WORK = {"both_lf", "lone_father_lf", "lone_mother_lf"}
PARENT_HOME = {"father_lf", "mother_lf", "neither_lf", "lone_father_nilf", "lone_mother_nilf"}
WORKFORCE = [("nanny", "Childcare workers employed by private households"), ("centre", "Childcare workers in child day care services"),
             ("elsewhere", "Childcare workers elsewhere (schools, religious and recreation organisations)"), ("preschool_teacher", "Preschool and kindergarten teachers in child day care services")]
RACES = ["White", "Black", "Hispanic", "Asian", "Other"]
HH_TYPES = [("under6", "A child under 6"), ("only6_17", "Children 6 to 17 only"), ("none", "No children")]
# Published parameters (sources.yaml: dccny_market_rates_2024 / ocfs_market_rates_2024_attachment_a; weekly rates, county group 4, effective 1 October 2024).
RATES_WEEKLY = {"centre_infant": 500, "centre_toddler": 450, "centre_preschool": 439, "centre_school_age": 325,
                "home_under2": 400, "home_2": 350, "home_preschool": 325, "home_school_age": 289}
WEEKS = 52
COST_CASES = [("centre_infant", "Infant in a day care centre ($500 a week)"), ("home_under2", "Infant in a family day care home ($400 a week)"),
              ("centre_toddler", "Toddler in a day care centre ($450 a week)"), ("centre_preschool", "Three- or four-year-old in a centre ($439 a week)")]
COST_THRESHOLDS = [("t7", "Within 7% of family income (the conventional benchmark)", 0.07), ("t10", "Within 10% of income", 0.10), ("t20", "Within 20% of income", 0.20), ("t30", "Within 30% of income", 0.30)]
# Live births, New York City, from Table PC1 of the Summary of Vital Statistics 2023 (sources.yaml: nyc_summary_of_vital_statistics_2023), transcribed.
BIRTHS = {2010: 124_791, 2011: 123_029, 2012: 123_231, 2013: 120_457, 2014: 122_084, 2015: 121_673, 2016: 120_367, 2017: 117_013, 2018: 114_296, 2019: 110_442,
          2020: 100_022, 2021: 99_262, 2022: 99_459, 2023: 98_389}


# ------------------------------------------------------------------ loading


PERSON_COLS = ["SERIALNO", "SPORDER", "PUMA", "STATE", "PWGTP", "AGEP", "SEX", "RAC1P", "HISP", "NATIVITY", "CIT", "SCHL", "SCH", "SCHG", "ESR", "ESP", "RELSHIPP",
               "WAGP", "PERNP", "PINCP", "WKHP", "WKWN", "OCCP", "NAICSP", "COW", "FER", "MIG", "MIGSP", "MIGPUMA", "POVPIP", "ADJINC", "MAR"]
HH_COLS = ["SERIALNO", "HINCP", "TEN", "NP", "TYPEHUGQ", "HUPAC", "HUPAOC", "NOC", "R18", "GRPIP", "OCPIP", "WGTP", "FINCP", "HHT", "BDSP", "PUMA", "STATE", "ADJINC"]


def child_band(a) -> np.ndarray:
    v = np.asarray(a, float)
    return np.select([v < 1, v < 2, v < 3, v < 4, v < 5, v < 6, v < 13, v < 18], ["a0", "a1", "a2", "a3", "a4", "a5", "a6_12", "a13_17"], "adult")


def leave_band(a) -> np.ndarray:
    v = np.asarray(a, float)
    return np.select([v < 1, v < 6, v < 13, v < 18, v < 25, v < 45], ["infant", "c1_5", "c6_12", "c13_17", "adults_18_24", "adults_25_44"], "adults_45_plus")


def _num(s):
    return pd.to_numeric(s, errors="coerce")


def _derive_person(h: pd.DataFrame, year: int, defl: float) -> pd.DataFrame:
    h["child"] = h["AGEP"] < 18
    h["age_band"] = child_band(h["AGEP"]); h["leave_band"] = leave_band(h["AGEP"])
    esp = _num(h["ESP"])
    h["esp"] = esp.map(ESP_CODES)
    gq = _num(h["TYPEHUGQ"]) != 1
    h["own_child"] = h["child"] & esp.notna()
    h["arrangement"] = np.where(h["own_child"], h["esp"].map(ARR_GROUP), np.where(h["child"] & ~gq, "no_parent", np.where(h["child"], "gq", "")))
    h["all_work"] = h["esp"].isin(ALL_WORK); h["parent_home"] = h["esp"].isin(PARENT_HOME)
    h["two_parents"] = esp.isin([1, 2, 3, 4]); h["one_parent"] = esp.isin([5, 6, 7, 8])
    h["race"] = race_label(h["RAC1P"], h["HISP"])
    h["foreign_born"] = _num(h["NATIVITY"]) == 2
    sch = _num(h["SCH"])
    h["enrolled"] = sch.isin([2, 3]); h["public_school"] = sch == 2; h["private_school"] = sch == 3
    h["age34"] = h["AGEP"].between(3, 4)
    adj = _num(h["ADJINC"]) / 1e6; hinc = _num(h["HINCP"])
    h["household"] = ~gq & hinc.notna()
    h["hh_income"] = np.where(h["household"], hinc * adj * defl, np.nan)
    h["income_band"] = income_band(h["hh_income"])
    pinc = _num(h["PINCP"]) * adj * defl
    h["other_income"] = np.where(h["household"], h["hh_income"] - pinc.fillna(0), np.nan)
    h["other_income_band"] = income_band(h["other_income"])
    ten = _num(h["TEN"]); h["owner"] = ten.isin([1, 2]); h["renter"] = ten.isin([3, 4])
    pov = _num(h["POVPIP"]); h["poverty"] = pov < 100; h["poverty_known"] = pov.notna()
    grp = _num(h["GRPIP"]); h["rent_burden"] = h["renter"] & (grp >= 30); h["rent_ratio_known"] = h["renter"] & grp.notna() & (grp < 101)
    esr = _num(h["ESR"]); h["in_lf"] = esr.isin([1, 2, 3, 4, 5]); h["employed"] = esr.isin([1, 2, 4, 5])
    h["fulltime"] = h["employed"] & (_num(h["WKHP"]) >= 35)
    h["adult16"] = h["AGEP"] >= 16
    # Parents: householder, spouse or partner living with a child of the householder under 18.
    rel = _num(h["RELSHIPP"])
    kids = h.loc[rel.isin([25, 26, 27]) & h["child"], ["SERIALNO", "AGEP"]].groupby("SERIALNO")["AGEP"].min()
    h["youngest"] = h["SERIALNO"].map(kids)
    h["parent"] = rel.isin([20, 21, 22, 23, 24]) & h["youngest"].notna() & h["adult16"]
    h["youngest_band"] = np.where(h["parent"], child_band(h["youngest"].fillna(99)), "")
    h["mother"] = h["parent"] & (_num(h["SEX"]) == 2); h["father"] = h["parent"] & (_num(h["SEX"]) == 1)
    h["new_mother"] = _num(h["FER"]) == 1
    # The child-care workforce.
    occ = _num(h["OCCP"]); ind = h["NAICSP"].astype(str).str.strip()
    cc = h["employed"] & (occ == 4600)
    h["workforce"] = np.select([cc & (ind == "814"), cc & (ind == "6244"), cc, h["employed"] & (occ == 2300) & (ind == "6244")], ["nanny", "centre", "elsewhere", "preschool_teacher"], "")
    earn = _num(h["PERNP"]) * adj * defl; hrs = _num(h["WKHP"]); wks = _num(h["WKWN"])
    ok = (earn > 0) & (hrs > 0) & (wks > 0)
    h["earnings"] = np.where(earn > 0, earn, np.nan)
    h["hourly"] = np.where(ok, earn / (hrs * wks), np.nan)
    h["hours"] = np.where(hrs > 0, hrs, np.nan)
    h["self_employed"] = _num(h["COW"]).isin([6, 7])
    h["female"] = _num(h["SEX"]) == 2
    h["year"] = year
    return h


def load_year(year: int) -> tuple[pd.DataFrame, np.ndarray, pd.DataFrame, np.ndarray]:
    """Residents of the five boroughs stacked with the leavers found in the national
    file (person replicate weights in the same order), and the city's households with
    their own replicate weights."""
    v = vintage(year); geo = _geo(v); codes = NYC_MIGPUMA[v]; defl = cpi_deflators().get(year, 1.0)
    p = pd.read_csv(INTERIM / "acs1" / f"metro_person_{year}.csv.gz", usecols=lambda c: c in PERSON_COLS, dtype={"STATE": str, "PUMA": str, "SERIALNO": str, "NAICSP": str, "OCCP": str})
    rw = np.load(INTERIM / "acs1" / f"metro_person_{year}_repwts.npy", mmap_mode="r")
    assert len(rw) == len(p)
    p["puma_geoid"] = p["STATE"].str.zfill(2) + p["PUMA"].str.zfill(5)
    p = p.merge(geo, on="puma_geoid", how="left")
    keep = (p["in_nyc"] == 1).to_numpy()
    r = p[keep].reset_index(drop=True); rw_r = np.asarray(rw[keep], dtype=np.float32)
    hh_all = pd.read_csv(INTERIM / "acs1" / f"metro_housing_{year}.csv.gz", usecols=lambda c: c in HH_COLS, dtype={"SERIALNO": str, "STATE": str, "PUMA": str})
    hrw = np.load(INTERIM / "acs1" / f"metro_housing_{year}_repwts.npy", mmap_mode="r")
    assert len(hrw) == len(hh_all)
    r = r.merge(hh_all.drop(columns=["PUMA", "STATE", "ADJINC", "WGTP"]), on="SERIALNO", how="left")
    mig = _num(r["MIG"]); sp = _num(r["MIGSP"]); mp = _num(r["MIGPUMA"])
    us_state = sp.between(1, 56)
    r["status"] = np.select([mig.isna(), mig == 1, (mig == 3) & (sp == 36) & mp.isin(list(codes)), (mig == 2) | ((mig == 3) & ~us_state)], ["infant", "stayed", "within", "arrived_abroad"], "arrived_domestic")
    r["destination"] = ""; r["dest_region"] = ""; r["now_owner"] = np.nan
    r = _derive_person(r, year, defl)

    # Leavers (post-move household attributes).
    l = pd.read_csv(MIG_DIR / f"leavers_{year}.csv.gz", dtype={"STATE": str, "PUMA": str, "SERIALNO": str, "POWPUMA": str, "OCCP": str, "NAICSP": str})
    rw_l = np.load(MIG_DIR / f"leavers_{year}_repwts.npy")
    assert len(rw_l) == len(l)
    keep = (~l["now_nyc"]).to_numpy()
    l = l[keep].reset_index(drop=True); rw_l = np.asarray(rw_l[keep], dtype=np.float32)
    lhh = pd.read_csv(MIG_DIR / f"leavers_{year}_hh.csv.gz", dtype={"SERIALNO": str})
    l = l.merge(lhh[[c for c in lhh.columns if c in HH_COLS and c not in ("ADJINC", "WGTP")]], on="SERIALNO", how="left")
    l["puma_geoid"] = l["STATE"].str.zfill(2) + l["PUMA"].str.zfill(5)
    l = l.merge(geo, on="puma_geoid", how="left")
    st = _num(l["STATE"]); suburb = l["region"].isin(list(SUBURB_REGIONS))
    l["status"] = "left"; l["destination"] = np.where(suburb, "suburbs", state_group(st)); l["dest_region"] = np.where(suburb, l["region"].fillna(""), "")
    l["borough"] = l["origin_borough"]
    for c in PERSON_COLS:
        if c not in l.columns: l[c] = np.nan
    for c in HH_COLS:
        if c not in l.columns: l[c] = np.nan
    l = _derive_person(l, year, defl)
    cols = sorted(set(r.columns) & set(l.columns))
    h = pd.concat([r[cols], l[cols]], ignore_index=True)
    h["resident"] = h["status"] != "left"

    # Households of the five boroughs.
    hh_all["puma_geoid"] = hh_all["STATE"].str.zfill(2) + hh_all["PUMA"].str.zfill(5)
    hh_all = hh_all.merge(geo, on="puma_geoid", how="left")
    hk = ((hh_all["in_nyc"] == 1) & (_num(hh_all["TYPEHUGQ"]) == 1) & (_num(hh_all["NP"]) > 0)).to_numpy()
    hh = hh_all[hk].reset_index(drop=True); hrw_k = np.asarray(hrw[hk], dtype=np.float32)
    adj = _num(hh["ADJINC"]) / 1e6
    hh["hh_income"] = _num(hh["HINCP"]) * adj * defl; hh["income_band"] = income_band(hh["hh_income"])
    hh["fam_income"] = _num(hh["FINCP"]) * adj * defl
    hup = _num(hh["HUPAC"]); hh["hh_type"] = np.select([hup.isin([1, 3]), hup == 2], ["under6", "only6_17"], "none")
    ten = _num(hh["TEN"]); hh["renter"] = ten.isin([3, 4]); hh["owner"] = ten.isin([1, 2])
    grp = _num(hh["GRPIP"]); hh["rent_burden"] = hh["renter"] & (grp >= 30); hh["rent_ratio_known"] = hh["renter"] & grp.notna() & (grp < 101)
    hh["crowded"] = _num(hh["NP"]) > 2 * _num(hh["BDSP"]).replace(0, 0.5)
    hh["youngest"] = hh["SERIALNO"].map(r.loc[r["parent"], ["SERIALNO", "youngest"]].drop_duplicates("SERIALNO").set_index("SERIALNO")["youngest"])
    hh["PWGTP"] = _num(hh["WGTP"]); hh["year"] = year
    return h, np.concatenate([rw_r, rw_l]), hh, hrw_k


# ------------------------------------------------------------------ published tables


CELLS = {"under6": 2, "u6_two_parents": 3, "u6_both_lf": 4, "u6_father_only_lf": 5, "u6_mother_only_lf": 6, "u6_neither_lf": 7, "u6_one_parent": 8,
         "u6_father_lf": 10, "u6_father_nilf": 11, "u6_mother_lf": 13, "u6_mother_nilf": 14, "age6_17": 15, "a6_17_both_lf": 17, "total": 1}
CELL_LABELS = {"total": "Own children under 18", "under6": "Own children under 6", "u6_two_parents": "Under 6, living with two parents", "u6_both_lf": "Under 6, two parents both in the labor force",
               "u6_father_only_lf": "Under 6, two parents, father only in the labor force", "u6_mother_only_lf": "Under 6, two parents, mother only in the labor force",
               "u6_neither_lf": "Under 6, two parents, neither in the labor force", "u6_one_parent": "Under 6, living with one parent", "u6_father_lf": "Under 6, with father, in the labor force",
               "u6_father_nilf": "Under 6, with father, not in the labor force", "u6_mother_lf": "Under 6, with mother, in the labor force", "u6_mother_nilf": "Under 6, with mother, not in the labor force",
               "age6_17": "Own children 6 to 17", "a6_17_both_lf": "6 to 17, two parents both in the labor force"}
ESP_OF_CELL = {"under6": None, "u6_two_parents": [1, 2, 3, 4], "u6_both_lf": [1], "u6_father_only_lf": [2], "u6_mother_only_lf": [3], "u6_neither_lf": [4], "u6_one_parent": [5, 6, 7, 8],
               "u6_father_lf": [5], "u6_father_nilf": [6], "u6_mother_lf": [7], "u6_mother_nilf": [8], "age6_17": None, "a6_17_both_lf": [1], "total": None}


def published_tables() -> dict:
    """B23008 as published for New York city (and the nation in 2024)."""
    out = {"table_note": "B23008: own children under 18 in families and subfamilies by living arrangement and parents' employment status; cell numbers as in the Bureau's table."}
    g = pd.read_csv(SF_2019 / "g20191ny.csv", header=None, dtype=str, encoding="latin-1"); logrec = dict(zip(g[48], g[4]))
    e = pd.read_csv(SF_2019 / "e20191ny0101000.txt", header=None, dtype=str, low_memory=False)
    m = pd.read_csv(SF_2019 / "m20191ny0101000.txt", header=None, dtype=str, low_memory=False)
    lr = logrec["16000US3651000"]; er = e[e[5] == lr].iloc[0]; mr = m[m[5] == lr].iloc[0]; start = 66 - 1
    out["2019"] = {k: {"estimate": float(er[start + c - 1]), "moe90": float(mr[start + c - 1])} for k, c in CELLS.items()}
    for y in (2021, 2022, 2023):
        d = pd.read_csv(REFERENCE / "tables" / f"acsdt1y{y}-b23008.dat", sep="|", dtype=str).set_index("GEO_ID")
        r = d.loc["1600000US3651000"]
        out[str(y)] = {k: {"estimate": float(r[f"B23008_E{c:03d}"]), "moe90": float(r[f"B23008_M{c:03d}"])} for k, c in CELLS.items()}
    j = json.loads(TABLE_2024.read_text(encoding="utf-8")); b = j["data"]["16000US3651000"]["B23008"]
    out["2024"] = {k: {"estimate": b["estimate"][f"B23008{c:03d}"], "moe90": b["error"][f"B23008{c:03d}"]} for k, c in CELLS.items()}
    bu = j["data"]["01000US"]["B23008"]
    out["2024_us"] = {k: {"estimate": bu["estimate"][f"B23008{c:03d}"], "moe90": bu["error"][f"B23008{c:03d}"]} for k, c in CELLS.items()}
    return out


def reproduction(pub: dict, censuses: dict) -> pd.DataFrame:
    rows = []
    for y in ("2019", "2021", "2022", "2023", "2024"):
        if y not in pub or int(y) not in censuses: continue
        cs = censuses[int(y)]; h = cs.h; res = h["resident"].to_numpy(); own = (h["own_child"] & res).to_numpy(); esp = _num(h["ESP"]).to_numpy()
        u6 = own & (h["AGEP"] < 6).to_numpy(); o = own & (h["AGEP"] >= 6).to_numpy()
        for k, c in CELLS.items():
            base = {"total": own, "age6_17": o, "a6_17_both_lf": o}.get(k, u6)
            codes = ESP_OF_CELL[k]
            mask = base if codes is None else base & np.isin(esp, codes)
            v, vm = cs.count(mask); p = pub[y][k]; band = float(np.sqrt(p["moe90"] ** 2 + vm ** 2))
            rows.append({"statistic": f"{CELL_LABELS[k]}, NYC, {y}", "published": p["estimate"], "published_moe90": p["moe90"], "reproduced": v, "reproduced_moe90": vm,
                         "difference": v - p["estimate"], "difference_moe90": band, "within_margin": bool(abs(v - p["estimate"]) <= band), "kind": "count", "year": int(y), "cell": k})
    return pd.DataFrame(rows)


# ------------------------------------------------------------------ census estimates


def year_estimates(cs: Census, hs: Census) -> tuple[dict, list[dict]]:
    h = cs.h; hh = hs.h; year = int(h["year"].iloc[0]); rows = []
    res = h["resident"].to_numpy(); child = (h["child"] & res).to_numpy(); own = (h["own_child"] & res).to_numpy(); age = h["AGEP"].to_numpy(dtype=float)
    band = h["age_band"].to_numpy(); arr = h["arrangement"].to_numpy(); esp = h["esp"].fillna("").to_numpy(); boro = h["borough"].fillna("").to_numpy()
    u6 = child & (age < 6); u3 = child & (age < 3); hhold = h["household"].to_numpy(); status = h["status"].to_numpy(); left = status == "left"
    inc = h["income_band"].to_numpy(); race = h["race"].to_numpy()

    def put(key, label, value, moe, n, kind):
        rows.append({"year": year, "key": key, "label": label, "value": value, "moe90": moe, "n_records": int(n), "kind": kind})

    out: dict = {"year": year, "n_records": {"residents": int(res.sum()), "children": int(child.sum()), "households": int(len(hh)), "leavers": int(left.sum())}, "experimental_weights": year == 2020}
    c, cm = cs.count(child); out["children"] = {"total": {"count": c, "moe90": cm}}
    for k, lab, mask in [("under6", "Children under 6", u6), ("under3", "Children under 3", u3), ("age3_5", "Children 3 to 5", child & (age >= 3) & (age < 6)), ("age6_17", "Children 6 to 17", child & (age >= 6)),
                         ("own_under6", "Own children under 6 (the published table's universe)", own & (age < 6))]:
        c, cm = cs.count(mask); out["children"][k] = {"count": c, "moe90": cm, "n_records": int(mask.sum())}; put(f"children_{k}", lab, c, cm, mask.sum(), "count")
    out["children"]["by_age"] = {k: dict(zip(("count", "moe90"), cs.count(child & (band == k))), label=l) for k, l, *_ in AGE_BANDS}
    out["children"]["by_borough"] = {}
    for b in BOROUGHS:
        mb = boro == b
        out["children"]["by_borough"][b] = {"under6": cs.count(u6 & mb)[0], "under6_moe90": cs.count(u6 & mb)[1], "age6_17": cs.count(child & (age >= 6) & mb)[0], "all": cs.count(child & mb)[0],
                                            "residents": cs.count(res & mb)[0], "under6_share_of_residents": cs.share(res & mb, u6)[0], "n_records": int((u6 & mb).sum())}
    out["children"]["under6_share_of_residents"] = cs.share(res, u6)[0]

    # Arrangements by age of child (own children; "no_parent" among all household children).
    out["arrangement"] = {"by_age": {}, "under6": {}, "age6_17": {}}
    for k, l, *_ in AGE_BANDS:
        mk = own & (band == k); blk = {"label": l, "n_records": int(mk.sum())}
        for ak, al in ARRANGEMENTS:
            if ak == "no_parent": v, vm = cs.share(child & hhold & (band == k), arr == "no_parent")
            else: v, vm = cs.share(mk, arr == ak)
            blk[ak] = {"share": v, "moe90": vm}
        blk["all_work"] = dict(zip(("share", "moe90"), cs.share(mk, h["all_work"].to_numpy()))); blk["parent_home"] = dict(zip(("share", "moe90"), cs.share(mk, h["parent_home"].to_numpy())))
        blk["two_parents"] = dict(zip(("share", "moe90"), cs.share(mk, h["two_parents"].to_numpy())))
        blk["mother_home_two_parents"] = dict(zip(("share", "moe90"), cs.share(mk, esp == "father_lf"))); blk["father_home_two_parents"] = dict(zip(("share", "moe90"), cs.share(mk, esp == "mother_lf")))
        out["arrangement"]["by_age"][k] = blk
    for key, mask in (("under6", own & (age < 6)), ("age6_17", own & (age >= 6))):
        blk = {"n_records": int(mask.sum()), "count": cs.count(mask)[0]}
        for ak, al in ARRANGEMENTS[:-1]:
            v, vm = cs.share(mask, arr == ak); c, cm = cs.count(mask & (arr == ak)); blk[ak] = {"share": v, "moe90": vm, "count": c, "count_moe90": cm, "label": al}
            put(f"arr_{key}_{ak}", f"{al}, {'under 6' if key == 'under6' else '6 to 17'}", v, vm, mask.sum(), "share")
        for ek, el in (("father_lf", "Two parents, mother not in the labor force"), ("mother_lf", "Two parents, father not in the labor force"), ("lone_mother_nilf", "Mother only, not in the labor force"),
                       ("lone_father_nilf", "Father only, not in the labor force"), ("lone_mother_lf", "Mother only, in the labor force"), ("lone_father_lf", "Father only, in the labor force")):
            v, vm = cs.share(mask, esp == ek); c, cm = cs.count(mask & (esp == ek)); blk[ek] = {"share": v, "moe90": vm, "count": c, "count_moe90": cm, "label": el}
        blk["all_work"] = dict(zip(("share", "moe90"), cs.share(mask, h["all_work"].to_numpy()))); blk["parent_home"] = dict(zip(("share", "moe90"), cs.share(mask, h["parent_home"].to_numpy())))
        base_all = child & hhold & ((age < 6) if key == "under6" else (age >= 6))
        blk["no_parent"] = dict(zip(("share", "moe90"), cs.share(base_all, arr == "no_parent"))); blk["no_parent_count"] = cs.count(base_all & (arr == "no_parent"))[0]
        out["arrangement"][key] = blk
    # Arrangement of under-sixes by household income band and borough.
    out["arrangement"]["under6_by_income"] = {}
    for k, l, *_ in INCOME_BANDS:
        mk = own & (age < 6) & (inc == k)
        out["arrangement"]["under6_by_income"][k] = {"label": l, "n_records": int(mk.sum()), "parent_home": cs.share(mk, h["parent_home"].to_numpy())[0], "parent_home_moe90": cs.share(mk, h["parent_home"].to_numpy())[1],
                                                     "all_work": cs.share(mk, h["all_work"].to_numpy())[0], "one_parent": cs.share(mk, h["one_parent"].to_numpy())[0], "count": cs.count(mk)[0]}
    out["arrangement"]["under6_by_borough"] = {}
    for b in BOROUGHS:
        mk = own & (age < 6) & (boro == b)
        out["arrangement"]["under6_by_borough"][b] = {"n_records": int(mk.sum()), "parent_home": cs.share(mk, h["parent_home"].to_numpy())[0], "all_work": cs.share(mk, h["all_work"].to_numpy())[0],
                                                      "one_parent": cs.share(mk, h["one_parent"].to_numpy())[0]}

    # Children's households: income, tenure, poverty, rent burden (child-weighted).
    out["children_households"] = {}
    for key, mask in (("under6", u6 & hhold), ("age6_17", child & (age >= 6) & hhold), ("all", child & hhold)):
        blk = {"n_records": int(mask.sum()), "median_hh_income": cs.median(mask, "hh_income"), "renter": cs.share(mask, h["renter"].to_numpy()), "owner": cs.share(mask, h["owner"].to_numpy()),
               "poverty": cs.share(mask & h["poverty_known"].to_numpy(), h["poverty"].to_numpy()), "rent_burden": cs.share(mask & h["rent_ratio_known"].to_numpy(), h["rent_burden"].to_numpy()),
               "income_lt50k": cs.share(mask, np.isin(inc, ["lt25k", "25k_50k"])), "income_150k_plus": cs.share(mask, np.isin(inc, ["150k_200k", "200k_300k", "300k_plus"])),
               "foreign_born_share": cs.share(mask, h["foreign_born"].to_numpy()), "income_bands": {k: cs.share(mask, inc == k)[0] for k, *_ in INCOME_BANDS}}
        blk["by_borough"] = {b: {"median_hh_income": cs.median(mask & (boro == b), "hh_income"), "renter": cs.share(mask & (boro == b), h["renter"].to_numpy())[0],
                                 "poverty": cs.share(mask & (boro == b) & h["poverty_known"].to_numpy(), h["poverty"].to_numpy())[0], "n_records": int((mask & (boro == b)).sum())} for b in BOROUGHS}
        out["children_households"][key] = blk
    for k, l, *_ in INCOME_BANDS:
        put(f"u6_income_{k}", f"Share of children under 6 in households {l}", *cs.share(u6 & hhold, inc == k), (u6 & hhold).sum(), "share")

    # Households (household-weighted, housing replicate weights).
    typ = hh["hh_type"].to_numpy(); hb = hh["borough"].fillna("").to_numpy(); hinc = hh["income_band"].to_numpy()
    out["households"] = {"total": hs.count(np.ones(len(hh), bool))[0]}
    for tk, tl in HH_TYPES:
        mk = typ == tk; c, cm = hs.count(mk)
        blk = {"label": tl, "count": c, "moe90": cm, "n_records": int(mk.sum()), "share_of_households": hs.share(np.ones(len(hh), bool), mk)[0], "median_hh_income": hs.median(mk, "hh_income"),
               "renter": hs.share(mk, hh["renter"].to_numpy()), "rent_burden": hs.share(mk & hh["rent_ratio_known"].to_numpy(), hh["rent_burden"].to_numpy()),
               "income_lt50k": hs.share(mk, np.isin(hinc, ["lt25k", "25k_50k"])), "income_150k_plus": hs.share(mk, np.isin(hinc, ["150k_200k", "200k_300k", "300k_plus"])),
               "crowded": hs.share(mk, hh["crowded"].to_numpy()), "income_bands": {k: hs.share(mk, hinc == k)[0] for k, *_ in INCOME_BANDS},
               "by_borough": {b: {"count": hs.count(mk & (hb == b))[0], "moe90": hs.count(mk & (hb == b))[1], "share_of_borough_households": hs.share(hb == b, mk)[0], "median_hh_income": hs.median(mk & (hb == b), "hh_income"),
                                  "renter": hs.share(mk & (hb == b), hh["renter"].to_numpy())[0], "n_records": int((mk & (hb == b)).sum())} for b in BOROUGHS}}
        out["households"][tk] = blk
        put(f"hh_{tk}", f"Households with {tl.lower()}", c, cm, mk.sum(), "count")
    # Households whose youngest child is under three: the cost-of-care universe. Family
    # income (FINCP) is the income the 7% benchmark is defined on; household income is a
    # sensitivity, since a household can hold people who are not in the child's family.
    yng = hh["youngest"].to_numpy(dtype=float); u3h = yng < 3; fi = hh["fam_income"].to_numpy(dtype=float); hi = hh["hh_income"].to_numpy(dtype=float)
    out["cost"] = {"universe": "households where the youngest child of the householder is under 3", "income_measure": "family income (FINCP), 2024 dollars", "n_records": int(u3h.sum()),
                   "households": hs.count(u3h)[0], "median_family_income": hs.median(u3h, "fam_income"), "median_hh_income": hs.median(u3h, "hh_income"),
                   "weeks": WEEKS, "rates_weekly": RATES_WEEKLY, "annual": {k: v * WEEKS for k, v in RATES_WEEKLY.items()}, "cases": {}, "sensitivity_household_income": {}}
    for ck, cl in COST_CASES:
        annual = RATES_WEEKLY[ck] * WEEKS; case = {"label": cl, "annual": annual, "income_needed_at_7pct": annual / 0.07, "thresholds": {}}; sens = {}
        for tk, tl, frac in COST_THRESHOLDS:
            ok = u3h & (fi > 0); v, vm = hs.share(ok, annual <= frac * np.nan_to_num(fi, nan=-1))
            case["thresholds"][tk] = {"label": tl, "share": v, "moe90": vm}
            put(f"cost_{ck}_{tk}", f"{cl}: share of households with a child under 3 for whom it is {tl.lower()} (family income)", v, vm, ok.sum(), "share")
            okh = u3h & (hi > 0); sens[tk] = hs.share(okh, annual <= frac * np.nan_to_num(hi, nan=-1))[0]
        case["share_of_median_income"] = annual / out["cost"]["median_family_income"] if out["cost"]["median_family_income"] else float("nan")
        out["cost"]["cases"][ck] = case; out["cost"]["sensitivity_household_income"][ck] = sens
    out["cost"]["by_borough"] = {b: {"median_family_income": hs.median(u3h & (hb == b), "fam_income"), "t7_centre_infant": hs.share(u3h & (hb == b) & (fi > 0), RATES_WEEKLY["centre_infant"] * WEEKS <= 0.07 * np.nan_to_num(fi, nan=-1))[0],
                                     "t20_centre_infant": hs.share(u3h & (hb == b) & (fi > 0), RATES_WEEKLY["centre_infant"] * WEEKS <= 0.20 * np.nan_to_num(fi, nan=-1))[0], "n_records": int((u3h & (hb == b)).sum())} for b in BOROUGHS}

    # Parents by the youngest child's age.
    par = (h["parent"] & res).to_numpy(); yb = h["youngest_band"].to_numpy(); mother = h["mother"].to_numpy(); father = h["father"].to_numpy()
    out["parents"] = {"by_youngest": {}, "n_records": int(par.sum()), "count": cs.count(par)[0], "mothers": cs.count(par & mother)[0], "fathers": cs.count(par & father)[0]}
    for k, l, *_ in YOUNGEST_BANDS:
        mk = par & (yb == k); blk = {"label": l, "n_records": int(mk.sum()), "count": cs.count(mk)[0]}
        for sk, sm in (("mothers", mother), ("fathers", father)):
            m2 = mk & sm; blk[sk] = {"n_records": int(m2.sum()), "count": cs.count(m2)[0]}
            for stat in ("employed", "in_lf", "fulltime"):
                v, vm = cs.share(m2, h[stat].to_numpy()); blk[sk][stat] = {"share": v, "moe90": vm}
                put(f"{sk}_{stat}_youngest_{k}", f"{sk.title()} {stat.replace('_', ' ')}, youngest child {l.lower()}", v, vm, m2.sum(), "share")
            blk[sk]["not_in_lf_count"] = cs.count(m2 & ~h["in_lf"].to_numpy())[0]
        out["parents"]["by_youngest"][k] = blk
    yu6 = par & np.isin(yb, ["a0", "a1", "a2", "a3", "a4", "a5"])
    out["parents"]["youngest_under6"] = {"count": cs.count(yu6)[0], "n_records": int(yu6.sum()),
                                         "mothers": {"count": cs.count(yu6 & mother)[0], "employed": cs.share(yu6 & mother, h["employed"].to_numpy()), "in_lf": cs.share(yu6 & mother, h["in_lf"].to_numpy()),
                                                     "fulltime": cs.share(yu6 & mother, h["fulltime"].to_numpy()), "not_in_lf_count": dict(zip(("count", "moe90"), cs.count(yu6 & mother & ~h["in_lf"].to_numpy())))},
                                         "fathers": {"count": cs.count(yu6 & father)[0], "employed": cs.share(yu6 & father, h["employed"].to_numpy()), "in_lf": cs.share(yu6 & father, h["in_lf"].to_numpy()),
                                                     "fulltime": cs.share(yu6 & father, h["fulltime"].to_numpy()), "not_in_lf_count": dict(zip(("count", "moe90"), cs.count(yu6 & father & ~h["in_lf"].to_numpy())))},
                                         "at_home_count": dict(zip(("count", "moe90"), cs.count(yu6 & ~h["in_lf"].to_numpy()))), "at_home_mother_share": cs.share(yu6 & ~h["in_lf"].to_numpy(), mother)[0]}
    # Parents at home by household income (youngest under 6): the salary question.
    out["parents"]["at_home_by_income"] = {}
    for k, l, *_ in INCOME_BANDS:
        for sk, sm in (("mothers", mother), ("fathers", father)):
            m2 = yu6 & sm & (inc == k); v, vm = cs.share(m2, ~h["in_lf"].to_numpy())
            out["parents"]["at_home_by_income"].setdefault(k, {"label": l})[sk] = {"share": v, "moe90": vm, "n_records": int(m2.sum()), "employed": cs.share(m2, h["employed"].to_numpy())[0]}
    oinc = h["other_income_band"].to_numpy(); out["parents"]["at_home_by_other_income"] = {}
    for k, l, *_ in INCOME_BANDS:
        for sk, sm in (("mothers", mother), ("fathers", father)):
            m2 = yu6 & sm & (oinc == k); v, vm = cs.share(m2, ~h["in_lf"].to_numpy())
            out["parents"]["at_home_by_other_income"].setdefault(k, {"label": l})[sk] = {"share": v, "moe90": vm, "n_records": int(m2.sum()), "employed": cs.share(m2, h["employed"].to_numpy())[0]}
    out["parents"]["at_home_by_borough"] = {b: {"mothers": cs.share(yu6 & mother & (boro == b), ~h["in_lf"].to_numpy())[0], "fathers": cs.share(yu6 & father & (boro == b), ~h["in_lf"].to_numpy())[0],
                                              "n_records": int((yu6 & (boro == b)).sum())} for b in BOROUGHS}
    out["parents"]["at_home_by_race"] = {rc: {"mothers": cs.share(yu6 & mother & (race == rc), ~h["in_lf"].to_numpy())[0], "n_records": int((yu6 & mother & (race == rc)).sum())} for rc in RACES}
    out["parents"]["at_home_by_nativity"] = {"us_born": cs.share(yu6 & mother & ~h["foreign_born"].to_numpy(), ~h["in_lf"].to_numpy())[0], "foreign_born": cs.share(yu6 & mother & h["foreign_born"].to_numpy(), ~h["in_lf"].to_numpy())[0]}
    nm = res & h["new_mother"].to_numpy()
    out["new_mothers"] = {"count": cs.count(nm)[0], "moe90": cs.count(nm)[1], "n_records": int(nm.sum()), "in_lf": dict(zip(("share", "moe90"), cs.share(nm, h["in_lf"].to_numpy()))),
                          "employed": dict(zip(("share", "moe90"), cs.share(nm, h["employed"].to_numpy()))), "fulltime": dict(zip(("share", "moe90"), cs.share(nm, h["fulltime"].to_numpy()))),
                          "median_age": cs.median(nm, "AGEP"), "by_income": {k: {"in_lf": cs.share(nm & (inc == k), h["in_lf"].to_numpy())[0], "n_records": int((nm & (inc == k)).sum())} for k, *_ in INCOME_BANDS}}

    # Preschool at three and four.
    out["preschool"] = {}
    for a in (3, 4, 5):
        mk = child & (age == a); blk = {"n_records": int(mk.sum()), "count": cs.count(mk)[0]}
        for stat in ("enrolled", "public_school", "private_school"):
            v, vm = cs.share(mk, h[stat].to_numpy()); blk[stat] = {"share": v, "moe90": vm}
            put(f"preschool_{stat}_age{a}", f"Age {a}: {stat.replace('_', ' ')}", v, vm, mk.sum(), "share")
        blk["by_income"] = {k: {"enrolled": cs.share(mk & (inc == k), h["enrolled"].to_numpy())[0], "public_school": cs.share(mk & (inc == k), h["public_school"].to_numpy())[0], "n_records": int((mk & (inc == k)).sum())} for k, *_ in INCOME_BANDS}
        blk["by_borough"] = {b: {"enrolled": cs.share(mk & (boro == b), h["enrolled"].to_numpy())[0], "public_school": cs.share(mk & (boro == b), h["public_school"].to_numpy())[0], "n_records": int((mk & (boro == b)).sum())} for b in BOROUGHS}
        blk["by_arrangement"] = {"all_work": cs.share(mk & h["all_work"].to_numpy(), h["enrolled"].to_numpy())[0], "parent_home": cs.share(mk & h["parent_home"].to_numpy(), h["enrolled"].to_numpy())[0]}
        out["preschool"][f"age{a}"] = blk

    # The child-care workforce.
    wf = h["workforce"].to_numpy(); out["workforce"] = {}
    for wk, wl in WORKFORCE:
        mk = res & (wf == wk); c, cm = cs.count(mk)
        blk = {"label": wl, "count": c, "moe90": cm, "n_records": int(mk.sum()), "median_hourly": cs.median(mk, "hourly"), "median_earnings": cs.median(mk, "earnings"), "median_hours": cs.median(mk, "hours"),
               "foreign_born": cs.share(mk, h["foreign_born"].to_numpy()), "female": cs.share(mk, h["female"].to_numpy()), "self_employed": cs.share(mk, h["self_employed"].to_numpy()),
               "fulltime": cs.share(mk, h["fulltime"].to_numpy()), "median_age": cs.median(mk, "AGEP"), "noncitizen": cs.share(mk, (_num(h["CIT"]) == 5).to_numpy()),
               "by_borough": {b: cs.count(mk & (boro == b))[0] for b in BOROUGHS}, "poverty": cs.share(mk & h["poverty_known"].to_numpy(), h["poverty"].to_numpy())[0],
               "lives_with_children": cs.share(mk & hhold, (_num(h["R18"]) == 1).to_numpy())[0], "median_hh_income": cs.median(mk & hhold, "hh_income")}
        out["workforce"][wk] = blk; put(f"workforce_{wk}", wl, c, cm, mk.sum(), "count")
    allcc = res & np.isin(wf, ["nanny", "centre", "elsewhere"]); out["workforce"]["all_childcare_workers"] = {"count": cs.count(allcc)[0], "moe90": cs.count(allcc)[1], "n_records": int(allcc.sum()), "median_hourly": cs.median(allcc, "hourly")}
    allw = res & h["employed"].to_numpy(); out["workforce"]["all_workers_median_hourly"] = cs.median(allw, "hourly")

    # Leaving: children and parents.
    lived_here = np.isin(status, ["stayed", "within", "left"]); lb = h["leave_band"].to_numpy(); dest = h["destination"].to_numpy(); dreg = h["dest_region"].to_numpy()
    out["leaving"] = {"by_band": {}, "children_1_17": {}, "under6_net": {}, "by_income_now": {}, "by_borough": {}, "destinations": {}, "by_race": {}}
    for k, l, *_ in LEAVE_BANDS:
        mk = lived_here & (lb == k); v, vm = cs.share(mk, left); c, cm = cs.count(mk & left)
        out["leaving"]["by_band"][k] = {"label": l, "out_rate": v, "moe90": vm, "left": c, "left_moe90": cm, "n_records": int(mk.sum()), "n_left": int((mk & left).sum())}
        put(f"out_rate_{k}", f"Out-migration rate, {l.lower()}", v, vm, mk.sum(), "share")
    kids = lived_here & (age >= 1) & (age < 18); v, vm = cs.share(kids, left)
    out["leaving"]["children_1_17"] = {"out_rate": v, "moe90": vm, "left": cs.count(kids & left)[0], "left_moe90": cs.count(kids & left)[1], "n_left": int((kids & left).sum())}
    adults = lived_here & (age >= 18); out["leaving"]["adults"] = dict(zip(("out_rate", "moe90"), cs.share(adults, left)))
    # With children now: adults 25-44 by whether their household now includes a child.
    r18 = (_num(h["R18"]) == 1).to_numpy(); a2544 = lived_here & (age >= 25) & (age < 45) & hhold
    out["leaving"]["adults_25_44_with_children_now"] = dict(zip(("out_rate", "moe90"), cs.share(a2544 & r18, left)))
    out["leaving"]["adults_25_44_no_children_now"] = dict(zip(("out_rate", "moe90"), cs.share(a2544 & ~r18, left)))
    # Net domestic flow of children (arrivals from elsewhere in the US minus leavers).
    for key, mask in (("under6", (age >= 1) & (age < 6)), ("age6_17", (age >= 6) & (age < 18)), ("all_1_17", (age >= 1) & (age < 18))):
        n_, nm_ = cs.diff(mask & (status == "arrived_domestic"), mask & left)
        out["leaving"]["under6_net"][key] = {"net": n_, "moe90": nm_, "arrived": cs.count(mask & (status == "arrived_domestic"))[0], "arrived_abroad": cs.count(mask & (status == "arrived_abroad"))[0], "left": cs.count(mask & left)[0]}
    yc = lived_here & (age >= 1) & (age < 6)
    for k, l, *_ in INCOME_BANDS:
        mk = yc & (inc == k); v, vm = cs.share(mk, left)
        out["leaving"]["by_income_now"][k] = {"label": l, "out_rate": v, "moe90": vm, "n_records": int(mk.sum()), "n_left": int((mk & left).sum())}
    for b in BOROUGHS:
        mk = yc & (boro == b); v, vm = cs.share(mk, left)
        out["leaving"]["by_borough"][b] = {"out_rate": v, "moe90": vm, "n_records": int(mk.sum()), "n_left": int((mk & left).sum()), "left": cs.count(mk & left)[0]}
    for rc in RACES:
        mk = yc & (race == rc); v, vm = cs.share(mk, left); out["leaving"]["by_race"][rc] = {"out_rate": v, "moe90": vm, "n_records": int(mk.sum()), "n_left": int((mk & left).sum())}
    lk = left & (age >= 1) & (age < 6); tot = cs.count(lk)[0]
    for dk, dl in [("suburbs", "The metro suburbs"), ("rest_ny", "Elsewhere in New York State"), ("rest_nj_ct", "Elsewhere in New Jersey or Connecticut"), ("florida", "Florida"), ("pennsylvania", "Pennsylvania"),
                   ("other_northeast", "The rest of the Northeast"), ("california", "California"), ("texas", "Texas"), ("other_south", "The rest of the South"), ("midwest", "The Midwest"), ("other_west", "The rest of the West")]:
        c, cm = cs.count(lk & (dest == dk)); out["leaving"]["destinations"][dk] = {"label": dl, "count": c, "moe90": cm, "share": c / tot if tot else float("nan"), "n_records": int((lk & (dest == dk)).sum())}
    out["leaving"]["suburb_regions"] = {k: {"label": l, "count": cs.count(lk & (dreg == k))[0]} for k, l in SUBURB_REGIONS.items()}
    out["leaving"]["young_leavers_now"] = {"n_records": int(lk.sum()), "owner": cs.share(lk & hhold, h["owner"].to_numpy())[0], "median_hh_income": cs.median(lk & hhold, "hh_income"),
                                           "income_150k_plus": cs.share(lk & hhold, np.isin(inc, ["150k_200k", "200k_300k", "300k_plus"]))[0], "income_lt50k": cs.share(lk & hhold, np.isin(inc, ["lt25k", "25k_50k"]))[0]}
    stay_young = np.isin(status, ["stayed", "within"]) & (age >= 1) & (age < 6) & hhold
    out["leaving"]["young_stayers_now"] = {"owner": cs.share(stay_young, h["owner"].to_numpy())[0], "median_hh_income": cs.median(stay_young, "hh_income"),
                                           "income_150k_plus": cs.share(stay_young, np.isin(inc, ["150k_200k", "200k_300k", "300k_plus"]))[0], "income_lt50k": cs.share(stay_young, np.isin(inc, ["lt25k", "25k_50k"]))[0]}
    return out, rows


# ------------------------------------------------------------------ administrative series


def under5_estimates() -> dict:
    """Population under five by borough, July to July, two vintages."""
    def load(path, year_of, vint):
        d = pd.read_csv(path, encoding="latin-1", dtype={"STATE": str, "COUNTY": str})
        n = d[(d["STATE"] == "36") & d["COUNTY"].isin(NYC_COUNTIES) & (d["AGEGRP"] == 1)]
        rows = []
        for code, year in year_of.items():
            s = n[n["YEAR"] == code]
            rows.append({"year": year, "vintage": vint, "total": int(s["TOT_POP"].sum()), "boroughs": {NYC_COUNTIES[c]: int(v) for c, v in zip(s["COUNTY"], s["TOT_POP"])}})
        return rows
    v20 = load(RAW / "popest" / "cc-est2020-alldata-36.csv", {y: 2007 + y for y in range(3, 14)}, "2020")          # YEAR 3 = July 2010 ... 13 = July 2020
    v25 = load(RAW / "popest" / "cc-est2025-alldata-36.csv", {y: 2018 + y for y in range(2, 8)}, "2025")           # YEAR 2 = July 2020 ... 7 = July 2025
    base = load(RAW / "popest" / "cc-est2025-alldata-36.csv", {1: 2020}, "2025 (April 2020 base)")[0]
    def chg(rows, a, b):
        ra = next(r for r in rows if r["year"] == a); rb = next(r for r in rows if r["year"] == b)
        return {"from": a, "to": b, "total": rb["total"] - ra["total"], "pct": rb["total"] / ra["total"] - 1, "boroughs": {k: rb["boroughs"][k] / ra["boroughs"][k] - 1 for k in ra["boroughs"]}}
    return {"note": "Residents under five, five counties, 1 July of each year. 2010-2020 from Vintage 2020 (built on the 2010 census), 2020-2025 from Vintage 2025 (built on the 2020 census); "
                    "the two vintages do not splice and no change is computed across them. The April 2020 census base is given separately.", "years": v20 + v25, "april_2020_base": base,
            "within_vintage_changes": {"v2020_2015_2020": chg(v20, 2015, 2020), "v2025_2020_2025": chg(v25, 2020, 2025), "v2025_2022_2025": chg(v25, 2022, 2025)}}


def school_enrolment() -> dict:
    out = {"note": "New York City public school enrolment by grade, audited register, from the Department of Education's demographic snapshot. Kindergarten is shown with and without charter schools; 3-K and pre-K are district programs (community-based 3-K and pre-K seats run by the DOE are included in its counts).", "years": []}
    x = pd.ExcelFile(DOE_SNAPSHOT)
    city = pd.read_excel(x, "Citywide"); charter = pd.read_excel(x, "Citywide with Charter Schools")
    for (_, a), (_, b) in zip(city.iterrows(), charter.iterrows()):
        out["years"].append({"year": str(a["Year"]), "three_k": int(a["Grade 3K"]), "pre_k": int(a["Grade PK (Half Day & Full Day)"]), "kindergarten": int(a["Grade K"]), "kindergarten_with_charters": int(b["Grade K"]),
                             "grade1_with_charters": int(b["Grade 1"]), "total": int(a["Total Enrollment"]), "total_with_charters": int(b["Total Enrollment"])})
    return out


# ------------------------------------------------------------------ figures (static twins)


def fig_from_cube(spec: dict, path: Path, height: float | None = None) -> None:
    """Draw the chart's default state as a grouped horizontal bar chart."""
    state = "|".join(str(c.get("default", c["options"][0][0])) for c in spec.get("controls", []))
    row = spec["cells"].get(state) or spec["cells"].get("") or {}
    mm = spec["metrics"][0]; fmt = {"pct": "{:.0f}%", "dollar": "${:,.0f}", "count": "{:,.0f}"}.get(mm["format"], "{:,.0f}")
    series = [s for s in spec["series"] if s[0] in row] if spec.get("series_by_cell") else spec["series"]
    if spec.get("paired"):
        labels = [p[1] for p in spec["paired"]]
        rows = [(s[1], [(row.get(s[0], {}).get(p[0]), row.get(s[0], {}).get(p[0] + "_moe")) for p in spec["paired"]]) for s in series]
    else:
        labels = [mm["label"]]
        rows = [(s[1], [(row.get(s[0], {}).get(mm["key"]), row.get(s[0], {}).get(mm["key"] + "_moe"))]) for s in series]
    rows = [(l, [(v if v is not None else float("nan"), m if m is not None else float("nan")) for v, m in vals]) for l, vals in rows]
    if mm["format"] == "count" and any(abs(v) >= 10_000 for _, vals in rows for v, _ in vals if v == v):
        rows = [(l, [(v / 1000 if v == v else v, m / 1000 if m == m else 0) for v, m in vals]) for l, vals in rows]; fmt = "{:,.0f}k"
    _grouped(rows, path, spec["title"], spec["subtitle"], labels, fmt=fmt, height=height or (1.6 + 0.3 * len(labels) * len(rows)), xlabel=mm.get("axis"))


def fig_under5(u5: dict, path: Path):
    plt = _mpl(); ys = u5["years"]
    fig, ax = plt.subplots(figsize=(FIG_W, 3.4), dpi=100); fig.patch.set_alpha(0)
    x = np.arange(len(ys)); ax.bar(x, [r["total"] / 1000 for r in ys], color=[BLUE if r["vintage"] == "2020" else GREEN for r in ys], width=0.8)
    ax.axvline(10.5, color=GRID, linewidth=1, linestyle=":"); ax.text(10.6, 590, "Vintage 2025 →", fontsize=8, color=INK); ax.text(10.4, 590, "← Vintage 2020", fontsize=8, color=INK, ha="right")
    ax.set_xticks(x); ax.set_xticklabels([str(r["year"]) for r in ys], fontsize=8, rotation=45)
    ax.set_ylim(0, 620); ax.set_yticks([0, 200, 400, 600]); ax.set_yticklabels(["0", "200k", "400k", "600k"])
    ax.grid(axis="y", color=GRID, linewidth=0.8); ax.grid(axis="x", visible=False)
    for s in ("top", "right"): ax.spines[s].set_visible(False)
    ax.set_title("Children under five, July of each year", loc="left", fontsize=12, color="#1a1a1a", pad=22)
    ax.annotate("Residents under five in the five boroughs, thousands (Census Bureau population estimates; administrative, no margins)", xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left")
    fig.tight_layout(); _save(fig, path); plt.close(fig)


def fig_births(path: Path):
    plt = _mpl(); ys = sorted(BIRTHS)
    fig, ax = plt.subplots(figsize=(FIG_W, 3.0), dpi=100); fig.patch.set_alpha(0)
    ax.bar(np.arange(len(ys)), [BIRTHS[y] / 1000 for y in ys], color=ORANGE, width=0.8)
    ax.set_xticks(np.arange(len(ys))); ax.set_xticklabels([str(y) for y in ys], fontsize=8, rotation=45)
    ax.set_ylim(0, 140); ax.set_yticks([0, 50, 100]); ax.set_yticklabels(["0", "50k", "100k"])
    ax.grid(axis="y", color=GRID, linewidth=0.8); ax.grid(axis="x", visible=False)
    for s in ("top", "right"): ax.spines[s].set_visible(False)
    ax.set_title("Live births, 2010 to 2023", loc="left", fontsize=12, color="#1a1a1a", pad=22)
    ax.annotate("Births to New York City residents, thousands (Health Department Summary of Vital Statistics; administrative)", xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left")
    fig.tight_layout(); _save(fig, path); plt.close(fig)


# ------------------------------------------------------------------- companion data


def companion_dataset(frames: list[pd.DataFrame], rng: np.random.Generator) -> pd.DataFrame:
    parts = []
    for h in frames:
        k = h[h["child"] | (h["parent"].fillna(False)) | (h["workforce"] != "")].copy()
        parts.append(pd.DataFrame({
            "year": k["year"], "weight": k["PWGTP"].astype(int), "record": np.where(k["child"], "child", np.where(k["workforce"] != "", "childcare_worker", "parent")),
            "status": k["status"], "borough": k["borough"], "destination": k["destination"].replace("", np.nan), "age": k["AGEP"], "sex": _num(k["SEX"]).map({1: "male", 2: "female"}),
            "race_ethnicity": k["race"], "foreign_born": k["foreign_born"], "own_child": k["own_child"].where(k["child"]), "parents_arrangement": k["arrangement"].replace("", np.nan),
            "enrolled_in_school_3_4": k["enrolled"].where(k["age34"]), "public_school_3_4": k["public_school"].where(k["age34"]),
            "parent": k["parent"], "youngest_child_age": k["youngest"], "in_labor_force": k["in_lf"].where(k["adult16"]), "employed": k["employed"].where(k["adult16"]), "works_35_hours_plus": k["fulltime"].where(k["adult16"]),
            "gave_birth_past_year": k["new_mother"].where(_num(k["FER"]).notna()), "childcare_workforce_group": k["workforce"].replace("", np.nan), "hourly_earnings_2024": k["hourly"].round(2),
            "household_income_2024": k["hh_income"].round(0), "household_income_excluding_own_2024": k["other_income"].where(k["parent"]).round(0), "income_band": k["income_band"].replace("na", np.nan), "owner": k["owner"].where(k["household"]), "renter": k["renter"].where(k["household"]),
            "below_poverty": k["poverty"].where(k["poverty_known"]), "rent_burdened": k["rent_burden"].where(k["rent_ratio_known"]),
        }))
    df = pd.concat(parts, ignore_index=True)
    for c in df.columns:
        if c in ("year", "weight", "age", "youngest_child_age", "hourly_earnings_2024", "household_income_2024"): continue
        if df[c].dropna().isin([True, False]).all() and df[c].notna().any():
            df[c] = df[c].map({True: 1, False: 0}).astype("Int64")
    df.insert(0, "row_id", rng.permutation(len(df)) + 1)
    return df.sort_values("row_id").reset_index(drop=True)


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


def main(argv=None) -> int:
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--seed", type=int, default=20240101)
    ap.add_argument("--years", nargs="*", type=int, default=YEARS)
    ap.add_argument("--manifest-only", action="store_true", help="rewrite build.json for the files already in place (run after article_parents_dashboard.py and article_precision.py)")
    args = ap.parse_args(argv)
    if args.manifest_only:
        write_build_manifest(); print(f"build manifest -> {SITE_OUT / 'build.json'}"); return 0
    rng = np.random.default_rng(args.seed); OUT.mkdir(parents=True, exist_ok=True)

    print("published tables ...")
    pub = published_tables()
    for y in ("2019", "2024"):
        p = pub[y]; print(f"  {y}: own children under 6 {p['under6']['estimate']:,.0f}, both parents in LF {p['u6_both_lf']['estimate']:,.0f}, father only {p['u6_father_only_lf']['estimate']:,.0f}")

    years: dict[int, dict] = {}; pums_rows = []; frames = {}; censuses = {}; households = {}
    for y in args.years:
        try:
            h, rw, hh, hrw = load_year(y)
        except FileNotFoundError as e:
            print(f"  {y}: not on disk ({e}), skipped"); continue
        cs = Census(h, rw); hs = Census(hh, hrw); res, rows = year_estimates(cs, hs); years[y] = res; pums_rows += rows; censuses[y] = cs; households[y] = hs
        if y in FOCUS: frames[y] = h
        a = res["arrangement"]["under6"]; w = res["workforce"]
        print(f"  {y}: children {res['children']['total']['count']:,.0f} (under 6 {res['children']['under6']['count']:,.0f}); under-6 both parents in LF {a['both_lf']['share']:.1%}, a parent at home {a['parent_home']['share']:.1%}; "
              f"nannies {w['nanny']['count']:,.0f} (${w['nanny']['median_hourly']:.0f}/h), centre workers {w['centre']['count']:,.0f}; young-child out-rate {res['leaving']['by_band']['c1_5']['out_rate']:.2%}; "
              f"mothers of under-6s employed {res['parents']['youngest_under6']['mothers']['employed'][0]:.1%}")

    print("reproduction ...")
    repro = reproduction(pub, censuses)
    print(repro[["statistic", "published", "reproduced", "within_margin"]].to_string(index=False))
    print("population estimates, births, schools ...")
    u5 = under5_estimates(); sch = school_enrolment()
    for r in u5["years"][-3:]: print(f"  under 5, July {r['year']}: {r['total']:,}")
    for r in sch["years"][-2:]: print(f"  {r['year']}: 3-K {r['three_k']:,}, pre-K {r['pre_k']:,}, K {r['kindergarten_with_charters']:,} (with charters)")

    y15, y24 = years.get(2015), years.get(2024)
    change = {}
    if y15 and y24:
        change["under6_2015_2024"] = y24["children"]["under6"]["count"] - y15["children"]["under6"]["count"]
        change["under6_pct_2015_2024"] = y24["children"]["under6"]["count"] / y15["children"]["under6"]["count"] - 1
        change["nannies_2015_2024"] = y24["workforce"]["nanny"]["count"] - y15["workforce"]["nanny"]["count"]
        change["preschool_public_age3_points_2015_2024"] = (y24["preschool"]["age3"]["public_school"]["share"] - y15["preschool"]["age3"]["public_school"]["share"]) * 100
        change["hh_under6_2015_2024"] = y24["households"]["under6"]["count"] - y15["households"]["under6"]["count"]
        change["children_total_2015_2024"] = y24["children"]["total"]["count"] - y15["children"]["total"]["count"]
    res = {"article": SLUG, "as_of": "2026-09-24",
           "definitions": {"universe": "residents of the five boroughs in each one-year file; children under 18, their parents (householder, spouse or partner living with a child of the householder), employed childcare workers; leavers from the national file",
                           "dollar_year": 2024, "deflator": "New York-metro CPI-U (all items) to 2024 dollars after the Bureau's within-year factor", "census_margin": "90% (successive difference replication, 80 replicate weights)",
                           "arrangements": dict(ARRANGEMENTS), "age_bands": {k: l for k, l, *_ in AGE_BANDS}, "workforce": dict(WORKFORCE), "cost_cases": dict(COST_CASES), "cost_thresholds": {k: l for k, l, _ in COST_THRESHOLDS},
                           "rates_source": "dccny_market_rates_2024 / ocfs_market_rates_2024_attachment_a (weekly maximum rates, county group 4, effective 1 October 2024) times 52 weeks, as in nyc_comptroller_child_care_affordability_2025"},
           "published": pub, "years": {str(y): r for y, r in years.items()}, "change": change, "under5": u5, "births": {"note": "Live births to city residents, transcribed from Table PC1 of the Summary of Vital Statistics 2023 (nyc_summary_of_vital_statistics_2023).",
                                                                                                                       "years": {str(k): v for k, v in BIRTHS.items()}}, "schools": sch,
           "context": {"comptroller_spotlight_2025": {"family_care_infant_toddler_2024": 18_200, "centre_infant_toddler_2024": 26_000, "centre_toddler_2024": 23_400, "income_for_7pct_toddler_centre": 334_000,
                                                     "vouchers_infants_toddlers_sept_2024": 31_600, "income_eligible_infants_toddlers": 108_000, "voucher_share": 0.29, "children_0_3_change_2020_2023": -0.193,
                                                     "households_with_children_change_2020_2023": -0.087, "employed_child_care_workers": 32_917, "source": "nyc_comptroller_child_care_affordability_2025"},
                       "osc_report_19_2026": {"under5_population_2024_v2024": 454_220, "under5_change_2020_2024": -0.154, "served_under5_fy2025": 159_348, "coverage_rate_under5": 0.35, "coverage_available_3_4": 0.86, "coverage_available_under3": 0.13,
                                              "three_k_participants_fy2025": 44_729, "three_k_share_of_3yo": 0.53, "three_k_unfilled": 7_988, "pre_k_participants_fy2025": 56_184, "pre_k_unfilled": 15_181,
                                              "three_year_olds_2020": 106_656, "three_year_olds_2024": 85_059, "vouchers_under5_fy2023": 29_615, "vouchers_under5_fy2025": 52_027, "vouchers_under5_fy2026_first_four_months": 64_129,
                                              "children_on_vouchers_close_fy2025": 110_000, "infant_toddler_participants": 6_378, "infant_toddler_unfilled": 4_105, "voucher_pause": "May 2025", "cost_per_child": 19_989,
                                              "family_child_care_monthly_enrolment_fy2025": 34_248, "source": "osc_child_care_under_five_2026"},
                       "two_k_fall_2026": {"seats": 2_000, "state_funding_first_year": 73_000_000, "communities": ["District 6 (Washington Heights, Inwood, Hamilton Heights)", "District 10 (Fordham, Belmont, Norwood)", "Districts 18 and 23 (Canarsie, Brownsville)", "District 27 (Ozone Park, Richmond Hill)"],
                                           "applications_opened": "2026-06-02", "offers": "2026-08-04", "eligible": "children born in 2024", "source": "nyc_mayor_2k_applications_open_2026"},
                       "ccdf_final_rule_2026": {"effective": "2026-07-13", "rescinds": "the March 2024 requirement that family co-payments be capped at 7% of family income", "source": "acf_ccdf_final_rule_2026"},
                       "nces_2019": {"under_1_in_weekly_nonparental_care": 0.42, "age_1_2": 0.55, "age_3_5": 0.74, "source": "nces_ecpp_2019"}},
           "reproduction": repro.to_dict(orient="records"), "seed": args.seed}

    print("interactive chart cubes ...")
    charts = CH.build(censuses, households, res["years"], u5, sch, BIRTHS)
    CH.write(OUT / "charts.json", charts)
    print(f"  {len(charts['charts'])} charts, {sum(len(c['cells']) for c in charts['charts'].values())} filter combinations, {(OUT / 'charts.json').stat().st_size / 1024:.0f} KB")
    print("figures ...")
    for n, cid in enumerate(["arrangement", "where", "households", "parents", "preschool", "cost", "workforce", "leaving", "destinations"], start=1):
        fig_from_cube(charts["charts"][cid], OUT / f"fig{n}_{cid}.svg")
    fig_under5(u5, OUT / "fig10_under5.svg"); fig_births(OUT / "fig11_births.svg")
    fig_from_cube(charts["charts"]["schools"], OUT / "fig12_schools.svg", height=3.4)

    companion_dataset([frames[y] for y in FOCUS if y in frames], rng).to_csv(OUT / "data.csv", index=False)
    pd.DataFrame(pums_rows).to_csv(OUT / "pums_direct.csv", index=False)
    repro.to_csv(OUT / "reproduction.csv", index=False)
    (OUT / "results.json").write_text(json.dumps(res, indent=1, default=float), encoding="utf-8")
    SITE_OUT.mkdir(parents=True, exist_ok=True)
    for f in OUT.glob("*"):
        if f.suffix in (".svg", ".csv", ".json") and f.name != "build.json": shutil.copy2(f, SITE_OUT / f.name)
    shutil.copy2(Path(__file__), SITE_OUT / "analysis.py")
    write_build_manifest()
    print(f"outputs -> {OUT.relative_to(P.root)} and {SITE_OUT.relative_to(P.root)}")
    return 0


def write_build_manifest() -> None:
    """One manifest per build: the article version this run was made for (read from
    the article's frontmatter), the checksum of this script and of every file served
    next to it. tests/test_article_parents.py refuses a build whose manifest does not
    match the served files or the article's version, so a reader who downloads the
    script can check it is the one that produced the numbers on the page."""
    import hashlib, re
    mdx = P.root / "site" / "src" / "content" / "articles" / f"{SLUG}.mdx"
    m = re.search(r'^version:\s*"([^"]+)"', mdx.read_text(encoding="utf-8"), re.M)
    files = sorted(f for f in SITE_OUT.glob("*") if f.is_file() and f.name != "build.json")
    manifest = {"article": SLUG, "article_version": m.group(1) if m else None, "seed": 20240101,
                "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
                "files": {f.name: hashlib.sha256(f.read_bytes()).hexdigest() for f in files},
                "note": "Every file listed was written by one run of analysis.py (the script's checksum is script_sha256). "
                        "The article's version box shows article_version; if the two disagree the page and the files are from different builds."}
    text = json.dumps(manifest, indent=1)
    (OUT / "build.json").write_text(text, encoding="utf-8"); (SITE_OUT / "build.json").write_text(text, encoding="utf-8")


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