"""What people do in New York: the jobs in the city, what its residents do for a living,
and how both have changed since 1990 — every number and figure in the article, with
the working.

    py src/article_jobs.py
    py src/article_jobs.py --manifest-only     # rewrite build.json after the dashboard and precision scripts

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
======================================================================================

"What do people do in New York?" has two honest answers that are usually run
together: the jobs located in the five boroughs, whoever holds them, and the work the
city's residents do, wherever it is. The first is counted by employers, the second by
the census. The article keeps them apart and reads each to its latest date.

* **Payroll jobs in the city, 1990 to August 2026.** The Bureau of Labor Statistics'
  Current Employment Statistics for New York city (data/reference/ces_nyc_monthly.csv,
  built by fetch_ces_nyc.py from the public API): jobs on employers' payrolls located
  in the five boroughs, by supersector, every month since January 1990. A sample
  survey of employers, benchmarked each year to unemployment-insurance records, with
  the latest months preliminary; BLS publishes its sampling error separately and no
  margin is drawn here. Not seasonally adjusted, so annual averages and same-month
  comparisons are used. Government here includes public schools and hospitals.

* **The published census tables.** ACS 2024 one-year tables C24030 (industry) and
  C24010 (occupation) for the civilian employed population of New York city, with
  B23025 for the totals (data/reference/tables/acs_c24030_c24010_nyc_us_2024_1yr.json).
  Every row is reproduced from the microdata before the article goes anywhere the
  tables do not.

* **The census microdata, six vintages.** Employed residents of the five boroughs from
  the 1990 and 2000 decennial 5% samples and the 2010 ACS (data/interim/acs1/
  commute_<year>.csv.gz, built by extract_commute_history.py, which keeps each year's
  industry, occupation and class-of-worker codes), and from the 2015-2024 one-year
  ACS files for the metro study area (metro_person_<year>.csv.gz), with replicate
  weights where the file has them (2010 on). The 2015-2024 files also hold every
  worker in the study area whose workplace is in the five boroughs, which is the
  census's own count of jobs in the city and of the commuters who hold them.

* **The 2026 view of residents.** The Current Population Survey basic monthly files,
  June 2024 to August 2026 (data/interim/cps_telework_nyc.csv.gz, built by
  extract_cps_telework.py): employed-at-work residents with the survey's major
  industry and occupation recodes, class of worker and hours. About 650 to 750
  records a month, pooled by period; adjacent months share most of their sample.

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

* **Employed residents**: civilian residents of the five boroughs aged 16 and over
  who were employed in the survey week (employment status 1 or 2). Armed forces
  excluded, as in the published tables.
* **Industry groups** (fourteen, harmonised across three coding schemes). 2000 and
  later use the NAICS-based census code, by its first two characters: agriculture and
  mining (11, 21); construction (23); manufacturing (31-33); wholesale trade (42);
  retail trade (44-45); transportation, warehousing and utilities (48-49, 22);
  information (51); finance, insurance and real estate (52-53); professional,
  scientific, management and administrative services (54-56); educational services
  (61); health care and social assistance (62); arts, entertainment, recreation,
  accommodation and food (71-72); other services (81); public administration (92).
  1990 uses the 1990 census industry codes mapped to the same groups; the mapping is
  the code's own, with the known seams named in the methods FAQ (publishing sat in
  manufacturing in 1990 and in information from 2000; social services were their own
  1990 heading and are health and social assistance from 2000). The CPS major recode
  combines wholesale with retail and education with health, so the 2026 view shows
  those pairs together.
* **Occupation groups** (eight): management, business and finance; professional
  (computing, engineering, science, community service, legal, education, arts and
  media, health practitioners and technicians); service (health support, protective,
  food, cleaning, personal care); sales; office and administrative support;
  construction, extraction, farming and installation and repair; production;
  transportation and material moving. 1990 codes are mapped from the 1990 scheme
  (technicians go with professionals, handlers and labourers with transportation).
* **Class of worker**: private (for-profit and non-profit), government (federal,
  state, local), self-employed (incorporated or not).
* **Place of work**: the workplace state and place-of-work area; the five boroughs are
  identified by their area codes in each vintage (1990 codes 5000-5400, 2000/2010
  codes 3700-4100, 2020 codes 4100-4500), Manhattan by its own code. Works from home
  is the usual means of transportation. The census records a home worker's place of
  work as their home address, not their employer's, so every workplace measure here
  (works in Manhattan, works in the five boroughs, jobs located in the city and the
  commuter shares) is computed among workers who do not usually work from home, and
  home workers are reported separately.
* **Earnings**: the person's own wage, salary and self-employment earnings in the
  income year, in 2024 dollars. The 2015-2024 extracts already carry the Bureau's
  within-year factor and the New York-metro CPI (extract_metro_pums.py applies both),
  so their dollars are used as they stand; the 1990, 2000 and 2010 extracts are
  nominal and are deflated here with the annual New York-metro CPI of the income year:
  1989 for the 1990 census, 1999 for the 2000 census, 2010 for the 2010 file, after the
  2010 file's within-year factor (ADJINC). Medians are among workers with positive
  earnings; earnings in the old censuses heap on round amounts, so several 1990 group
  medians share a value.
* **Degree**: bachelor's or higher, workers 25 and over; **foreign-born**: born outside
  the United States and its territories (1990: naturalised or not a citizen).

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

PUMS files from 2010 carry 80 replicate weights: 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 at 90%. The 1990 and 2000 5% samples
carry no replicate weights, so no margin is calculated for them (they are samples,
of 125,000 and 132,000 employed city residents, not censuses of every worker).
Medians from the replicate-weight files carry a Woodruff-style approximate 90%
interval (median_moe): the replicate variation of the share below the median, read
back as dollars on the group's earnings distribution. CES estimates carry sampling error that BLS estimates
separately; none is drawn here. CPS shares carry a tick of 1.645 times the standard
error of the monthly estimates within a period: a measure of month-to-month
dispersion, not a formal survey margin, since adjacent months share three quarters of
their sample and real monthly change is mixed in.

Outputs (all regenerable; nothing is edited by hand):
    output/articles/what_people_do/results.json      every number quoted in the prose
    output/articles/what_people_do/pums_direct.csv   census-side estimates with margins
    output/articles/what_people_do/reproduction.csv  the published-table reproduction
    output/articles/what_people_do/charts.json       cubes behind the interactive charts
    output/articles/what_people_do/fig*.svg          static twins of the same figures
    output/articles/what_people_do/data.csv          companion dataset, 2010 and 2024 records (re-randomised ids)
    output/articles/what_people_do/data_1990_2000.csv.gz  the 1990 and 2000 5% sample records, gzipped (the Worker serves assets up to 25 MiB)
    output/articles/what_people_do/build.json        version and checksums of the served files
    site/public/articles/what-people-do/             the same, copied for the website
"""

from __future__ import annotations

import argparse
import hashlib
import json
import re
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, NYC_MIGPUMA, Z90, MIN_RECORDS, _geo, vintage, cpi_deflators, race_label
from article_parents import fig_from_cube
import article_jobs_charts as CH

SLUG = "what-people-do"
OUT = OUTPUT / "articles" / "what_people_do"
SITE_OUT = P.root / "site" / "public" / "articles" / SLUG
CENSUS_YEARS = [1990, 2000, 2010, 2015, 2019, 2024]
FOCUS = (1990, 2000, 2010, 2019, 2024)
TABLES_2024 = REFERENCE / "tables" / "acs_c24030_c24010_nyc_us_2024_1yr.json"
CES_CSV = REFERENCE / "ces_nyc_monthly.csv"
CPS_CSV = INTERIM / "cps_telework_nyc.csv.gz"

INDUSTRIES = [("agri_mining", "Agriculture and mining"), ("construction", "Construction"), ("manufacturing", "Manufacturing"), ("wholesale", "Wholesale trade"),
              ("retail", "Retail trade"), ("transport_utilities", "Transportation, warehousing and utilities"), ("information", "Information"),
              ("financial", "Finance, insurance and real estate"), ("professional_business", "Professional, management and administrative services"),
              ("education", "Educational services"), ("health", "Health care and social assistance"), ("leisure_hospitality", "Arts, entertainment, hotels and restaurants"),
              ("other_services", "Other services"), ("public_admin", "Public administration")]
OCCUPATIONS = [("management_business", "Management, business and finance"), ("professional", "Professional (computing, engineering, science, law, teaching, arts, health)"),
               ("service", "Service (health support, protective, food, cleaning, personal care)"), ("sales", "Sales"), ("office_admin", "Office and administrative support"),
               ("construction_repair", "Construction, extraction, farming, installation and repair"), ("production", "Production"), ("transportation", "Transportation and material moving")]
CLASSES = [("private", "Private employer"), ("government", "Government"), ("self_employed", "Self-employed")]
AGE_BANDS = [("16_24", "16 to 24", 16, 25), ("25_34", "25 to 34", 25, 35), ("35_44", "35 to 44", 35, 45), ("45_54", "45 to 54", 45, 55), ("55_64", "55 to 64", 55, 65), ("65_plus", "65 and over", 65, 200)]
EARN_BANDS = [("lt25k", "Under $25k", -np.inf, 25_000), ("25k_50k", "$25k–50k", 25_000, 50_000), ("50k_75k", "$50k–75k", 50_000, 75_000), ("75k_100k", "$75k–100k", 75_000, 100_000),
              ("100k_150k", "$100k–150k", 100_000, 150_000), ("150k_plus", "$150k and over", 150_000, np.inf)]
TRAITS = [("female", "Women"), ("foreign_born", "Born abroad"), ("ba_plus", "Bachelor's degree or higher (25+)"), ("self_employed", "Self-employed"), ("government", "Government employee"),
          ("works_manhattan", "Works in Manhattan"), ("works_nyc", "Works in the five boroughs"), ("wfh", "Works from home"), ("fulltime", "Works 35 hours or more"), ("age_55_plus", "Aged 55 and over")]
# CES supersectors: series key -> (label, census industry groups that correspond)
CES_SECTORS = [("construction", "Construction (with mining and logging)", ["construction", "agri_mining"]), ("manufacturing", "Manufacturing", ["manufacturing"]),
               ("wholesale", "Wholesale trade", ["wholesale"]), ("retail", "Retail trade", ["retail"]), ("transport_utilities", "Transportation, warehousing and utilities", ["transport_utilities"]),
               ("information", "Information", ["information"]), ("financial", "Financial activities", ["financial"]), ("professional_business", "Professional and business services", ["professional_business"]),
               ("education", "Educational services (private)", ["education"]), ("health", "Health care and social assistance (private)", ["health"]),
               ("leisure_hospitality", "Leisure and hospitality", ["leisure_hospitality"]), ("other_services", "Other services", ["other_services"]),
               ("government", "Government (all levels, including public schools and hospitals)", ["public_admin"])]
CPS_INDUSTRY = {1: "agri_mining", 2: "agri_mining", 3: "construction", 4: "manufacturing", 5: "trade", 6: "transport_utilities", 7: "information", 8: "financial", 9: "professional_business",
                10: "education_health", 11: "leisure_hospitality", 12: "other_services", 13: "public_admin"}
CPS_INDUSTRIES = [("agri_mining", "Agriculture and mining"), ("construction", "Construction"), ("manufacturing", "Manufacturing"), ("trade", "Wholesale and retail trade"),
                  ("transport_utilities", "Transportation and utilities"), ("information", "Information"), ("financial", "Financial activities"), ("professional_business", "Professional and business services"),
                  ("education_health", "Education and health services"), ("leisure_hospitality", "Leisure and hospitality"), ("other_services", "Other services"), ("public_admin", "Public administration")]
CPS_OCCUPATION = {1: "management_business", 2: "professional", 3: "service", 4: "sales", 5: "office_admin", 6: "construction_repair", 7: "construction_repair", 8: "construction_repair", 9: "production", 10: "transportation"}
CPS_CLASS = {1: "government", 2: "government", 3: "government", 4: "private", 5: "self_employed", 6: "self_employed"}   # PRCOW1: 1-3 government, 4 private, 5-6 self-employed, 7 without pay
CPS_PERIODS = [("2024_jun_dec", "June–December 2024", lambda y, m: y == 2024), ("2025", "January–December 2025 (no October)", lambda y, m: y == 2025), ("2026_jan_aug", "January–August 2026", lambda y, m: y == 2026)]
POW_NYC = {"1990": {5000: "Bronx", 5100: "Manhattan", 5200: "Staten Island", 5300: "Brooklyn", 5400: "Queens"}, "2000": {3700: "Bronx", 3800: "Manhattan", 3900: "Staten Island", 4000: "Brooklyn", 4100: "Queens"},
           "2010": {3700: "Bronx", 3800: "Manhattan", 3900: "Staten Island", 4000: "Brooklyn", 4100: "Queens"}, "2020": NYC_MIGPUMA["2020"]}


# ------------------------------------------------------------------ code mapping


def industry_from_naics(naics) -> np.ndarray:
    s = pd.Series(naics).fillna("").astype(str).str.strip().str.upper()
    p2 = s.str[:2]
    return np.select([s == "", p2.isin(["11", "21"]), p2 == "23", p2.isin(["31", "32", "33", "3M"]), p2 == "42", p2.isin(["44", "45", "4M"]), p2.isin(["48", "49", "22"]), p2 == "51",
                      p2.isin(["52", "53"]), p2.isin(["54", "55", "56"]), p2 == "61", p2 == "62", p2.isin(["71", "72"]), p2 == "81", s.str.startswith("928110"), p2 == "92"],
                     ["", "agri_mining", "construction", "manufacturing", "wholesale", "retail", "transport_utilities", "information", "financial", "professional_business", "education", "health",
                      "leisure_hospitality", "other_services", "military", "public_admin"], "")


def industry_from_1990(code, variant: str = "main") -> np.ndarray:
    """1990 census industry codes to the fourteen groups. `variant="alt"` is the sensitivity
    mapping for the two seams the harmonisation cannot close: newspapers and other publishing
    (171, 172) go to information rather than manufacturing, as the 2000 scheme moved publishing,
    and the 1990 social-service headings (861 to 871) stay with other services rather than
    joining health care as NAICS sector 62 folds them."""
    c = pd.to_numeric(pd.Series(code), errors="coerce").to_numpy(dtype=float)
    def rng(a, b): return (c >= a) & (c <= b)
    if variant == "alt":
        out = industry_from_1990(code, "main").copy()
        out[rng(171, 172)] = "information"; out[rng(861, 871)] = "other_services"
        return out
    # 742 (automotive rental and leasing) and 801 (video rental) are rental and leasing, which
    # NAICS places in sector 53 with real estate; 770 (lodging places other than hotels) is
    # accommodation, beside 762 (hotels). All three sit inside the 1990 business, repair and
    # personal-service ranges that otherwise map to other services and to arts and recreation.
    conds = [np.isnan(c) | (c == 0), rng(10, 50), c == 60, rng(100, 392), rng(400, 432) | rng(450, 472), rng(440, 442), rng(500, 571), c == 641, rng(580, 691),
             rng(700, 712) | (c == 742) | (c == 801),
             rng(721, 741) | (c == 841) | rng(882, 893), rng(750, 761) | rng(771, 791) | (c == 873) | rng(880, 881), (c == 762) | (c == 770), (c == 800) | rng(802, 810) | (c == 872),
             rng(812, 840) | rng(861, 871), rng(842, 860), rng(900, 932), rng(940, 960)]
    vals = ["", "agri_mining", "construction", "manufacturing", "transport_utilities", "information", "wholesale", "leisure_hospitality", "retail", "financial",
            "professional_business", "other_services", "leisure_hospitality", "leisure_hospitality", "health", "education", "public_admin", "military"]
    return np.select(conds, vals, "")


def occupation_from_code(code, scheme: str, variant: str = "main") -> np.ndarray:
    """Occupation codes to the eight groups. For the 1990 scheme, `variant="alt"` follows the
    1990 census's own grouping where the harmonisation departs from it: technicians (203 to
    235) sit with office and administrative support, as the 1990 major group "technical, sales
    and administrative support" had them, rather than with professionals; and handlers, helpers
    and labourers (863 to 889) sit with production, as the 1990 "operators, fabricators and
    labourers" group had them, rather than with transportation and material moving."""
    c = pd.to_numeric(pd.Series(code), errors="coerce").to_numpy(dtype=float)
    if scheme == "1990" and variant == "alt":
        out = occupation_from_code(code, "1990", "main").copy()
        out[(c >= 203) & (c <= 235)] = "office_admin"; out[(c >= 863) & (c <= 889)] = "production"
        return out
    if scheme == "2000": c = c * 10          # the 2000 census stores the 2000 codes as three digits
    def rng(a, b): return (c >= a) & (c <= b)
    if scheme == "1990":
        conds = [np.isnan(c) | (c == 0), rng(3, 37), rng(43, 235), rng(243, 285), rng(303, 389), rng(403, 469), rng(473, 617), rng(628, 799), rng(803, 889), rng(903, 905)]
    else:
        conds = [np.isnan(c) | (c == 0), rng(10, 950), rng(1000, 3599), rng(4700, 4999), rng(5000, 5999), rng(3600, 4699), rng(6000, 7699), rng(7700, 8999), rng(9000, 9799), c >= 9800]
    vals = ["", "management_business", "professional", "sales", "office_admin", "service", "construction_repair", "production", "transportation", "military"]
    return np.select(conds, vals, "")


def class_from_cow(cow) -> np.ndarray:
    c = pd.to_numeric(pd.Series(cow), errors="coerce").to_numpy(dtype=float)
    return np.select([np.isin(c, [1, 2]), np.isin(c, [3, 4, 5]), np.isin(c, [6, 7])], ["private", "government", "self_employed"], "")


def earn_band(v) -> np.ndarray:
    v = np.asarray(v, float)
    out = np.select([v < 25_000, v < 50_000, v < 75_000, v < 100_000, v < 150_000], ["lt25k", "25k_50k", "50k_75k", "75k_100k", "100k_150k"], "150k_plus")
    return np.where(np.isnan(v), "na", out)


def age_band(a) -> np.ndarray:
    v = np.asarray(a, float)
    return np.select([v < 25, v < 35, v < 45, v < 55, v < 65], ["16_24", "25_34", "35_44", "45_54", "55_64"], "65_plus")


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


def _finish(h: pd.DataFrame, year: int) -> pd.DataFrame:
    h["employed"] = h["esr"].isin([1, 2]) & (h["age"] >= 16)
    h["female"] = h["sex"] == 2
    h["adult25"] = h["age"] >= 25
    h["ba_plus"] = h["ba_plus"].astype(bool) & h["adult25"]
    h["foreign_born"] = h["foreign_born"].astype(bool)
    h["self_employed"] = h["cls"] == "self_employed"; h["government"] = h["cls"] == "government"; h["private"] = h["cls"] == "private"
    h["wfh"] = h["mode"] == 11
    # The journey-to-work questions are asked of workers at work in the survey week; those
    # with a job but absent that week have no answer, so working from home is measured
    # among workers at work, the universe of the Bureau's own commuting tables.
    h["at_work"] = h["mode"].notna()
    h["fulltime"] = h["hours"] >= 35
    # Full-time, full-year: usual hours of 35 or more and 50 to 52 weeks worked, both in the
    # earnings reference year, so the comparison holds work intensity roughly constant.
    h["ftfy"] = h["fulltime"] & (h["fullyear"] == 1)
    h["age_55_plus"] = h["age"] >= 55
    h["age_band"] = age_band(h["age"]); h["earn_band"] = earn_band(h["earnings_2024"])
    h["year"] = year
    return h


def load_census_year(year: int) -> tuple[pd.DataFrame, np.ndarray | None]:
    """Employed residents of the five boroughs in one census year (and, for 2015 on, every
    worker in the study area whose workplace is in the city, flagged `resident` False)."""
    defl = cpi_deflators()
    if year <= 2010:
        d = pd.read_csv(INTERIM / "acs1" / f"commute_{year}.csv.gz", dtype={"ind": str, "naics": str, "occ": str, "puma": str})
        rw = np.load(INTERIM / "acs1" / f"commute_{year}_repwts.npy", mmap_mode="r") if (INTERIM / "acs1" / f"commute_{year}_repwts.npy").exists() else None
        v = "1990" if year == 1990 else "2000" if year == 2000 else "2010"
        # The 1990 5% file codes SEX 0 male / 1 female; 2000 and the ACS code 1 male / 2 female.
        h = pd.DataFrame({"borough": d["borough"], "weight": d["weight"].astype(float), "age": d["age"], "sex": d["sex"] + (1 if year == 1990 else 0), "esr": d["esr"], "mode": d["mode"], "hours": d["hours"],
                          "race": d["race"], "foreign_born": d["foreign_born"], "ba_plus": d["ba_plus"], "pow_state": d["pow_state"], "pow_puma": d["pow_puma"]})
        # Usual weekly hours and weeks worked in the earnings year. For 1990 the extract's `hours`
        # is hours worked last week, so the usual-hours field (HOUR89) replaces it; for 2000 and
        # 2010 the two are the same field.
        h["hours"] = d["hours_usual"]; h["fullyear"] = d["fullyear"].astype(float)
        h["industry"] = industry_from_1990(d["ind"]) if year == 1990 else industry_from_naics(d["naics"])
        if year == 1990: h["ind_code"] = d["ind"].to_numpy(); h["occ_code"] = d["occ"].to_numpy()
        h["occupation"] = occupation_from_code(d["occ"], "1990" if year == 1990 else "2000" if year == 2000 else "acs")
        h["cls"] = class_from_cow(d["cow"])
        income_year = {1990: 1989, 2000: 1999, 2010: 2010}[year]      # the 1990 census asked about 1989 income, the 2000 census about 1999
        # The 2010 ACS file carries the Bureau's within-year factor (ADJINC), which the 2015-2024
        # extracts already have applied; the decennial samples asked about a calendar year and need none.
        within = d["adjinc"].fillna(1.0) if "adjinc" in d else 1.0
        h["earnings_2024"] = np.where(d["earnings"] > 0, d["earnings"] * within * defl.get(income_year, np.nan), np.nan)
        h["resident"] = True
        codes = POW_NYC[v]
    else:
        v = vintage(year); geo = _geo(v); codes = POW_NYC[v]
        cols = ["SERIALNO", "PUMA", "STATE", "PWGTP", "AGEP", "SEX", "RAC1P", "HISP", "NATIVITY", "SCHL", "ESR", "COW", "PERNP", "ADJINC", "WKHP", "JWTRNS", "OCCP", "NAICSP", "POWSP", "POWPUMA"]
        header = pd.read_csv(INTERIM / "acs1" / f"metro_person_{year}.csv.gz", nrows=0).columns
        wk = "WKWN" if "WKWN" in header else "WKW"      # weeks worked: a count from 2019, bands before (1 = 50-52)
        cols = cols + [wk]
        p = pd.read_csv(INTERIM / "acs1" / f"metro_person_{year}.csv.gz", usecols=cols, dtype={"STATE": str, "PUMA": str, "SERIALNO": str, "OCCP": str, "NAICSP": str, "POWPUMA": str})
        rw_all = np.load(INTERIM / "acs1" / f"metro_person_{year}_repwts.npy", mmap_mode="r")
        p["puma_geoid"] = p["STATE"].str.zfill(2) + p["PUMA"].str.zfill(5)
        p = p.merge(geo, on="puma_geoid", how="left")
        pow_state = pd.to_numeric(p["POWSP"], errors="coerce"); pow_puma = pd.to_numeric(p["POWPUMA"], errors="coerce")
        in_city_work = (pow_state == 36) & pow_puma.isin(list(codes))
        keep = ((p["in_nyc"] == 1) | in_city_work).to_numpy()
        p = p[keep].reset_index(drop=True); rw = np.asarray(rw_all[keep], dtype=np.float32)
        # extract_metro_pums.py has already applied the Bureau's within-year factor and the
        # New York-metro CPI to every dollar column of the 2015-2024 files, so PERNP here
        # is in 2024 dollars as it stands; nothing is applied again.
        pernp = pd.to_numeric(p["PERNP"], errors="coerce")
        h = pd.DataFrame({"borough": np.where(p["in_nyc"] == 1, p["borough"], ""), "weight": p["PWGTP"].astype(float), "age": p["AGEP"], "sex": pd.to_numeric(p["SEX"], errors="coerce"),
                          "esr": pd.to_numeric(p["ESR"], errors="coerce"), "mode": pd.to_numeric(p["JWTRNS"], errors="coerce"), "hours": pd.to_numeric(p["WKHP"], errors="coerce"),
                          "race": race_label(p["RAC1P"], p["HISP"]), "foreign_born": (pd.to_numeric(p["NATIVITY"], errors="coerce") == 2).astype(int),
                          "ba_plus": (pd.to_numeric(p["SCHL"], errors="coerce") >= 21).astype(int), "pow_state": pow_state[keep].to_numpy(), "pow_puma": pow_puma[keep].to_numpy()})
        h["industry"] = industry_from_naics(p["NAICSP"]); h["occupation"] = occupation_from_code(p["OCCP"], "acs"); h["cls"] = class_from_cow(p["COW"])
        h["earnings_2024"] = np.where(pernp > 0, pernp, np.nan)
        h["resident"] = (p["in_nyc"] == 1).to_numpy()
        weeks = pd.to_numeric(p[wk], errors="coerce")
        h["fullyear"] = ((weeks >= 50) if wk == "WKWN" else (weeks == 1)).astype(float).to_numpy()
    h["works_nyc"] = (h["pow_state"] == 36) & h["pow_puma"].isin(list(codes))
    man = [k for k, b in codes.items() if b == "Manhattan"][0]
    h["works_manhattan"] = (h["pow_state"] == 36) & (h["pow_puma"] == man)
    # A home worker's place of work is their home address in the census, so the
    # workplace measures are computed among workers who do not usually work from home.
    h["pow_known"] = h["pow_state"].notna() & (h["mode"] != 11)
    h["pow_borough"] = np.where(h["works_nyc"], h["pow_puma"].map(codes), "")
    h = _finish(h, year)
    return h, (np.asarray(rw, dtype=np.float32) if rw is not None else None)


class Plain:
    """Weighted estimates without replicate weights (the 1990 and 2000 5% samples)."""

    def __init__(self, h: pd.DataFrame):
        self.h, self.w = h, h["weight"].to_numpy(dtype=float)

    def count(self, mask): return float(self.w[np.asarray(mask)].sum()), float("nan")
    def share(self, mask, num):
        mask = np.asarray(mask); w = self.w[mask]
        return (float((w * np.asarray(num)[mask]).sum() / w.sum()) if w.sum() else float("nan")), float("nan")
    def diff(self, a, b): return self.count(a)[0] - self.count(b)[0], float("nan")
    def median(self, mask, col):
        mask = np.asarray(mask); v = self.h.loc[mask, col].to_numpy(dtype=float); ok = ~np.isnan(v)
        return weighted_median(v[ok], self.w[mask][ok]) if ok.sum() >= MIN_RECORDS else float("nan")


def median_moe(est, mask, col: str) -> tuple[float, float]:
    """Weighted median and, where the file has replicate weights, an approximate 90% margin
    in the Woodruff style (not the Bureau's design-factor recipe): the replicate-weight standard error
    of the share of the group below the median, turned into a dollar interval by reading
    50% minus and plus 1.645 standard errors off the group's (interpolated) earnings
    distribution. Recomputing the median itself on each replicate does not work here:
    earnings are reported in round amounts, so the median sits on the same figure in every
    replicate and the margin would come out as zero."""
    mask = np.asarray(mask); v = est.h.loc[mask, col].to_numpy(dtype=float); ok = ~np.isnan(v)
    if ok.sum() < MIN_RECORDS: return float("nan"), float("nan")
    x, w = v[ok], est.w[mask][ok]
    full = weighted_median(x, w)
    if isinstance(est, Plain): return float(full), float("nan")
    rw = np.asarray(est.rw[mask][ok], dtype=float)
    below = x < full
    p_full = w[below].sum() / w.sum(); p_reps = rw[below].sum(0) / rw.sum(0)
    se = np.sqrt(4 / 80 * ((p_reps - p_full) ** 2).sum())
    # The distribution is read on $2,500 intervals with linear interpolation inside each, as the
    # Bureau's table medians are: earnings heap on round amounts, and an unbinned distribution
    # would put the whole interval inside one heap.
    edges = np.arange(0, max(x.max(), 1) + 2500, 2500.0)
    cdf = np.concatenate([[0.0], np.cumsum(np.histogram(x, bins=edges, weights=w)[0]) / w.sum()])
    q = lambda p: float(np.interp(p, cdf, edges))
    lo, hi = q(max(0.5 - Z90 * se, 0.0)), q(min(0.5 + Z90 * se, 1.0))
    return float(full), float((hi - lo) / 2)


def estimator(h: pd.DataFrame, rw):
    if rw is None: return Plain(h)
    h = h.copy(); h["PWGTP"] = h["weight"]
    return Census(h, rw)


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


def published_tables() -> dict:
    j = json.loads(TABLES_2024.read_text(encoding="utf-8"))
    def cell(geo, table, male, female):
        t = j["data"][geo][table]
        est = t["estimate"][f"{table}{male:03d}"] + t["estimate"][f"{table}{female:03d}"]
        moe = float(np.sqrt(t["error"][f"{table}{male:03d}"] ** 2 + t["error"][f"{table}{female:03d}"] ** 2))
        return {"estimate": est, "moe90": moe}
    # C24030 male columns 3..28, female = male + 27; C24010 male 3..37, female = male + 35.
    ind = {"agri_mining": 3, "construction": 6, "manufacturing": 7, "wholesale": 8, "retail": 9, "transport_utilities": 10, "information": 13, "financial": 14, "professional_business": 17,
           "education": 22, "health": 23, "leisure_hospitality": 24, "other_services": 27, "public_admin": 28}
    # C24010 columns: 3 management/business/science/arts (4 management/business/financial, 7 computer/engineering/science, 11 education/legal/community/arts/media,
    # 16 health practitioners), 19 service, 27 sales and office (28 sales, 29 office), 30 natural resources/construction/maintenance, 34 production/transportation
    # (35 production, 36 transportation, 37 material moving); female columns are the male ones plus 36.
    occ = {"management_business": [4], "professional": [7, 11, 16], "service": [19], "sales": [28], "office_admin": [29], "construction_repair": [30], "production": [35], "transportation": [36, 37]}
    out = {"table_note": "C24030 (industry) and C24010 (occupation), civilian employed population 16 and over; sexes summed, margins combined in quadrature.",
           "total": {"estimate": j["data"]["16000US3651000"]["C24030"]["estimate"]["C24030001"], "moe90": j["data"]["16000US3651000"]["C24030"]["error"]["C24030001"]},
           "industry": {k: cell("16000US3651000", "C24030", c, c + 27) for k, c in ind.items()}, "occupation": {}}
    t = j["data"]["16000US3651000"]["C24010"]
    for k, cols in occ.items():
        est = sum(t["estimate"][f"C24010{c:03d}"] + t["estimate"][f"C24010{c + 36:03d}"] for c in cols)
        moe = float(np.sqrt(sum(t["error"][f"C24010{c:03d}"] ** 2 + t["error"][f"C24010{c + 36:03d}"] ** 2 for c in cols)))
        out["occupation"][k] = {"estimate": est, "moe90": moe}
    b = j["data"]["16000US3651000"]["B23025"]
    out["employment"] = {"population_16_plus": b["estimate"]["B23025001"], "labor_force": b["estimate"]["B23025002"], "employed": b["estimate"]["B23025004"], "unemployed": b["estimate"]["B23025005"], "not_in_lf": b["estimate"]["B23025007"]}
    bu = j["data"]["01000US"]["C24030"]
    out["us_industry_shares"] = {k: (bu["estimate"][f"C24030{c:03d}"] + bu["estimate"][f"C24030{c + 27:03d}"]) / bu["estimate"]["C24030001"] for k, c in ind.items()}
    return out


def reproduction(pub: dict, est, h: pd.DataFrame) -> pd.DataFrame:
    rows = []; res = (h["resident"] & h["employed"]).to_numpy(); ind = h["industry"].to_numpy(); occ = h["occupation"].to_numpy()
    def put(label, p, mask, key):
        v, vm = est.count(mask); band = float(np.sqrt(p["moe90"] ** 2 + vm ** 2))
        rows.append({"statistic": label, "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": 2024, "cell": key})
    put("Employed residents 16 and over, NYC, 2024", pub["total"], res, "total")
    for k, l in INDUSTRIES: put(f"{l}, NYC, 2024", pub["industry"][k], res & (ind == k), f"ind_{k}")
    for k, l in OCCUPATIONS: put(f"{l.split(' (')[0]} occupations, NYC, 2024", pub["occupation"][k], res & (occ == k), f"occ_{k}")
    return pd.DataFrame(rows)


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


def year_estimates(est, h: pd.DataFrame) -> tuple[dict, list[dict]]:
    year = int(h["year"].iloc[0]); rows = []
    res = (h["resident"] & h["employed"]).to_numpy(); ind = h["industry"].to_numpy(); occ = h["occupation"].to_numpy(); cls = h["cls"].to_numpy(); boro = h["borough"].to_numpy()
    traits = {t: h[t].to_numpy() for t, _ in TRAITS}
    bases = {"ba_plus": h["adult25"].to_numpy(), "works_manhattan": h["pow_known"].to_numpy(), "works_nyc": h["pow_known"].to_numpy(), "wfh": h["at_work"].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})

    def profile(mask) -> dict:
        prof = {"n_records": int(mask.sum())}
        for t, l in TRAITS:
            base = mask & bases.get(t, np.ones(len(h), bool)); v, vm = est.share(base, traits[t]); prof[t] = {"share": v, "moe90": vm}
        prof["median_earnings"] = est.median(mask, "earnings_2024"); prof["median_age"] = est.median(mask, "age")
        return prof

    tot, totm = est.count(res)
    out: dict = {"year": year, "n_records": int(res.sum()), "employed": {"count": tot, "moe90": totm}, "has_margins": not isinstance(est, Plain)}
    put("employed", "Employed residents 16 and over", tot, totm, res.sum(), "count")
    out["industry"] = {}
    for k, l in INDUSTRIES:
        mk = res & (ind == k); c, cm = est.count(mk); s, sm = est.share(res, ind == k)
        out["industry"][k] = {"label": l, "count": c, "moe90": cm, "share": s, "share_moe90": sm, **profile(mk)}
        put(f"industry_{k}", l, s, sm, res.sum(), "share")
    out["occupation"] = {}
    for k, l in OCCUPATIONS:
        mk = res & (occ == k); c, cm = est.count(mk); s, sm = est.share(res, occ == k)
        out["occupation"][k] = {"label": l, "count": c, "moe90": cm, "share": s, "share_moe90": sm, **profile(mk)}
        put(f"occupation_{k}", l, s, sm, res.sum(), "share")
    out["class"] = {k: {"label": l, "count": est.count(res & (cls == k))[0], "share": est.share(res, cls == k)[0], "share_moe90": est.share(res, cls == k)[1]} for k, l in CLASSES}
    out["all"] = profile(res)
    out["unclassified_industry_share"] = est.share(res, ind == "")[0]
    out["borough"] = {}
    for b in BOROUGHS:
        mb = res & (boro == b); c, cm = est.count(mb)
        top = sorted(((k, est.share(mb, ind == k)[0]) for k, _ in INDUSTRIES), key=lambda x: -x[1])[:3]
        out["borough"][b] = {"employed": c, "moe90": cm, "n_records": int(mb.sum()), "industry_shares": {k: est.share(mb, ind == k)[0] for k, _ in INDUSTRIES}, "top_industries": top,
                             "median_earnings": est.median(mb, "earnings_2024"), "works_manhattan": est.share(mb & bases["works_manhattan"], traits["works_manhattan"])[0],
                             "wfh": est.share(mb & bases["wfh"], traits["wfh"])[0], "foreign_born": est.share(mb, traits["foreign_born"])[0]}
    # Earnings by attachment to work: annual earnings mix pay with hours and weeks worked, so the
    # median is also taken among full-time and among full-time, full-year workers.
    earn = ~np.isnan(h["earnings_2024"].to_numpy()); ft = h["fulltime"].to_numpy(); ftfy = h["ftfy"].to_numpy()
    out["earnings_attachment"] = {
        "all": dict(zip(("median", "moe90"), median_moe(est, res & earn, "earnings_2024"))),
        "fulltime": dict(zip(("median", "moe90"), median_moe(est, res & earn & ft, "earnings_2024"))),
        "ftfy": dict(zip(("median", "moe90"), median_moe(est, res & earn & ftfy, "earnings_2024"))),
        "ftfy_share": est.share(res & earn, ftfy), "fulltime_share": est.share(res & earn, ft),
        "note": "median own earnings in 2024 dollars among employed residents with positive earnings; full-time is usual hours of 35 or more, full-year 50 to 52 weeks, both in the earnings year"}
    for k, _ in INDUSTRIES:
        out["industry"][k]["median_earnings_ftfy"] = est.median(res & (ind == k) & ftfy, "earnings_2024")
    for k, _ in OCCUPATIONS:
        out["occupation"][k]["median_earnings_ftfy"] = est.median(res & (occ == k) & ftfy, "earnings_2024")
    hl = res & (ind == "health")
    out["health_by_occupation"] = {k: {"share": est.share(hl, occ == k)[0], "median": est.median(hl & (occ == k), "earnings_2024"),
                                       "median_ftfy": est.median(hl & (occ == k) & ftfy, "earnings_2024"), "n_records": int((hl & (occ == k)).sum())} for k, _ in OCCUPATIONS}
    # Where residents work.
    known = res & bases["works_nyc"]
    out["place_of_work"] = {"universe": "workplace shares among employed residents who do not usually work from home and report a workplace; working from home among workers at work in the survey week", "works_nyc": est.share(known, traits["works_nyc"]), "works_manhattan": est.share(known, traits["works_manhattan"]),
                            "wfh": est.share(res & bases["wfh"], traits["wfh"]), "at_work": est.count(res & bases["wfh"]), "outside_city": est.share(known, ~traits["works_nyc"]), "n_records": int(known.sum())}
    # Jobs located in the city (2015 on): every worker in the study area with a city workplace.
    if (~h["resident"]).any():
        jobs = h["employed"].to_numpy() & traits["works_nyc"] & bases["works_nyc"]
        c, cm = est.count(jobs); rc, _ = est.count(jobs & h["resident"].to_numpy())
        out["jobs_in_city"] = {"total": c, "moe90": cm, "held_by_residents": rc, "commuter_share": est.share(jobs, ~h["resident"].to_numpy()), "n_records": int(jobs.sum()),
                               "note": "workers living in the study area (the five boroughs, Long Island, the Lower Hudson Valley, northern New Jersey, southwestern Connecticut) who do not usually work from home and whose workplace is in the five boroughs; commuters from beyond the study area are not in the files",
                               "resident_home_workers": est.count(res & traits["wfh"])[0],
                               "industry": {}, "by_pow_borough": {}}
        for k, l in INDUSTRIES:
            mk = jobs & (ind == k); jc, jm = est.count(mk)
            out["jobs_in_city"]["industry"][k] = {"label": l, "jobs": jc, "moe90": jm, "commuter_share": est.share(mk, ~h["resident"].to_numpy()), "median_earnings": est.median(mk, "earnings_2024"),
                                                  "median_earnings_commuters": est.median(mk & ~h["resident"].to_numpy(), "earnings_2024"), "median_earnings_residents": est.median(mk & h["resident"].to_numpy(), "earnings_2024"), "n_records": int(mk.sum())}
        pb = h["pow_borough"].to_numpy()
        for b in BOROUGHS:
            mk = jobs & (pb == b); jc, jm = est.count(mk)
            out["jobs_in_city"]["by_pow_borough"][b] = {"jobs": jc, "moe90": jm, "commuter_share": est.share(mk, ~h["resident"].to_numpy())[0], "n_records": int(mk.sum())}
        out["jobs_in_city"]["commuters_median_earnings"] = est.median(jobs & ~h["resident"].to_numpy(), "earnings_2024")
        out["jobs_in_city"]["residents_median_earnings"] = est.median(jobs & h["resident"].to_numpy(), "earnings_2024")
    return out, rows


# ------------------------------------------------------------------ administrative and CPS series


def ces_series() -> dict:
    t = pd.read_csv(CES_CSV)
    out = {"note": "Jobs on employers' payrolls located in the five boroughs, thousands, not seasonally adjusted, from the BLS Current Employment Statistics survey of employers (benchmarked annually to "
                   "unemployment-insurance records; the latest months are preliminary). Annual figures are means of the published months (12, or 11 for a year missing one); the current year is the "
                   "January-August mean, flagged. Government includes public schools and hospitals.", "annual": {}, "latest": {}, "monthly_total": [], "recovery": {}}
    last = t[t["series"] == "total_nonfarm"].sort_values(["year", "month"]).iloc[-1]; ly, lm = int(last["year"]), int(last["month"])
    out["latest_month"] = f"{ly}-{lm:02d}"; out["latest_is_preliminary"] = True
    rec = t[t["series"] == "total_nonfarm"].sort_values("thousands", ascending=False).iloc[0]
    out["record"] = {"month": f"{int(rec['year'])}-{int(rec['month']):02d}", "thousands": float(rec["thousands"]), "note": "the highest month in the series since January 1990, not seasonally adjusted"}
    for key, g in t.groupby("series"):
        ann = g.groupby("year")["thousands"].agg(["mean", "size"]).reset_index()
        out["annual"][key] = {int(r["year"]): {"thousands": float(r["mean"]), "months": int(r["size"]), "partial": bool(int(r["year"]) == ly and int(r["size"]) < 12)} for _, r in ann.iterrows()}
        gl = g[(g["year"] == ly) & (g["month"] == lm)]
        out["latest"][key] = float(gl["thousands"].iloc[0]) if len(gl) else None
    tot = t[t["series"] == "total_nonfarm"].sort_values(["year", "month"])
    out["monthly_total"] = [{"label": f"{int(r.year)}-{int(r.month):02d}", "thousands": float(r.thousands)} for r in tot[tot["year"] >= 2019].itertuples()]
    def at(key, y, m):
        g = t[(t["series"] == key) & (t["year"] == y) & (t["month"] == m)]; return float(g["thousands"].iloc[0]) if len(g) else None
    out["recovery"] = {"feb_2020": at("total_nonfarm", 2020, 2), "apr_2020": at("total_nonfarm", 2020, 4), "aug_2019": at("total_nonfarm", 2019, 8), "aug_2024": at("total_nonfarm", 2024, 8),
                       "aug_2025": at("total_nonfarm", 2025, 8), "latest": at("total_nonfarm", ly, lm), "latest_by_sector": {k: (at(k, ly, lm), at(k, 2019, lm), at(k, 2000, lm), at(k, 1990, lm)) for k, *_ in CES_SECTORS}}
    return out


def cps_series() -> dict | None:
    if not CPS_CSV.exists(): return None
    c = pd.read_csv(CPS_CSV)
    c["industry"] = c["ind_major"].map(CPS_INDUSTRY).fillna(""); c["occupation"] = c["occ_major"].map(CPS_OCCUPATION).fillna(""); c["cls"] = c["class_of_worker"].map(CPS_CLASS).fillna("")
    out = {"note": "Employed-at-work residents of the five boroughs from the CPS basic monthly files, pooled by period; the moe90 fields are 1.645 times the standard error of the monthly estimates, "
                   "a measure of month-to-month dispersion rather than a formal survey margin (adjacent months share about three quarters of their sample, and real monthly change is mixed in). No October 2025 file exists.", "periods": {}}
    for key, label, f in CPS_PERIODS:
        p = c[[f(y, m) for y, m in zip(c["year"], c["month"])]]
        if p.empty: continue
        months = sorted(set(zip(p["year"], p["month"])))
        def monthly_share(col, val):
            vals = []
            for y, m in months:
                q = p[(p["year"] == y) & (p["month"] == m)]; w = q["weight"]; vals.append(float((w * (q[col] == val)).sum() / w.sum()))
            return vals
        blk = {"label": label, "months": len(months), "n": int(len(p)), "employed_at_work": float(p["weight"].sum() / len(months)), "industry": {}, "occupation": {}, "class": {}}
        for col, keys, dest in (("industry", [k for k, _ in CPS_INDUSTRIES], "industry"), ("occupation", [k for k, _ in OCCUPATIONS], "occupation"), ("cls", [k for k, _ in CLASSES], "class")):
            for k in keys:
                v = np.array(monthly_share(col, k)); w = p["weight"]
                blk[dest][k] = {"share": float((w * (p[col] == k)).sum() / w.sum()), "moe90": float(Z90 * v.std(ddof=1) / np.sqrt(len(v))) if len(v) > 1 else float("nan"), "n_records": int((p[col] == k).sum())}
        blk["telework_share"] = float((p["weight"] * p["telework"]).sum() / p["weight"].sum())
        blk["full_time_share"] = float((p["weight"] * p["full_time"]).sum() / p["weight"].sum())
        blk["ba_share"] = float((p["weight"] * p["ba_plus"]).sum() / p["weight"].sum())
        out["periods"][key] = blk
    out["first_month"] = f"{int(c['year'].min())}-{int(c[c['year'] == c['year'].min()]['month'].min()):02d}"; out["last_month"] = f"{int(c['year'].max())}-{int(c[c['year'] == c['year'].max()]['month'].max()):02d}"
    return out


# ------------------------------------------------------------------ figures


def fig_payroll_monthly(ces: dict, path: Path):
    import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt
    from article_migration import INK, GRID, BLUE, FIG_W, _save
    plt.rcParams.update({"font.family": "sans-serif", "font.size": 10, "axes.edgecolor": GRID, "axes.labelcolor": INK, "xtick.color": INK, "ytick.color": INK, "text.color": INK, "svg.fonttype": "none"})
    ms = ces["monthly_total"]; fig, ax = plt.subplots(figsize=(FIG_W, 3.2), dpi=100); fig.patch.set_alpha(0)
    x = np.arange(len(ms)); ax.plot(x, [m["thousands"] / 1000 for m in ms], color=BLUE, linewidth=2)
    ticks = [i for i, m in enumerate(ms) if m["label"].endswith("-01")]
    ax.set_xticks(ticks); ax.set_xticklabels([ms[i]["label"][:4] for i in ticks], fontsize=8.5)
    ax.set_ylim(3.6, 5.0); ax.set_yticks([3.75, 4.0, 4.25, 4.5, 4.75, 5.0]); ax.set_yticklabels(["3.75m", "4.0m", "4.25m", "4.5m", "4.75m", "5.0m"])
    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(f"Payroll jobs in the city, month by month, to {CH.month_label(ces['latest_month'])}", loc="left", fontsize=12, color="#1a1a1a", pad=22)
    ax.annotate("Jobs on payrolls located in the five boroughs, total nonfarm, not seasonally adjusted (BLS Current Employment 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["employed"] & (h["resident"] | h["works_nyc"])].copy()
        parts.append(pd.DataFrame({
            "year": k["year"], "weight": k["weight"].round().astype(int), "resident_of_nyc": k["resident"], "borough": k["borough"].replace("", np.nan), "works_in_nyc": k["works_nyc"].where(k["pow_known"]),
            "works_in_manhattan": k["works_manhattan"].where(k["pow_known"]), "workplace_borough": k["pow_borough"].replace("", np.nan), "works_from_home": k["wfh"].where(k["at_work"]), "at_work_in_survey_week": k["at_work"],
            "industry_group": k["industry"].replace("", np.nan), "occupation_group": k["occupation"].replace("", np.nan), "class_of_worker": k["cls"].replace("", np.nan),
            "age": k["age"], "sex": k["sex"].map({1: "male", 2: "female"}), "race_ethnicity": k["race"], "foreign_born": k["foreign_born"], "bachelors_or_higher_25plus": k["ba_plus"].where(k["adult25"]),
            "usual_hours": k["hours"], "worked_50_52_weeks": k["fullyear"], "full_time_full_year": k["ftfy"], "earnings_2024_dollars": k["earnings_2024"].round(0),
        }))
    df = pd.concat(parts, ignore_index=True)
    for c in df.columns:
        if c in ("year", "weight", "age", "usual_hours", "earnings_2024_dollars"): 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)


def harmonisation_sensitivity(est, h: pd.DataFrame, y2024: dict) -> dict:
    """How far the 1990 shares, and the 1990-to-2024 changes, move when the ambiguous 1990
    codes are mapped the other reasonable way (see industry_from_1990 and occupation_from_code,
    variant "alt"). Reported for every group; the article quotes the largest movements."""
    res = (h["resident"] & h["employed"]).to_numpy()
    out = {"note": "1990 shares of employed residents under the article's mapping (main) and the alternative (alt), and the 1990-to-2024 change in points under each; "
                   "alt moves newspapers and other publishing to information, keeps the 1990 social-service headings in other services, puts technicians with office work and handlers and labourers with production",
           "industry": {}, "occupation": {}}
    for kind, main, alt, groups in (("industry", h["industry"].to_numpy(), industry_from_1990(h["ind_code"], "alt"), INDUSTRIES),
                                    ("occupation", h["occupation"].to_numpy(), occupation_from_code(h["occ_code"], "1990", "alt"), OCCUPATIONS)):
        for k, _ in groups:
            sm = est.share(res, main == k)[0] * 100; sa = est.share(res, alt == k)[0] * 100; s24 = y2024[kind][k]["share"] * 100
            out[kind][k] = {"main_1990": sm, "alt_1990": sa, "change_main": s24 - sm, "change_alt": s24 - sa, "shift": sa - sm}
    return out


def deflator_sensitivity(e1990: dict, e2024: dict) -> dict:
    """The 1990-to-2024 change in median earnings under the New York-metro CPI (the article's
    deflator) and under the Bureau of Labor Statistics' research series R-CPI-U-RS, which the
    Census Bureau uses for its own constant-dollar comparisons. The 1990 medians are taken
    back to 1989 dollars with the metro index and forward with the research series."""
    import openpyxl
    wb = openpyxl.load_workbook(RAW / "cpi" / "r-cpi-u-rs-allitems.xlsx", read_only=True, data_only=True)
    rs = {}
    for row in wb.active.iter_rows(values_only=True):
        if row and isinstance(row[0], (int, float)) and 1977 <= row[0] <= 2030 and isinstance(row[-1], (int, float)):
            rs[int(row[0])] = float(row[-1])      # last column is the annual average
    metro = cpi_deflators()
    out = {"note": "median earnings, 1990 (income year 1989) and 2024, in 2024 dollars under two deflators", "r_cpi_u_rs": {"1989": rs.get(1989), "2024": rs.get(2024)}, "groups": {}}
    for g in ("all", "fulltime", "ftfy"):
        m1990_metro = e1990[g]["median"]; nominal_1989 = m1990_metro / metro[1989]
        m1990_rs = nominal_1989 * rs[2024] / rs[1989]
        out["groups"][g] = {"nominal_1989": nominal_1989, "metro_cpi_1990": m1990_metro, "r_cpi_u_rs_1990": m1990_rs, "median_2024": e2024[g]["median"],
                            "change_pct_metro": (e2024[g]["median"] / m1990_metro - 1) * 100, "change_pct_r_cpi_u_rs": (e2024[g]["median"] / m1990_rs - 1) * 100}
    return out


def write_build_manifest() -> None:
    mdx = P.root / "site" / "src" / "content" / "articles" / f"{SLUG}.mdx"
    m = re.search(r'^version:\s*"([^"]+)"', mdx.read_text(encoding="utf-8"), re.M) if mdx.exists() else None
    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")


# ------------------------------------------------------------------------- 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=CENSUS_YEARS)
    ap.add_argument("--manifest-only", action="store_true")
    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()
    print(f"  2024: employed {pub['total']['estimate']:,.0f}; health {pub['industry']['health']['estimate']:,.0f}, finance/insurance/real estate {pub['industry']['financial']['estimate']:,.0f}")
    years: dict[int, dict] = {}; pums_rows = []; frames = {}; ests = {}
    for y in args.years:
        try:
            h, rw = load_census_year(y)
        except FileNotFoundError as e:
            print(f"  {y}: not on disk ({e}), skipped"); continue
        est = estimator(h, rw); res, rows = year_estimates(est, h); years[y] = res; pums_rows += rows; frames[y] = h; ests[y] = (est, h)
        top = sorted(res["industry"].items(), key=lambda kv: -kv[1]["share"])[:3]
        print(f"  {y}: employed residents {res['employed']['count']:,.0f} (n {res['n_records']:,}); top industries " + ", ".join(f"{k} {v['share']:.1%}" for k, v in top)
              + f"; unclassified {res['unclassified_industry_share']:.2%}; median earnings ${res['all']['median_earnings']:,.0f}; WFH {res['place_of_work']['wfh'][0]:.1%}"
              + (f"; jobs in city {res['jobs_in_city']['total']:,.0f}, commuter share {res['jobs_in_city']['commuter_share'][0]:.1%}" if "jobs_in_city" in res else ""))
        ea = res["earnings_attachment"]; print(f"        median earnings all ${ea['all']['median']:,.0f}, full-time ${ea['fulltime']['median']:,.0f}, full-time full-year ${ea['ftfy']['median']:,.0f} (share {ea['ftfy_share'][0]:.1%})")
    print("reproduction ...")
    repro = reproduction(pub, *ests[2024]) if 2024 in ests else pd.DataFrame()
    if len(repro): print(repro[["statistic", "published", "reproduced", "within_margin"]].to_string(index=False))
    print("payroll and CPS ...")
    ces = ces_series(); cps = cps_series()
    r = ces["recovery"]; print(f"  CES total nonfarm: Feb 2020 {r['feb_2020']:,.1f}k, Apr 2020 {r['apr_2020']:,.1f}k, latest {ces['latest_month']} {r['latest']:,.1f}k (preliminary), record {ces['record']['month']} {ces['record']['thousands']:,.1f}k")
    if cps:
        for k, b in cps["periods"].items(): print(f"  CPS {b['label']}: n {b['n']:,}, education & health {b['industry']['education_health']['share']:.1%}, professional {b['occupation']['professional']['share']:.1%}")

    res = {"article": SLUG, "as_of": "2026-09-26",
           "definitions": {"universe": "civilian employed residents of the five boroughs aged 16 and over; jobs in the city are workers in the study area with a workplace in the five boroughs",
                           "dollar_year": 2024, "deflator": "New York-metro CPI-U (all items) to 2024 dollars after the Bureau's within-year factor where the file has one; 1989 income (1990 census) deflated with the 1990 index",
                           "census_margin": "90% (successive difference replication, 80 replicate weights) from 2010; none for the 1990 and 2000 5% samples", "industries": dict(INDUSTRIES), "occupations": dict(OCCUPATIONS),
                           "classes": dict(CLASSES), "cps_industries": dict(CPS_INDUSTRIES), "ces_sectors": {k: l for k, l, _ in CES_SECTORS}},
           "published": pub, "years": {str(y): r for y, r in years.items()}, "ces": ces, "cps": cps, "reproduction": repro.to_dict(orient="records"), "seed": args.seed}
    if 1990 in years and 2024 in years:
        res["deflator_sensitivity"] = deflator_sensitivity(years[1990]["earnings_attachment"], years[2024]["earnings_attachment"])
        res["harmonisation_sensitivity"] = harmonisation_sensitivity(ests[1990][0], ests[1990][1], years[2024])
        res["change"] = {"employed_1990_2024": years[2024]["employed"]["count"] - years[1990]["employed"]["count"],
                         "industry_points_1990_2024": {k: (years[2024]["industry"][k]["share"] - years[1990]["industry"][k]["share"]) * 100 for k, _ in INDUSTRIES},
                         "occupation_points_1990_2024": {k: (years[2024]["occupation"][k]["share"] - years[1990]["occupation"][k]["share"]) * 100 for k, _ in OCCUPATIONS},
                         "median_earnings_pct_1990_2024": {g: (years[2024]["earnings_attachment"][g]["median"] / years[1990]["earnings_attachment"][g]["median"] - 1) * 100
                                                           for g in ("all", "fulltime", "ftfy")}}

    print("interactive chart cubes ...")
    charts = CH.build(ests, res["years"], ces, cps)
    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(["payroll", "residents_industry", "residents_occupation", "who", "earnings", "jobs_in_city", "cps", "borough"], start=1):
        fig_from_cube(charts["charts"][cid], OUT / f"fig{n}_{cid}.svg")
    fig_payroll_monthly(ces, OUT / "fig9_payroll_monthly.svg")

    companion_dataset([frames[y] for y in (2010, 2024) if y in frames], rng).to_csv(OUT / "data.csv", index=False)
    companion_dataset([frames[y] for y in (1990, 2000) if y in frames], rng).to_csv(OUT / "data_1990_2000.csv.gz", index=False, compression="gzip")
    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", ".gz") and f.name != "build.json": shutil.copy2(f, SITE_OUT / f.name)
    shutil.copy2(Path(__file__), SITE_OUT / "analysis.py")
    # The extraction script is published too, so the 1990/2000/2010 definitions (degree codes,
    # hours, weeks, place of work) can be checked without the raw census files.
    shutil.copy2(Path(__file__).with_name("extract_commute_history.py"), SITE_OUT / "extract_history.py")
    write_build_manifest()
    print(f"outputs -> {OUT.relative_to(P.root)} and {SITE_OUT.relative_to(P.root)}")
    return 0


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