"""Who still goes to the office — every number and figure in the article, with the working.

    py src/article_office.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
======================================================================================

Six years after the offices emptied, who in New York City works from home, who still
travels to a job, what did that do to the commute and to Manhattan, and where does it
stand in 2026?

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

* **The published journey-to-work table.** Census table B08301 (workers 16 and over by
  means of transportation to work, whose last row is "worked from home") for New York
  city, its five counties and the United States, from the one-year ACS: 2019 from the
  Bureau's sequence-format summary file (data/raw/acs_sf/), 2021-2023 from the
  table-based summary files (data/reference/tables/acsdt1y*-b08301.dat), 2024 from the
  Census Reporter mirror (data/reference/tables/acs_commute_nyc_us_2024_1yr.json and
  acs_b08301_nyc_counties_2024_1yr.json). No 2020 one-year table was published. These
  are the numbers a reader could find on their own; the article reproduces each of
  them before going anywhere the table does not.

* **The census microdata.** The 2015-2024 one-year PUMS files for the study area
  (data/interim/acs1/metro_person_<year>.csv.gz with their 80 replicate weights),
  restricted to residents of the five boroughs who were at work in the survey week.
  Every count, share and margin labelled "census" is a weighted tabulation of those
  records. 2019 and 2024 carry the article; the other years give the curve. The same
  files hold the suburban residents of the metro study area, which is how the
  inbound Manhattan commuters are counted. The 2020 file carries the Bureau's
  experimental weights and is shown with that label.

* **The Current Population Survey, June 2024 to August 2026.** Since June 2024 the CPS
  asks every employed person whether they teleworked or worked at home for pay at any
  time last week, and for how many hours. The five NYC counties are identified in the
  public file (data/interim/cps_telework_nyc.csv.gz, built by extract_cps_telework.py).
  It is a small sample and has no replicate weights; the article reports it with an
  approximate margin and labels it as the 2026 view. It also asks a different
  question from the census: any hours at home, not the usual place of work.

* **MTA ridership.** Daily subway, bus, railroad and bridge-and-tunnel counts through
  the first days of September 2026 (data/raw/transit/), administrative data that
  counts trips rather than commuters, used for the shape of the office week in 2026.

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

* **Workers**: residents of the five boroughs, 16 and over, employed and at work in the
  survey reference week (civilian or armed forces), so with a means of transportation
  recorded. This is the universe of table B08301.
* **Worked from home**: the census's means-of-transportation answer "worked from home",
  which the questionnaire asks as the way the person *usually* got to work last week.
  A hybrid worker who was in the office three days and at home two is a commuter by
  that rule. Everything the census says about working from home here is about the
  usual place of work, and the CPS section says what the rule leaves out.
* **Commuters**: workers who did not usually work from home.
* **Earnings**: the worker's own earnings, converted to 2024 dollars with the New
  York-metro CPI-U (data/reference/cpi_annual.csv) after the Bureau's own within-year
  adjustment, in fixed bands.
* **Office-based industries**: the worker's industry is information, finance and
  insurance, real estate, professional and technical services, management of
  companies, or administrative and support services (NAICS sectors 51-56). This is the
  conventional "office-using" grouping; it is an industry test, not a test of whether
  the job is at a desk.
* **Occupation groups**: the census occupation codes collapsed into the Bureau's major
  groups. The 2015-2017 files use the 2010 occupation codes; the major-group ranges are
  the same to within a handful of codes.
* **Place of work Manhattan**: the place-of-work area for Manhattan (3800 in the 2010
  vintage used through 2021, 4100 in the 2020 vintage from 2022). A worker "physically
  commutes to Manhattan" if their place of work is Manhattan and they did not usually
  work from home.
* **Race/ethnicity**: Hispanic of any race; otherwise White, Black, Asian, Other.
* **Degree**: a bachelor's degree or higher.

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

PUMS ships 80 replicate weights per person. 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. Published margins are the Bureau's own; a share derived from two
published cells is compared to our estimate using our margin.

The CPS has no replicate weights in its public file. Its margins here are approximate:
each month is treated as a separate sample, the pooled estimate is the weighted mean
of the monthly estimates, and its standard error is the spread of the monthly
estimates divided by the square root of the number of months. Because the CPS
re-interviews the same households in consecutive months, this understates the true
error somewhat; the article says "about" and rounds accordingly.

Outputs (all regenerable; nothing is edited by hand):
    output/articles/who_still_goes_to_the_office/results.json    every number quoted in the prose
    output/articles/who_still_goes_to_the_office/pums_direct.csv  the census-side estimates with margins
    output/articles/who_still_goes_to_the_office/reproduction.csv the published-table reproduction
    output/articles/who_still_goes_to_the_office/charts.json      cubes behind the interactive charts
    output/articles/who_still_goes_to_the_office/fig*.svg         static twins of the same figures
    output/articles/who_still_goes_to_the_office/data.csv         companion dataset (re-randomised ids)
    site/public/articles/who-still-goes-to-the-office/            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
import article_office_charts as CH
from article_commute import mta_series, weighted_median

SLUG = "who-still-goes-to-the-office"
OUT = OUTPUT / "articles" / "who_still_goes_to_the_office"
SITE_OUT = P.root / "site" / "public" / "articles" / SLUG
Z90 = 1.645
MIN_RECORDS = 30
YEARS = list(range(2015, 2025))
FOCUS = (2019, 2024)

BOROUGHS = ["Manhattan", "Brooklyn", "Queens", "Bronx", "Staten Island"]
REGIONS = [("nyc", "New York City"), ("long_island", "Long Island"), ("lower_hudson", "Lower Hudson Valley"),
           ("north_jersey", "Northern New Jersey"), ("southwest_ct", "Southwestern Connecticut")]
RESIDENCES = [(b, b if b != "Bronx" else "The Bronx") for b in BOROUGHS] + [r for r in REGIONS if r[0] != "nyc"]
MAN_POW = {"2010": 3800, "2020": 4100}
MODE_LABEL = {1: "car", 2: "bus", 3: "subway", 4: "commuter rail", 5: "light rail", 6: "ferry", 7: "taxi",
              8: "motorcycle", 9: "bicycle", 10: "walked", 11: "worked from home", 12: "other"}


def mode_group(code: pd.Series) -> np.ndarray:
    c = code.to_numpy(dtype=float)
    return np.select([c == 11, c == 3, c == 2, np.isin(c, [1, 8]), c == 4, c == 10, c == 9],
                     ["home", "subway", "bus", "car", "commuter rail", "walked", "bicycle"], "other")


MODE_GROUPS = [("home", "Worked from home"), ("subway", "Subway"), ("bus", "Bus"), ("car", "Car, truck or van"),
               ("commuter rail", "Commuter rail"), ("walked", "Walked"), ("bicycle", "Bicycle"), ("other", "Taxi, ferry and other")]
EARN_BANDS = [("lt30k", "Under $30k", -np.inf, 30_000), ("30k_60k", "$30k–60k", 30_000, 60_000),
              ("60k_100k", "$60k–100k", 60_000, 100_000), ("100k_150k", "$100k–150k", 100_000, 150_000),
              ("150k_plus", "$150k and over", 150_000, np.inf)]
RACES = ["White", "Black", "Hispanic", "Asian", "Other"]
AGE_BANDS = [("16_24", "16–24", 16, 25), ("25_34", "25–34", 25, 35), ("35_44", "35–44", 35, 45), ("45_54", "45–54", 45, 55),
             ("55_64", "55–64", 55, 65), ("65_plus", "65 and over", 65, 200)]
SECTORS = [("office", "Office-based industries"), ("health", "Health care and social assistance"), ("education", "Education"),
           ("retail_trade", "Retail and wholesale"), ("food_accommodation", "Restaurants and hotels"), ("construction", "Construction"),
           ("transport", "Transportation and warehousing"), ("manufacturing", "Manufacturing"), ("government", "Public administration"),
           ("arts_other", "Arts, personal services and other")]
OCC_GROUPS = [(10, 440, "Management"), (500, 750, "Business and finance"), (800, 960, "Business and finance"),
              (1005, 1240, "Computer and math"), (1305, 1560, "Architecture and engineering"), (1600, 1980, "Science"),
              (2001, 2060, "Community and social service"), (2100, 2180, "Legal"), (2205, 2555, "Education"),
              (2600, 2920, "Arts, media and entertainment"), (3000, 3550, "Health practitioners"),
              (3601, 3655, "Health support"), (3700, 3960, "Protective service"), (4000, 4160, "Food service"),
              (4200, 4255, "Cleaning and maintenance"), (4330, 4655, "Personal care and service"), (4700, 4965, "Sales"),
              (5000, 5940, "Office and administrative"), (6005, 6130, "Farming"), (6200, 6765, "Construction"),
              (6800, 6950, "Extraction"), (7000, 7640, "Installation and repair"), (7700, 8990, "Production"),
              (9005, 9760, "Transportation and moving"), (9800, 9830, "Military")]
CLASSES = [("private", "Private employer"), ("government", "Government"), ("self_employed", "Self-employed")]
INK = "#6f6e69"; GRID = "#d5d4cd"; BLUE = "#2a78d6"; ORANGE = "#eb6834"; TICK = "#1a1a1a"; FIG_W = 6.4

TABLES_2024_1YR = REFERENCE / "tables" / "acs_commute_nyc_us_2024_1yr.json"
TABLES_2024_COUNTIES = REFERENCE / "tables" / "acs_b08301_nyc_counties_2024_1yr.json"
SF_2019 = RAW / "acs_sf" / "NewYork_All_Geographies_2019_1yr"
CPS_FILE = INTERIM / "cps_telework_nyc.csv.gz"
CPS_SUMMARY = INTERIM / "cps_telework_summary.json"
GEOS = {"01000US": "United States", "16000US3651000": "New York City", "05000US36005": "Bronx", "05000US36047": "Brooklyn",
        "05000US36061": "Manhattan", "05000US36081": "Queens", "05000US36085": "Staten Island"}


# ------------------------------------------------------------------ derived variables


def occ_group(s: pd.Series) -> pd.Series:
    c = pd.to_numeric(s, errors="coerce"); out = pd.Series("Other", index=s.index)
    for lo, hi, lab in OCC_GROUPS:
        out[(c >= lo) & (c <= hi)] = lab
    return out


def sector(naics: pd.Series) -> np.ndarray:
    s = naics.fillna("").astype(str).str[:2]
    return np.select([s.isin(["51", "52", "53", "54", "55", "56"]), s == "62", s == "61", s.isin(["44", "45", "42", "4M"]),
                      s == "72", s == "23", s.isin(["48", "49"]), s.isin(["31", "32", "33", "3M"]), s == "92"],
                     ["office", "health", "education", "retail_trade", "food_accommodation", "construction", "transport", "manufacturing", "government"],
                     "arts_other")


def race_label(rac1p, hisp) -> np.ndarray:
    r = pd.to_numeric(rac1p, errors="coerce").to_numpy(); h = pd.to_numeric(hisp, errors="coerce").to_numpy()
    return np.where(h > 1, "Hispanic", np.where(r == 1, "White", np.where(r == 2, "Black", np.where(r == 6, "Asian", "Other"))))


def earn_band(e: pd.Series) -> np.ndarray:
    v = e.to_numpy(dtype=float)
    return np.select([v < 30_000, v < 60_000, v < 100_000, v < 150_000], ["lt30k", "30k_60k", "60k_100k", "100k_150k"], "150k_plus")


def age_band(a: pd.Series) -> np.ndarray:
    v = a.to_numpy(dtype=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")


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


def cpi_deflators() -> dict[int, float]:
    d = pd.read_csv(REFERENCE / "cpi_annual.csv")
    d = d[d["name"] == "cpi_ny_all_items"]
    return dict(zip(d["year"].astype(int), d["deflator_to_2024"].astype(float)))


# ------------------------------------------------------------------ census side


def _geo(vintage: str) -> pd.DataFrame:
    f = "geography_pumas_2010.csv" if vintage == "2010" else "geography_pumas.csv"
    g = pd.read_csv(REFERENCE / f, dtype={"state_fips": str, "puma": str})
    g["puma_geoid"] = g["state_fips"].str.zfill(2) + g["puma"].str.zfill(5)
    g["in_nyc"] = pd.to_numeric(g["in_nyc"]).astype(int)
    return g[["puma_geoid", "in_nyc", "borough", "region"]]


def load_year(year: int) -> tuple[pd.DataFrame, np.ndarray]:
    """Every worker at work in the one-year study-area file (NYC and suburbs), with
    replicate weights and the article's derived variables."""
    cols = ["PUMA", "STATE", "PWGTP", "AGEP", "SEX", "ESR", "JWTRNS", "JWMNP", "POWSP", "POWPUMA", "PERNP", "ADJINC",
            "WKHP", "RAC1P", "HISP", "SCHL", "OCCP", "NAICSP", "NATIVITY", "COW"]
    p = pd.read_csv(INTERIM / "acs1" / f"metro_person_{year}.csv.gz", usecols=cols, dtype={"STATE": str, "PUMA": str, "OCCP": str, "NAICSP": str})
    rw = np.load(INTERIM / "acs1" / f"metro_person_{year}_repwts.npy", mmap_mode="r")
    assert len(rw) == len(p)
    vintage = "2010" if year <= 2021 else "2020"
    p["puma_geoid"] = p["STATE"].str.zfill(2) + p["PUMA"].str.zfill(5)
    p = p.merge(_geo(vintage), on="puma_geoid", how="left")
    keep = (p["in_nyc"].notna() & p["JWTRNS"].notna()).to_numpy()
    h = p[keep].reset_index(drop=True); rw = np.asarray(rw[keep])
    h["year"] = year
    h["nyc"] = h["in_nyc"] == 1
    h["mode"] = mode_group(h["JWTRNS"])
    h["wfh"] = h["mode"] == "home"
    h["commuter"] = ~h["wfh"]
    defl = cpi_deflators().get(year, 1.0)
    h["earnings"] = h["PERNP"].fillna(0) * h["ADJINC"] / 1e6 * defl
    h["earn_band"] = earn_band(h["earnings"])
    h["race"] = race_label(h["RAC1P"], h["HISP"])
    h["ba_plus"] = h["SCHL"] >= 21
    h["female"] = h["SEX"] == 2
    h["foreign_born"] = h["NATIVITY"] == 2
    h["age_band"] = age_band(h["AGEP"])
    h["occ"] = occ_group(h["OCCP"])
    h["sector"] = sector(h["NAICSP"])
    h["office_industry"] = h["sector"] == "office"
    h["cow"] = class_of_worker(h["COW"])
    h["full_time"] = h["WKHP"] >= 35
    pow_state = pd.to_numeric(h["POWSP"], errors="coerce"); pow_puma = pd.to_numeric(h["POWPUMA"], errors="coerce")
    h["pow_manhattan"] = (pow_state == 36) & (pow_puma == MAN_POW[vintage])
    h["to_manhattan"] = h["pow_manhattan"] & h["commuter"]
    h["residence"] = np.where(h["nyc"], h["borough"], h["region"])
    return h, rw


class Census:
    """Weighted estimates with successive-difference-replication margins (90%)."""

    def __init__(self, h: pd.DataFrame, rw: np.ndarray):
        self.h, self.rw, self.w = h, rw, h["PWGTP"].to_numpy(dtype=float)

    def share(self, mask, num) -> tuple[float, float]:
        mask = np.asarray(mask); num = np.asarray(num)
        w, rw = self.w[mask], self.rw[mask]; x = num[mask].astype(float)
        if w.sum() == 0: return float("nan"), float("nan")
        full = (w * x).sum() / w.sum(); reps = (rw * x[:, None]).sum(0) / rw.sum(0)
        return float(full), float(Z90 * np.sqrt(4 / 80 * ((reps - full) ** 2).sum()))

    def count(self, mask) -> tuple[float, float]:
        mask = np.asarray(mask); full = self.w[mask].sum(); reps = self.rw[mask].sum(0)
        return float(full), float(Z90 * np.sqrt(4 / 80 * ((reps - full) ** 2).sum()))

    def median(self, mask, col: str) -> float:
        mask = np.asarray(mask); v = self.h.loc[mask, col].to_numpy(dtype=float); w = self.w[mask]
        return weighted_median(v, w)


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


def published_tables() -> dict:
    """B08301 as published, by year: workers, worked from home, subway, with margins."""
    out = {"table_note": "B08301: workers 16+ by means of transportation to work; cell 21 is worked from home, cell 12 subway or elevated rail (2019 coding)."}

    def row(total, total_m, home, home_m, subway=None, subway_m=None):
        r = {"workers": total, "workers_moe90": total_m, "worked_from_home": home, "worked_from_home_moe90": home_m, "wfh_share": home / total}
        if subway is not None: r.update({"subway": subway, "subway_moe90": subway_m, "subway_share_workers": subway / total})
        return r

    # 2019: sequence-format summary file, sequence 36, start position 124, 21 cells.
    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 / "e20191ny0036000.txt", header=None, dtype=str, low_memory=False)
    m = pd.read_csv(SF_2019 / "m20191ny0036000.txt", header=None, dtype=str, low_memory=False)
    start = 124 - 1
    y2019 = {}
    for geo, label in GEOS.items():
        if geo == "01000US": continue
        lr = logrec.get(geo)
        if lr is None: continue
        er = e[e[5] == lr].iloc[0]; mr = m[m[5] == lr].iloc[0]
        cell = lambda k: float(er[start + k - 1]); cm = lambda k: float(mr[start + k - 1])
        y2019[label] = row(cell(1), cm(1), cell(21), cm(21), cell(12), cm(12))
    out["2019"] = y2019

    # 2021-2023: table-based summary files.
    for y in (2021, 2022, 2023):
        d = pd.read_csv(REFERENCE / "tables" / f"acsdt1y{y}-b08301.dat", sep="|", dtype=str)
        d = d.set_index("GEO_ID")
        yy = {}
        for geo, label in GEOS.items():
            key = geo.replace("01000US", "0100000US").replace("16000US", "1600000US").replace("05000US", "0500000US")
            if key not in d.index: continue
            r = d.loc[key]
            f = lambda c: float(r[c])
            yy[label] = row(f("B08301_E001"), f("B08301_M001"), f("B08301_E021"), f("B08301_M021"), f("B08301_E012"), f("B08301_M012"))
        out[str(y)] = yy

    # 2024: Census Reporter JSON (city and nation; counties in a second file).
    y2024 = {}
    for path in (TABLES_2024_1YR, TABLES_2024_COUNTIES):
        d = json.loads(path.read_text(encoding="utf-8"))
        for geo, label in GEOS.items():
            if geo not in d["data"]: continue
            b = d["data"][geo]["B08301"]; est, err = b["estimate"], b["error"]
            y2024[label] = row(est["B08301001"], err["B08301001"], est["B08301021"], err["B08301021"], est["B08301012"], err["B08301012"])
    out["2024"] = y2024
    return out


# ------------------------------------------------------------------ estimates


def year_estimates(cs: Census, full: bool) -> tuple[dict, list[dict]]:
    """Everything for one year's NYC workers. `full` adds the two-way cuts the 2019
    and 2024 charts need; the other years get the citywide and borough series only."""
    h = cs.h; nyc = h["nyc"].to_numpy(); rows = []
    wfh = h["wfh"].to_numpy(); cm = h["commuter"].to_numpy(); mode = h["mode"].to_numpy(); boro = h["borough"].to_numpy()
    eb = h["earn_band"].to_numpy(); race = h["race"].to_numpy(); occ = h["occ"].to_numpy(); sec = h["sector"].to_numpy()
    office = h["office_industry"].to_numpy(); ba = h["ba_plus"].to_numpy(); fem = h["female"].to_numpy(); ab = h["age_band"].to_numpy()
    fb = h["foreign_born"].to_numpy(); cow = h["cow"].to_numpy(); powm = h["pow_manhattan"].to_numpy(); toman = h["to_manhattan"].to_numpy()
    res: dict = {"year": int(h["year"].iloc[0])}

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

    def block(mask, extra=False):
        n = int(mask.sum()); c, cmo = cs.count(mask); p, pm = cs.share(mask, wfh); hc, hcm = cs.count(mask & wfh)
        d = {"workers": c, "workers_moe90": cmo, "wfh_share": p, "wfh_share_moe90": pm, "wfh": hc, "wfh_moe90": hcm,
             "commuters": c - hc, "n_records": n}
        if extra and n >= MIN_RECORDS:
            d["subway_share_commuters"] = cs.share(mask & cm, mode == "subway")[0]
            d["car_share_commuters"] = cs.share(mask & cm, mode == "car")[0]
            d["to_manhattan_share"] = cs.share(mask, toman)[0]
            d["median_earnings"] = cs.median(mask, "earnings")
            d["median_earnings_wfh"] = cs.median(mask & wfh, "earnings") if (mask & wfh).sum() >= MIN_RECORDS else float("nan")
            d["median_earnings_commuters"] = cs.median(mask & cm, "earnings")
        return d

    res["all"] = block(nyc, True)
    put("workers", "NYC resident workers at work", res["all"]["workers"], res["all"]["workers_moe90"], nyc.sum(), "count")
    put("wfh_share", "Share of workers who usually worked from home", res["all"]["wfh_share"], res["all"]["wfh_share_moe90"], nyc.sum(), "share")
    put("wfh", "Workers who usually worked from home", res["all"]["wfh"], res["all"]["wfh_moe90"], nyc.sum(), "count")
    res["mode_share_workers"] = {}
    for k, label in MODE_GROUPS:
        p, pm = cs.share(nyc, mode == k); res["mode_share_workers"][k] = {"share": p, "share_moe90": pm, "label": label}
        put(f"mode_{k}", f"Share of workers, {label}", p, pm, nyc.sum(), "share")
    res["mode_share_commuters"] = {k: cs.share(nyc & cm, mode == k)[0] for k, _ in MODE_GROUPS if k != "home"}
    res["boroughs"] = {b: block(nyc & (boro == b), True) for b in BOROUGHS}
    for b in BOROUGHS:
        put(f"wfh_{b}", f"Share working from home, {b}", res["boroughs"][b]["wfh_share"], res["boroughs"][b]["wfh_share_moe90"], (nyc & (boro == b)).sum(), "share")

    # Manhattan as a destination: residents of the city and of the suburbs.
    res["manhattan"] = {"workplace_nyc_residents": cs.count(nyc & powm), "physical_nyc_residents": cs.count(nyc & toman),
                        "wfh_share_nyc_residents": cs.share(nyc & powm, wfh)}
    resid = h["residence"].to_numpy()
    res["manhattan"]["by_residence"] = {}
    for key, label in RESIDENCES:
        mk = (resid == key) & powm
        c, cmo = cs.count(mk); pc, pcm = cs.count(mk & cm); p, pm = cs.share(mk, wfh)
        res["manhattan"]["by_residence"][key] = {"label": label, "workplace": c, "workplace_moe90": cmo, "physical": pc, "physical_moe90": pcm,
                                                 "wfh_share": p, "wfh_share_moe90": pm, "n_records": int(mk.sum()),
                                                 "subway_share": cs.share(mk & cm, mode == "subway")[0] if mk.sum() >= MIN_RECORDS else float("nan"),
                                                 "rail_share": cs.share(mk & cm, mode == "commuter rail")[0] if mk.sum() >= MIN_RECORDS else float("nan"),
                                                 "car_share": cs.share(mk & cm, mode == "car")[0] if mk.sum() >= MIN_RECORDS else float("nan")}
    sub = ~nyc & powm
    res["manhattan"]["physical_suburban"] = cs.count(sub & cm); res["manhattan"]["workplace_suburban"] = cs.count(sub)
    res["manhattan"]["physical_all"] = cs.count(powm & cm); res["manhattan"]["workplace_all"] = cs.count(powm)
    put("manhattan_physical_nyc", "NYC residents physically commuting to a Manhattan workplace", *res["manhattan"]["physical_nyc_residents"], (nyc & toman).sum(), "count")
    put("manhattan_physical_all", "Study-area residents physically commuting to a Manhattan workplace", *res["manhattan"]["physical_all"], (powm & cm).sum(), "count")

    if not full:
        return res, rows

    # One-way cuts.
    res["earnings"] = {k: dict(block(nyc & (eb == k), True), label=l) for k, l, *_ in EARN_BANDS}
    res["occupations"] = {}
    for lab in sorted(set(occ)):
        mk = nyc & (occ == lab)
        if mk.sum() < 300: continue
        res["occupations"][lab] = block(mk, True)
    res["sectors"] = {k: dict(block(nyc & (sec == k), True), label=l) for k, l in SECTORS}
    res["race"] = {r: block(nyc & (race == r), True) for r in RACES}
    res["degree"] = {"ba_plus": block(nyc & ba, True), "no_degree": block(nyc & ~ba, True)}
    res["sex"] = {"women": block(nyc & fem, True), "men": block(nyc & ~fem, True)}
    res["age"] = {k: dict(block(nyc & (ab == k), True), label=l) for k, l, *_ in AGE_BANDS}
    res["nativity"] = {"foreign_born": block(nyc & fb, True), "us_born": block(nyc & ~fb, True)}
    res["class"] = {k: dict(block(nyc & (cow == k), True), label=l) for k, l in CLASSES}
    for k, l, *_ in EARN_BANDS:
        put(f"wfh_earn_{k}", f"Share working from home, own earnings {l}", res["earnings"][k]["wfh_share"], res["earnings"][k]["wfh_share_moe90"], res["earnings"][k]["n_records"], "share")
    for r in RACES:
        put(f"wfh_race_{r}", f"Share working from home, {r}", res["race"][r]["wfh_share"], res["race"][r]["wfh_share_moe90"], res["race"][r]["n_records"], "share")
    put("wfh_ba_plus", "Share working from home, bachelor's or higher", res["degree"]["ba_plus"]["wfh_share"], res["degree"]["ba_plus"]["wfh_share_moe90"], res["degree"]["ba_plus"]["n_records"], "share")
    put("wfh_no_degree", "Share working from home, no bachelor's", res["degree"]["no_degree"]["wfh_share"], res["degree"]["no_degree"]["wfh_share_moe90"], res["degree"]["no_degree"]["n_records"], "share")

    # Office-based industries: who still goes, and how.
    o = nyc & office
    res["office"] = block(o, True)
    res["office"]["modes"] = {k: {"share": cs.share(o, mode == k)[0], "share_moe90": cs.share(o, mode == k)[1], "label": l} for k, l in MODE_GROUPS}
    res["office"]["earnings"] = {k: dict(block(o & (eb == k), True), label=l) for k, l, *_ in EARN_BANDS}
    res["office"]["boroughs"] = {b: block(o & (boro == b), True) for b in BOROUGHS}
    res["office"]["to_manhattan_share_commuters"] = cs.share(o & cm, powm)[0]
    res["office"]["share_of_workers"] = cs.share(nyc, office)[0]
    res["office"]["share_of_wfh"] = cs.share(nyc & wfh, office)[0]
    res["not_office"] = block(nyc & ~office, True)
    put("office_wfh", "Share working from home, office-based industries", res["office"]["wfh_share"], res["office"]["wfh_share_moe90"], o.sum(), "share")
    put("office_commuting", "Share still commuting, office-based industries", 1 - res["office"]["wfh_share"], res["office"]["wfh_share_moe90"], o.sum(), "share")

    # Two-way: earnings within borough, earnings within office industries, race within degree.
    res["earnings_x_borough"] = {b: {k: cs.share(nyc & (boro == b) & (eb == k), wfh)[0] for k, *_ in EARN_BANDS} for b in BOROUGHS}
    res["race_x_degree"] = {r: {"ba_plus": cs.share(nyc & (race == r) & ba, wfh)[0], "no_degree": cs.share(nyc & (race == r) & ~ba, wfh)[0]} for r in RACES}
    res["race_x_office"] = {r: cs.share(nyc & (race == r) & office, wfh)[0] for r in RACES}

    # The remote worker against the commuter, in one profile.
    traits = [("office", "Office-based industry", office), ("ba_plus", "Bachelor's or higher", ba), ("over_100k", "Earns $100k or more", np.isin(eb, ["100k_150k", "150k_plus"])),
              ("under_60k", "Earns under $60k", np.isin(eb, ["lt30k", "30k_60k"])), ("white", "White", race == "White"), ("black_hispanic", "Black or Hispanic", np.isin(race, ["Black", "Hispanic"])),
              ("manhattan", "Lives in Manhattan", boro == "Manhattan"), ("brooklyn", "Lives in Brooklyn", boro == "Brooklyn"), ("women", "Women", fem),
              ("foreign_born", "Foreign-born", fb), ("self_employed", "Self-employed", cow == "self_employed"), ("under_35", "Under 35", np.isin(ab, ["16_24", "25_34"])),
              ("full_time", "Works 35+ hours", h["full_time"].to_numpy())]
    res["profile"] = []
    for key, label, num in traits:
        pw, pwm = cs.share(nyc & wfh, num); pc, pcm = cs.share(nyc & cm, num)
        res["profile"].append({"trait": key, "label": label, "wfh": pw, "wfh_moe90": pwm, "commuters": pc, "commuters_moe90": pcm,
                               "distinguishable": bool(abs(pw - pc) > pwm + pcm)})
    return res, rows


# ------------------------------------------------------------------ the 2026 view: CPS


def cps_view() -> dict | None:
    """NYC residents' telework in the CPS, June 2024 to the latest month: any hours at
    home last week, all of the week's hours at home, some but not all. Pooled periods
    with approximate margins (see the docstring)."""
    if not CPS_FILE.exists():
        return None
    d = pd.read_csv(CPS_FILE)
    d = d[d["weight"] > 0].copy()
    d["period"] = pd.to_datetime(dict(year=d["year"], month=d["month"], day=1))
    hours = d["actual_hours_all"] if "actual_hours_all" in d.columns else d["actual_hours"]   # all jobs, as the hours question is asked
    ah = hours.where(hours > 0)
    th = d["telework_hours"].where(d["telework"], 0).fillna(0)
    d["all_hours"] = d["telework"] & (th >= ah)
    d["some_hours"] = d["telework"] & ~d["all_hours"]
    d["office_occ"] = d["occ_major"].isin([1, 2])            # management, business and financial; professional and related
    d["manhattan"] = d["borough"] == "Manhattan"
    months = sorted(d["period"].unique())
    summary = json.loads(CPS_SUMMARY.read_text(encoding="utf-8")) if CPS_SUMMARY.exists() else {}

    def pooled(mask, num):
        """Weighted mean over months, with the spread of monthly estimates as the error."""
        sub = d[mask]; est = []
        for m in months:
            s = sub[sub["period"] == m]
            if s["weight"].sum() == 0: continue
            est.append((s["weight"] * s[num]).sum() / s["weight"].sum())
        est = np.array(est, float)
        if len(est) == 0: return float("nan"), float("nan"), 0
        se = est.std(ddof=1) / np.sqrt(len(est)) if len(est) > 1 else float("nan")
        return float(est.mean()), float(Z90 * se), int(len(sub))

    periods = {"2024_jun_dec": ("June–December 2024", (d["year"] == 2024)),
               "2025_jan_aug": ("January–August 2025", (d["year"] == 2025) & (d["month"] <= 8)),
               "2025_sep_dec": ("September–December 2025", (d["year"] == 2025) & (d["month"] >= 9)),
               "2026_jan_aug": ("January–August 2026", (d["year"] == 2026))}
    groups = [("all", "All employed NYC residents", np.ones(len(d), bool)), ("ba_plus", "Bachelor's or higher", d["ba_plus"].to_numpy()),
              ("no_degree", "No bachelor's", ~d["ba_plus"].to_numpy()), ("office_occ", "Management, business and professional occupations", d["office_occ"].to_numpy()),
              ("other_occ", "All other occupations", ~d["office_occ"].to_numpy()), ("manhattan", "Lives in Manhattan", d["manhattan"].to_numpy()),
              ("outer", "Lives in the other boroughs", ~d["manhattan"].to_numpy()), ("full_time", "Full-time", d["full_time"].to_numpy()),
              ("private", "Private employer", d["class_of_worker"].isin([4, 5]).to_numpy()), ("government", "Government", d["class_of_worker"].isin([1, 2, 3]).to_numpy()),
              ("self_employed", "Self-employed", d["class_of_worker"].isin([6, 7]).to_numpy())]
    out = {"months": [], "periods": {}, "groups": {}, "question": "At any time last week, did you telework or work at home for pay? (CPS, asked of everyone employed and at work since June 2024)",
           "first_month": months[0].strftime("%Y-%m"), "last_month": months[-1].strftime("%Y-%m"), "n_records": int(len(d)),
           "boroughs_identified": sorted(d["borough"].dropna().unique().tolist()), "summary_notes": summary.get("notes")}
    for m in months:
        s = d[d["period"] == m]; w = s["weight"]
        row = {"month": m.strftime("%Y-%m"), "n": int(len(s)), "employed": float(w.sum()),
               "telework": float((w * s["telework"]).sum() / w.sum()), "all_hours": float((w * s["all_hours"]).sum() / w.sum()),
               "some_hours": float((w * s["some_hours"]).sum() / w.sum())}
        nat = next((x for x in summary.get("months", []) if x.get("label") == row["month"]), None) if isinstance(summary.get("months"), list) else None
        if nat and "national" in nat: row["us_telework"] = nat["national"].get("telework_share"); row["us_all_hours"] = nat["national"].get("all_hours_share")
        out["months"].append(row)
    for pk, (plabel, pmask) in periods.items():
        if pmask.sum() == 0: continue
        block = {"label": plabel, "months": int(d.loc[pmask, "period"].nunique())}
        for stat in ("telework", "all_hours", "some_hours"):
            v, m_, n = pooled(pmask.to_numpy(), stat); block[stat] = v; block[f"{stat}_moe90"] = m_; block["n"] = n
        out["periods"][pk] = block
    for gk, glabel, gmask in groups:
        out["groups"][gk] = {"label": glabel}
        for pk, (plabel, pmask) in periods.items():
            mk = pmask.to_numpy() & gmask
            if mk.sum() < MIN_RECORDS: continue
            v, m_, n = pooled(mk, "telework"); a, am, _ = pooled(mk, "all_hours"); s_, sm, _ = pooled(mk, "some_hours")
            out["groups"][gk][pk] = {"telework": v, "telework_moe90": m_, "all_hours": a, "all_hours_moe90": am, "some_hours": s_, "some_hours_moe90": sm, "n": n}
    return out


# ------------------------------------------------------------------ the 2026 view: MTA by day of week


def mta_day_of_week() -> dict:
    """Average count by day of week and year, against the MTA's own 2019 baseline for
    that day of the week (recovered from the frozen file over 2024)."""
    cur = pd.read_csv(RAW / "transit" / "mta_daily_ridership.csv"); cur["Date"] = pd.to_datetime(cur["Date"])
    old = pd.read_csv(RAW / "transit" / "mta_daily_ridership_2020_2025.csv"); old["Date"] = pd.to_datetime(old["Date"])
    modes = {"Subway": ("Subways: Total Estimated Ridership", "Subways: % of Comparable Pre-Pandemic Day"),
             "Bus": ("Buses: Total Estimated Ridership", "Buses: % of Comparable Pre-Pandemic Day"),
             "LIRR": ("LIRR: Total Estimated Ridership", "LIRR: % of Comparable Pre-Pandemic Day"),
             "MNR": ("Metro-North: Total Estimated Ridership", "Metro-North: % of Comparable Pre-Pandemic Day"),
             "BT": ("Bridges and Tunnels: Total Traffic", "Bridges and Tunnels: % of Comparable Pre-Pandemic Day")}
    old["dow"] = old["Date"].dt.dayofweek
    base = {}
    for m, (c, pc) in modes.items():
        o = old[(old["Date"].dt.year == 2024) & (old[pc] > 0)]
        implied = o[c] / o[pc]
        base[m] = {int(k): float(v) for k, v in implied.groupby(o["dow"]).mean().items()}
    cur["dow"] = cur["Date"].dt.dayofweek; cur["year"] = cur["Date"].dt.year
    cur = cur[cur["Mode"].isin(modes) & (cur["year"] >= 2022)]
    days = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"]
    out = {"last_date": cur["Date"].max().strftime("%Y-%m-%d"), "baseline_note": "2019 baseline per day of week recovered from the MTA's comparable-day percentages over 2024.",
           "years": {}, "days": days, "baseline": {m: {days[k]: v for k, v in b.items()} for m, b in base.items()}}
    for y, g in cur.groupby("year"):
        row = {}
        for m in modes:
            gm = g[g["Mode"] == m]
            by = gm.groupby("dow")["Count"].mean()
            row[m] = {days[int(k)]: {"count": float(v), "vs_2019": float(v / base[m][int(k)])} for k, v in by.items()}
        out["years"][str(int(y))] = row
    # The shape of the office week: Tuesday-Thursday against Monday and Friday, subway, 2026 to date.
    y = out["years"][str(cur["year"].max())]["Subway"]
    out["office_week_2026"] = {"tue_thu": float(np.mean([y[d]["vs_2019"] for d in ("Tuesday", "Wednesday", "Thursday")])),
                               "mon": y["Monday"]["vs_2019"], "fri": y["Friday"]["vs_2019"],
                               "tue_thu_count": float(np.mean([y[d]["count"] for d in ("Tuesday", "Wednesday", "Thursday")])), "fri_count": y["Friday"]["count"]}
    b19 = out["baseline"]["Subway"]
    out["office_week_2019"] = {"tue_thu_count": float(np.mean([b19[d] for d in ("Tuesday", "Wednesday", "Thursday")])), "fri_count": b19["Friday"], "mon_count": b19["Monday"]}
    return out


# ------------------------------------------------------------------ reproduction


def reproduction(pub: dict, years: dict) -> pd.DataFrame:
    rows = []

    def add(label, published, pub_moe, ours, ours_moe, kind):
        band = max(pub_moe if pub_moe == pub_moe else 0, ours_moe if ours_moe == ours_moe else 0)
        rows.append({"statistic": label, "published": published, "published_moe90": pub_moe, "reproduced": ours, "reproduced_moe90": ours_moe,
                     "difference": ours - published, "within_margin": bool(abs(ours - published) <= band or abs(ours - published) <= 0.002), "kind": kind})

    for y in ("2019", "2021", "2022", "2023", "2024"):
        if y not in pub or int(y) not in years: continue
        p = pub[y]; r = years[int(y)]
        if "New York City" in p:
            c = p["New York City"]
            add(f"Workers, NYC, {y}", c["workers"], c["workers_moe90"], r["all"]["workers"], r["all"]["workers_moe90"], "count")
            add(f"Worked from home, NYC, {y}", c["worked_from_home"], c["worked_from_home_moe90"], r["all"]["wfh"], r["all"]["wfh_moe90"], "count")
            add(f"Share working from home, NYC, {y}", c["wfh_share"], np.nan, r["all"]["wfh_share"], r["all"]["wfh_share_moe90"], "share")
            if "subway_share_workers" in c:
                add(f"Subway share of workers, NYC, {y}", c["subway_share_workers"], np.nan, r["mode_share_workers"]["subway"]["share"], r["mode_share_workers"]["subway"]["share_moe90"], "share")
        for b in BOROUGHS:
            if b in p and y in ("2019", "2024"):
                add(f"Share working from home, {b}, {y}", p[b]["wfh_share"], np.nan, r["boroughs"][b]["wfh_share"], r["boroughs"][b]["wfh_share_moe90"], "share")
    return pd.DataFrame(rows)


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


def _mpl():
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    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"})
    return plt


def _style(ax, title, sub):
    ax.set_title(title, loc="left", fontsize=12, color="#1a1a1a", pad=22)
    ax.annotate(sub, xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left")
    for s in ("top", "right", "left"): ax.spines[s].set_visible(False)
    ax.grid(axis="x", color=GRID, linewidth=0.8); ax.set_axisbelow(True); ax.tick_params(length=0)


def _save(fig, path: Path):
    path.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(path, format="svg", bbox_inches="tight"); (path.parent / "pdf").mkdir(exist_ok=True)
    fig.savefig(path.parent / "pdf" / path.with_suffix(".pdf").name, format="pdf", bbox_inches="tight")


def _paired(rows, path, title, sub, xmax, fmt="{:.0f}%", labels=("2019", "2024"), height=None, ticks=None):
    """rows: (label, v2019, m2019, v2024, m2024)."""
    plt = _mpl(); n = len(rows)
    fig, ax = plt.subplots(figsize=(FIG_W, height or (1.2 + 0.5 * n)), dpi=100); fig.patch.set_alpha(0)
    y = np.arange(n)[::-1]; hgt = 0.36
    for j, (col, lab) in enumerate(((BLUE, labels[0]), (ORANGE, labels[1]))):
        v = np.array([r[1 + 2 * j] for r in rows], float); m = np.array([r[2 + 2 * j] if r[2 + 2 * j] == r[2 + 2 * j] else 0 for r in rows], float)
        yy = y + (0.5 - j) * (hgt + 0.03)
        ax.barh(yy, v, height=hgt, color=col, label=lab); ax.errorbar(v, yy, xerr=m, fmt="none", ecolor=TICK, elinewidth=1.2)
        for yi, vi in zip(yy, v): ax.text(vi + xmax * 0.012, yi, fmt.format(vi), va="center", fontsize=8.5, color=INK)
    ax.set_yticks(y); ax.set_yticklabels([r[0] for r in rows]); ax.set_xlim(0, xmax)
    if ticks: ax.set_xticks(ticks[0]); ax.set_xticklabels(ticks[1])
    _style(ax, title, sub); ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.06), ncol=2, frameon=False, fontsize=9, labelcolor=INK)
    fig.tight_layout(); _save(fig, path); plt.close(fig)


def fig_history(years: dict, path):
    plt = _mpl(); ys = [years[y] for y in YEARS if y in years]
    fig, ax = plt.subplots(figsize=(FIG_W, 3.4), dpi=100); fig.patch.set_alpha(0)
    x = [r["year"] for r in ys]; v = [r["all"]["wfh_share"] * 100 for r in ys]; m = [r["all"]["wfh_share_moe90"] * 100 for r in ys]
    ax.plot(x, v, color=BLUE, linewidth=2, marker="o", markersize=4)
    ax.fill_between(x, np.array(v) - np.array(m), np.array(v) + np.array(m), color=BLUE, alpha=0.15, linewidth=0)
    for xi, vi in zip(x, v):
        if xi in (2019, 2021, 2024): ax.text(xi, vi + 1.2, f"{vi:.1f}%", ha="center", fontsize=9, color=INK)
    ax.set_ylim(0, 30); ax.set_yticks([0, 10, 20, 30]); ax.set_yticklabels(["0%", "10%", "20%", "30%"])
    ax.set_xticks(x); ax.set_xticklabels([str(xi) if xi != 2020 else "2020*" for xi in x], fontsize=8.5)
    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("Working from home before, during and after", loc="left", fontsize=12, color="#1a1a1a", pad=22)
    ax.annotate("Share of NYC resident workers who usually worked from home (census, annual; band: 90% margin; *experimental weights)", 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_earnings(y19, y24, path):
    rows = [(l, y19["earnings"][k]["wfh_share"] * 100, y19["earnings"][k]["wfh_share_moe90"] * 100, y24["earnings"][k]["wfh_share"] * 100, y24["earnings"][k]["wfh_share_moe90"] * 100) for k, l, *_ in EARN_BANDS]
    _paired(rows, path, "A pay gradient that did not exist before", "Share of NYC workers who usually worked from home, by own earnings in 2024 dollars (census)", 30,
            ticks=([0, 10, 20, 30], ["0%", "10%", "20%", "30%"]))


def fig_occupations(y19, y24, path):
    occs = sorted(y24["occupations"], key=lambda k: -y24["occupations"][k]["wfh_share"])
    rows = [(k, y19["occupations"].get(k, {}).get("wfh_share", np.nan) * 100, y19["occupations"].get(k, {}).get("wfh_share_moe90", np.nan) * 100,
             y24["occupations"][k]["wfh_share"] * 100, y24["occupations"][k]["wfh_share_moe90"] * 100) for k in occs]
    _paired(rows, path, "Who can work from home", "Share of NYC workers who usually worked from home, by occupation group (census)", 45,
            height=7.5, ticks=([0, 10, 20, 30, 40], ["0%", "10%", "20%", "30%", "40%"]))


def fig_boroughs(y19, y24, path):
    rows = [(b, y19["boroughs"][b]["wfh_share"] * 100, y19["boroughs"][b]["wfh_share_moe90"] * 100, y24["boroughs"][b]["wfh_share"] * 100, y24["boroughs"][b]["wfh_share_moe90"] * 100) for b in BOROUGHS]
    _paired(rows, path, "Where the remote workers live", "Share of resident workers who usually worked from home, by borough (census)", 25,
            ticks=([0, 5, 10, 15, 20, 25], ["0%", "5%", "10%", "15%", "20%", "25%"]))


def fig_office(y19, y24, path):
    modes = [k for k, _ in MODE_GROUPS if k not in ("bicycle", "commuter rail")]
    rows = [(y24["office"]["modes"][k]["label"], y19["office"]["modes"][k]["share"] * 100, y19["office"]["modes"][k]["share_moe90"] * 100,
             y24["office"]["modes"][k]["share"] * 100, y24["office"]["modes"][k]["share_moe90"] * 100) for k in modes]
    _paired(rows, path, "How office-industry workers get to work", "NYC residents in office-based industries, by usual means of transportation (census)", 60,
            ticks=([0, 20, 40, 60], ["0%", "20%", "40%", "60%"]))


def fig_manhattan(y19, y24, path):
    rows = [(y24["manhattan"]["by_residence"][k]["label"], y19["manhattan"]["by_residence"][k]["physical"] / 1000, y19["manhattan"]["by_residence"][k]["physical_moe90"] / 1000,
             y24["manhattan"]["by_residence"][k]["physical"] / 1000, y24["manhattan"]["by_residence"][k]["physical_moe90"] / 1000) for k, _ in RESIDENCES]
    _paired(rows, path, "Who still travels to a Manhattan job", "Workers whose workplace is Manhattan and who did not usually work from home, by where they live (census, thousands)", 800,
            fmt="{:,.0f}k", ticks=([0, 200, 400, 600, 800], ["0", "200k", "400k", "600k", "800k"]))


def fig_transit_dow(dow: dict, path):
    plt = _mpl(); days = dow["days"][:5]; years = [y for y in ("2023", "2024", "2025", "2026") if y in dow["years"]]
    fig, ax = plt.subplots(figsize=(FIG_W, 3.2), dpi=100); fig.patch.set_alpha(0)
    x = np.arange(len(days)); w = 0.8 / len(years); cols = ["#b9c8dd", "#7fa6d6", BLUE, ORANGE]
    for j, y in enumerate(years):
        v = [dow["years"][y]["Subway"][d]["vs_2019"] * 100 for d in days]
        ax.bar(x + (j - (len(years) - 1) / 2) * w, v, width=w, color=cols[j + 4 - len(years)], label=y + (" (to Sep)" if y == "2026" else ""))
    ax.axhline(100, color=TICK, linewidth=0.8, linestyle="--")
    ax.set_xticks(x); ax.set_xticklabels(days); ax.set_ylim(0, 110); ax.set_yticks([0, 25, 50, 75, 100]); ax.set_yticklabels(["0%", "25%", "50%", "75%", "100%"])
    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("The office week, in turnstile counts", loc="left", fontsize=12, color="#1a1a1a", pad=22)
    ax.annotate("Average subway ridership by day of the week as a share of the MTA's 2019 baseline for that day (administrative)", xy=(0, 1), xycoords="axes fraction",
                xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left")
    ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.12), ncol=4, frameon=False, fontsize=9, labelcolor=INK)
    fig.tight_layout(); _save(fig, path); plt.close(fig)


def fig_cps(cps: dict, path):
    plt = _mpl(); ms = cps["months"]
    fig, ax = plt.subplots(figsize=(FIG_W, 3.2), dpi=100); fig.patch.set_alpha(0)
    x = np.arange(len(ms)); a = np.array([m["all_hours"] * 100 for m in ms]); s = np.array([m["some_hours"] * 100 for m in ms])
    ax.bar(x, a, color=BLUE, label="All of last week's hours at home"); ax.bar(x, s, bottom=a, color=ORANGE, label="Some hours at home")
    ax.set_xticks(x[::3]); ax.set_xticklabels([ms[i]["month"] for i in range(0, len(ms), 3)], fontsize=8.5, rotation=0)
    ax.set_ylim(0, 40); ax.set_yticks([0, 10, 20, 30, 40]); ax.set_yticklabels(["0%", "10%", "20%", "30%", "40%"])
    ax.grid(axis="y", color=GRID, linewidth=0.8); ax.grid(axis="x", visible=False)
    for sp in ("top", "right"): ax.spines[sp].set_visible(False)
    ax.set_title("Any hours at home, month by month, to 2026", loc="left", fontsize=12, color="#1a1a1a", pad=22)
    ax.annotate("Employed NYC residents who teleworked or worked at home for pay at any time last week (Current Population Survey, monthly)", xy=(0, 1),
                xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left")
    ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.1), ncol=2, frameon=False, fontsize=9, labelcolor=INK)
    fig.tight_layout(); _save(fig, path); plt.close(fig)


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


def companion_dataset(frames: list[pd.DataFrame], rng: np.random.Generator) -> pd.DataFrame:
    """The census records the 2019-and-2024 comparisons are computed from: one row per
    NYC resident worker at work in each one-year file, with the survey weight and the
    derived variables. Public-use microdata, already disclosure-protected; serial
    numbers are not carried and rows are shuffled."""
    parts = []
    for h in frames:
        h = h[h["nyc"]]
        parts.append(pd.DataFrame({
            "year": h["year"], "weight": h["PWGTP"].astype(int), "borough": h["borough"], "worked_from_home": h["wfh"],
            "mode": h["JWTRNS"].map(MODE_LABEL), "mode_group": h["mode"], "one_way_minutes": h["JWMNP"],
            "workplace_manhattan": h["pow_manhattan"], "physically_commutes_to_manhattan": h["to_manhattan"],
            "own_earnings_2024": h["earnings"].round(0), "earnings_band": h["earn_band"].map({k: l for k, l, *_ in EARN_BANDS}),
            "occupation_group": h["occ"], "industry_sector": h["sector"].map(dict(SECTORS)), "office_based_industry": h["office_industry"],
            "class_of_worker": h["cow"], "usual_hours": h["WKHP"], "race_ethnicity": h["race"], "foreign_born": h["foreign_born"],
            "bachelors_or_higher": h["ba_plus"], "age": h["AGEP"], "sex": h["SEX"].map({1: "male", 2: "female"}),
        }))
    df = pd.concat(parts, ignore_index=True)
    df.insert(0, "row_id", [f"wgo-{i:06d}" for i in 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)
    args = ap.parse_args(argv)
    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"):
        c = pub[y]["New York City"]; print(f"  {y}: workers {c['workers']:,.0f}, worked from home {c['worked_from_home']:,.0f} ({c['wfh_share']:.1%})")

    years: dict[int, dict] = {}; pums_rows = []; frames = {}; censuses = {}
    for y in args.years:
        try:
            h, rw = load_year(y)
        except FileNotFoundError:
            print(f"  {y}: not on disk, skipped"); continue
        cs = Census(h, rw); full = y in FOCUS
        res, rows = year_estimates(cs, full); res["experimental_weights"] = (y == 2020); res["n_records"] = int(h["nyc"].sum())
        years[y] = res; pums_rows += rows
        if full: frames[y] = h; censuses[y] = cs
        a = res["all"]; m = res["manhattan"]
        print(f"  {y}: workers {a['workers']:,.0f}, wfh {a['wfh_share']:.1%} ± {a['wfh_share_moe90']:.1%}, physical Manhattan commuters (NYC residents) "
              f"{m['physical_nyc_residents'][0]:,.0f}, (all study area) {m['physical_all'][0]:,.0f}")

    print("reproduction ...")
    repro = reproduction(pub, years)
    print(repro[["statistic", "published", "reproduced", "within_margin"]].to_string(index=False))
    print("CPS telework, 2024-2026 ...")
    cps = cps_view()
    if cps:
        for pk, b in cps["periods"].items(): print(f"  {b['label']}: any {b['telework']:.1%} ± {b['telework_moe90']:.1%}, all hours {b['all_hours']:.1%}, some {b['some_hours']:.1%} (n {b['n']:,})")
    else:
        print("  (no CPS extract on disk; the 2026 CPS section is skipped)")
    print("MTA by day of week ...")
    mta = mta_series(); dow = mta_day_of_week()
    ow = dow["office_week_2026"]; print(f"  2026 subway vs 2019: Tue-Thu {ow['tue_thu']:.0%}, Monday {ow['mon']:.0%}, Friday {ow['fri']:.0%}")

    y19, y24 = years[2019], years[2024]
    change = {"wfh_share_points": (y24["all"]["wfh_share"] - y19["all"]["wfh_share"]) * 100,
              "wfh_count": y24["all"]["wfh"] - y19["all"]["wfh"],
              "commuters": y24["all"]["commuters"] - y19["all"]["commuters"],
              "subway_share_points": (y24["mode_share_workers"]["subway"]["share"] - y19["mode_share_workers"]["subway"]["share"]) * 100,
              "manhattan_physical_nyc": y24["manhattan"]["physical_nyc_residents"][0] - y19["manhattan"]["physical_nyc_residents"][0],
              "manhattan_physical_all": y24["manhattan"]["physical_all"][0] - y19["manhattan"]["physical_all"][0],
              "manhattan_workplace_all": y24["manhattan"]["workplace_all"][0] - y19["manhattan"]["workplace_all"][0],
              "manhattan_physical_all_pct": y24["manhattan"]["physical_all"][0] / y19["manhattan"]["physical_all"][0] - 1,
              "manhattan_physical_nyc_pct": y24["manhattan"]["physical_nyc_residents"][0] / y19["manhattan"]["physical_nyc_residents"][0] - 1,
              "top_band_points": (y24["earnings"]["150k_plus"]["wfh_share"] - y19["earnings"]["150k_plus"]["wfh_share"]) * 100,
              "bottom_band_points": (y24["earnings"]["lt30k"]["wfh_share"] - y19["earnings"]["lt30k"]["wfh_share"]) * 100}
    res = {"article": SLUG, "as_of": "2026-09-09",
           "definitions": {"universe": "NYC resident workers 16+ at work in the reference week; worked from home is the census's usual means of transportation",
                           "dollar_year": 2024, "deflator": "New York-metro CPI-U (all items) to 2024 dollars", "census_margin": "90% (successive difference replication, 80 replicate weights)",
                           "office_industries": "NAICS sectors 51-56: information, finance and insurance, real estate, professional and technical services, management of companies, administrative and support",
                           "modes": dict(MODE_GROUPS), "earnings_bands": {k: l for k, l, *_ in EARN_BANDS}, "sectors": dict(SECTORS), "residences": dict(RESIDENCES)},
           "published": pub, "years": {str(y): r for y, r in years.items()}, "change_2019_2024": change, "cps": cps, "mta": mta, "mta_day_of_week": dow,
           "context": {"partnership_march_2025": {"average_weekday_attendance": 0.57, "share_of_pre_pandemic": 0.76, "fully_remote": 0.08, "five_days": 0.10, "three_days": 0.30,
                                                  "source": "pfnyc_return_to_office_survey_2025"},
                       "kastle_dec_2025": {"nyc_weekly_occupancy": 0.595, "ten_city": 0.563, "tuesday_ten_city": 0.66, "source": "kastle_back_to_work_barometer_dec_2025"},
                       "manhattan_office_vacancy": {"early_2020": "under 8%", "april_2024": 0.16, "march_2026": 0.131, "sources": ["nyc_comptroller_office_market_spotlight_2024", "propertyshark_nyc_office_q1_2026"]},
                       "atus_2025": {"worked_at_home_on_days_worked": 0.345, "ba_plus": 0.514, "high_school": 0.190, "source": "bls_atus_2025_table6"}},
           "reproduction": repro.to_dict(orient="records"), "seed": args.seed}

    print("figures ...")
    fig_history(years, OUT / "fig1_history.svg")
    fig_earnings(y19, y24, OUT / "fig2_earnings.svg")
    fig_occupations(y19, y24, OUT / "fig3_occupations.svg")
    fig_boroughs(y19, y24, OUT / "fig4_boroughs.svg")
    fig_office(y19, y24, OUT / "fig5_office_modes.svg")
    fig_manhattan(y19, y24, OUT / "fig6_manhattan.svg")
    fig_transit_dow(dow, OUT / "fig7_office_week.svg")
    if cps: fig_cps(cps, OUT / "fig8_cps.svg")

    print("interactive chart cubes ...")
    charts = CH.build(censuses, res["years"], cps, dow)
    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")

    companion_dataset([frames[y] for y in FOCUS], 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"): shutil.copy2(f, SITE_OUT / f.name)
    shutil.copy2(Path(__file__), SITE_OUT / "analysis.py")
    print(f"outputs -> {OUT.relative_to(P.root)} and {SITE_OUT.relative_to(P.root)}")
    return 0


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