diff --git a/packages/populace-build/src/populace/build/us_runtime/US_PUMA_LADDER.md b/packages/populace-build/src/populace/build/us_runtime/US_PUMA_LADDER.md new file mode 100644 index 00000000..edfa00c2 --- /dev/null +++ b/packages/populace-build/src/populace/build/us_runtime/US_PUMA_LADDER.md @@ -0,0 +1,78 @@ +# US PUMA-anchored geography ladder + +The US ACS artifact's geography spine is anchored at the **2020 Public Use +Microdata Area** — the finest geography the ACS PUMS publishes. This is the US +analogue of the UK output-area ladder (populace #354) and a sibling of the US +block ladder (`geography_ladder.py`): one national dataset, filterable at any +grain, never per-area files (populace #275). + +## What it does + +- An **ACS-spine** record already carries its `puma`; the ladder keeps it. +- An **ASEC-spine** record knows only `state_fips`; the ladder draws a PUMA + within that state proportional to 2020 PUMA population. +- Every coarser layer then derives from the PUMA by a population-weighted draw + over the PUMA's block-population overlap with that layer: + - **congressional district (119th)** — `(PUMA, CD)` overlap + - **county** — `(PUMA, county)` overlap + - **tract** (behind `assign_tract`) — `(PUMA, tract)` overlap; when assigned, + the county derives structurally from the tract so the two never disagree. + +Congressional district, county and state are the launch requirement; tract is +the flagged-off stretch rung. Draws come from one seeded generator consumed in +a fixed sorted order (states, then PUMAs), so an assignment is reproducible +from its `seed`. + +The dense ACS-spine artifact is filtered to state / congressional district / +county for local analysis because every record carries all three. + +## Why one block pass + +2020 tabulation blocks nest in 2020 census tracts, and 2020 tracts nest in 2020 +PUMAs. So `tract = block // 10**4`, `county = block // 10**10`, and +`tract → PUMA` is a lookup. Summing 2020 block populations by `(PUMA, CD)`, +`(PUMA, county)`, `(PUMA, tract)` and by PUMA yields the three overlap tables +and the anchor population from a single block pass. Every populated block +contributes to exactly one PUMA, one CD, one county and one tract, so **each +overlap table conserves its PUMA's population exactly** — the loader refuses an +artifact that violates this. + +## Sources (pinned in `us_puma_ladder.provenance.json`) + +| Role | Census source | Vintage | +|---|---|---| +| tract → PUMA (and county via the tract prefix) | 2020 Census Tract to 2020 PUMA relationship file (`2020_Census_Tract_to_2020_PUMA.txt`) | 2020 | +| block → 119th CD | 119th Congressional District BEF (`NationalCD119.txt` in `cd119.zip`) | 119th Congress | +| block → population (the weight) | 2020 P.L. 94-171 geographic headers, per state (`{usps}geo2020.pl`) | 2020 | + +Every source's URL + SHA-256, the reproducible artifact SHA-256, and the +national validation totals are recorded in `us_puma_ladder.provenance.json`. +The `block → CD` and `block → population` parsers are reused unchanged from the +block ladder (`block_ladder_sources.py`); the only source this ladder adds is +the small tract-to-PUMA relationship file. + +## Build + +```bash +uv run python tools/build_us_puma_ladder_artifact.py \ + --out build/us/us_puma_ladder_2020.npz \ + --cache-dir ~/.cache/populace-us-geography +``` + +The tool downloads (and caches) the sources, runs the block pass, aggregates, +writes one national NPZ, and self-checks it by loading it back through +`load_us_puma_ladder` before writing the provenance summary. The artifact is +byte-reproducible from the pinned sources. + +National build (all 51 states): **2,462 PUMAs**, **331,449,281 people** (the +2020 apportionment population, conserved exactly), **436** congressional +districts, **3,143** counties, 83,848 tract overlaps. + +## Wire-up (documented, not rewired) + +The runtime — `load_us_puma_ladder` → `assign_us_puma_ladder` (or +`with_household_us_puma_ladder`) → `us_puma_ladder_gate` — is the drop-in for +the ACS/ASEC spine's geography stage. Wiring it behind the country-spec +geography-spine schema (a PUMA-anchored spine method, and reading the built +artifact from the release bundle) is the follow-up, kept out of this change to +stay reviewable — the same posture #354 took for the UK ladder. diff --git a/packages/populace-build/src/populace/build/us_runtime/__init__.py b/packages/populace-build/src/populace/build/us_runtime/__init__.py index 29a76241..d677677a 100644 --- a/packages/populace-build/src/populace/build/us_runtime/__init__.py +++ b/packages/populace-build/src/populace/build/us_runtime/__init__.py @@ -200,6 +200,25 @@ support_clone_index_column, support_source_id_column, ) +from populace.build.us_runtime.puma_ladder import ( + PUMA_LADDER_ARTIFACT_SHA256_ATTR, + PUMA_LADDER_VINTAGES_ATTR, + US_PUMA_LADDER_COLUMNS, + US_PUMA_LADDER_DERIVED_LAYERS, + US_PUMA_LADDER_KIND, + US_PUMA_LADDER_SCHEMA_VERSION, + US_PUMA_LADDER_TRACT_COLUMN, + UsPumaLadder, + assign_us_puma_ladder, + load_us_puma_ladder, + us_puma_ladder_assignment_summary, + us_puma_ladder_gate, + with_household_us_puma_ladder, +) +from populace.build.us_runtime.puma_ladder_sources import ( + assemble_us_puma_ladder, + parse_tract_to_puma_relationship, +) from populace.build.us_runtime.reform_coverage_smoke import ( us_reform_coverage_smoke_gate, ) @@ -545,6 +564,21 @@ "validation_only_family_ids", "translate_congressional_district_facts_to_current_vintage", "with_household_congressional_districts", + "PUMA_LADDER_ARTIFACT_SHA256_ATTR", + "PUMA_LADDER_VINTAGES_ATTR", + "US_PUMA_LADDER_COLUMNS", + "US_PUMA_LADDER_DERIVED_LAYERS", + "US_PUMA_LADDER_KIND", + "US_PUMA_LADDER_SCHEMA_VERSION", + "US_PUMA_LADDER_TRACT_COLUMN", + "UsPumaLadder", + "assemble_us_puma_ladder", + "assign_us_puma_ladder", + "load_us_puma_ladder", + "parse_tract_to_puma_relationship", + "us_puma_ladder_assignment_summary", + "us_puma_ladder_gate", + "with_household_us_puma_ladder", ] diff --git a/packages/populace-build/src/populace/build/us_runtime/puma_ladder.py b/packages/populace-build/src/populace/build/us_runtime/puma_ladder.py new file mode 100644 index 00000000..9ce5d65b --- /dev/null +++ b/packages/populace-build/src/populace/build/us_runtime/puma_ladder.py @@ -0,0 +1,901 @@ +"""US PUMA-anchored geography-ladder assignment. + +The US ACS artifact's geography spine is anchored at the 2020 Public Use +Microdata Area — the finest geography the ACS PUMS publishes. A household that +already knows its PUMA (an ACS-spine record) keeps it; a household that knows +only its state (an ASEC-spine record) is assigned a PUMA within that state, +sampled proportional to 2020 PUMA population. Every coarser layer of the +ladder then derives from the PUMA by a population-weighted draw over the +PUMA's overlap with that layer: + +- **congressional district (119th)** — sampled within the PUMA proportional to + the ``(PUMA, CD)`` block-population overlap (a 2020 PUMA can span several + districts and a district several PUMAs; they do not nest). +- **county** — sampled within the PUMA proportional to the ``(PUMA, county)`` + block-population overlap. +- **tract** (behind ``assign_tract``) — sampled within the PUMA proportional to + the ``(PUMA, tract)`` block-population overlap; when tract is assigned the + county derives structurally from it (``county = tract // 10**6``) so the two + never disagree. Congressional district, county and state are the launch + requirement; tract is the stretch rung. + +One national dataset, filter by geography, at any grain (populace #275; no +per-area files, the standing rule): the dense ACS-spine artifact is filtered to +state / congressional district / county for local analysis because every record +carries them. + +Vintage discipline follows the country-spec geography-spine schema +(``vintage_policy: "error"``): the artifact records one vintage per derived +layer, the loader refuses an artifact that omits any of them, and the +assignment refuses a ladder whose congressional-district vintage differs from +the vintage the caller expects — a silent partial join on mismatched geography +vintages is the failure these checks exist to forbid (the #205 lesson). + +Column names follow the policyengine-us household input surface +(``county_fips``, ``tract_geoid``, ``congressional_district_geoid``); the +2020-PUMA geoid rides as a ``puma`` data column (policyengine-us has no PUMA +input, exactly as the UK ladder carries ``*_code`` columns policyengine-uk has +no input for). The exported artifact carries engine inputs, never formula +outputs, so ``in_nyc``/``nyc_income_tax`` recompute from the county rung +instead of being persisted (the #34 regression). +""" + +from __future__ import annotations + +import json +from collections.abc import Mapping +from dataclasses import dataclass +from pathlib import Path +from typing import Any + +import numpy as np +import pandas as pd + +from populace.build.gates import GateResult +from populace.build.us_runtime.congressional_district_geography import ( + CONGRESSIONAL_DISTRICT_GEOID_COLUMN, +) +from populace.build.us_runtime.geography_ladder import US_NYC_COUNTY_FIPS +from populace.build.us_runtime.puma_ladder_sources import ( + COUNTY_FROM_TRACT_DIVISOR, + PUMA_GEOID_STATE_DIVISOR, +) +from populace.frame import Frame + +#: Derived layers the ladder artifact must carry a vintage for. +US_PUMA_LADDER_DERIVED_LAYERS = ( + "congressional_district", + "county", + "tract", +) + +#: Household spine columns the launch assignment writes (state / CD / county +#: filtering surface), in write order. ``puma`` is the anchor; the other two +#: are policyengine-us household inputs. +US_PUMA_LADDER_COLUMNS = ( + "puma", + CONGRESSIONAL_DISTRICT_GEOID_COLUMN, + "county_fips", +) + +#: The finer tract rung, written only when ``assign_tract=True``. +US_PUMA_LADDER_TRACT_COLUMN = "tract_geoid" + +US_PUMA_LADDER_SCHEMA_VERSION = 1 +US_PUMA_LADDER_KIND = "us_puma_ladder" + +PUMA_LADDER_ARTIFACT_SHA256_ATTR = "populace_puma_ladder_artifact_sha256" +PUMA_LADDER_VINTAGES_ATTR = "populace_puma_ladder_vintages" + +_ANCHOR_KEYS = ("puma", "puma_population") +_OVERLAP_KEYS = { + "congressional_district": ( + "cd_overlap_puma", + "cd_overlap_cd", + "cd_overlap_population", + ), + "county": ( + "county_overlap_puma", + "county_overlap_county", + "county_overlap_population", + ), + "tract": ( + "tract_overlap_puma", + "tract_overlap_tract", + "tract_overlap_population", + ), +} +_REQUIRED_ARRAY_KEYS = ( + *_ANCHOR_KEYS, + *(key for keys in _OVERLAP_KEYS.values() for key in keys), +) + + +@dataclass(frozen=True) +class UsPumaLadder: + """The national PUMA ladder: one anchor row per populated 2020 PUMA plus + three ``(PUMA, layer_value, population)`` overlap tables sorted by + ``(puma, layer_value)``. + + Attributes: + puma: 2020 PUMA geoids (``state_fips * 10**5 + PUMA5CE``) as ``int64``, + unique and sorted — the anchor rung and the ASEC state -> PUMA + sampling frame. + puma_population: 2020 census population per PUMA (the state -> PUMA + sampling weight), strictly positive. + cd_overlap_puma / cd_overlap_cd / cd_overlap_population: the + ``(PUMA, congressional district)`` population overlap. CD geoids use + the populace ``state_fips * 100 + district`` convention; at-large + and delegate districts are district ``00``. + county_overlap_puma / county_overlap_county / county_overlap_population: + the ``(PUMA, county)`` population overlap; county fips are the + 5-digit ``state+county`` integer. + tract_overlap_puma / tract_overlap_tract / tract_overlap_population: + the ``(PUMA, tract)`` population overlap; tract geoids are the + 11-digit ``state+county+tract`` integer. + metadata: the artifact's embedded metadata, including one vintage per + derived layer. + """ + + puma: np.ndarray + puma_population: np.ndarray + cd_overlap_puma: np.ndarray + cd_overlap_cd: np.ndarray + cd_overlap_population: np.ndarray + county_overlap_puma: np.ndarray + county_overlap_county: np.ndarray + county_overlap_population: np.ndarray + tract_overlap_puma: np.ndarray + tract_overlap_tract: np.ndarray + tract_overlap_population: np.ndarray + metadata: Mapping[str, Any] + + def __len__(self) -> int: + return len(self.puma) + + @property + def layer_vintages(self) -> dict[str, str]: + """Vintage per derived layer, plus the PUMA anchor's own vintage.""" + layers = self.metadata["layers"] + vintages = { + layer: str(layers[layer]["vintage"]) + for layer in US_PUMA_LADDER_DERIVED_LAYERS + } + vintages["puma"] = str(self.metadata["puma_vintage"]) + return vintages + + def _overlap(self, layer: str) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + puma_key, value_key, population_key = _OVERLAP_KEYS[layer] + return ( + getattr(self, puma_key), + getattr(self, value_key), + getattr(self, population_key), + ) + + +def load_us_puma_ladder(path: str | Path) -> UsPumaLadder: + """Load and validate a US PUMA-ladder artifact (NPZ). + + A single national file (never per-area files) built by + ``tools/build_us_puma_ladder_artifact.py`` from primary Census sources. + Every validation failure raises — a ladder that loads is a ladder every + assignment invariant holds for, including exact per-PUMA population + conservation of each overlap table. + """ + + source = Path(path) + if not source.exists(): + raise FileNotFoundError(f"US PUMA ladder artifact not found: {source}") + with np.load(source, allow_pickle=False) as payload: + missing = [key for key in _REQUIRED_ARRAY_KEYS if key not in payload.files] + if missing: + raise ValueError( + f"US PUMA ladder artifact is missing required key(s): {missing}." + ) + if "metadata_json" not in payload.files: + raise ValueError( + "US PUMA ladder artifact is missing required key 'metadata_json'." + ) + arrays = {key: np.asarray(payload[key]) for key in _REQUIRED_ARRAY_KEYS} + metadata = _metadata_from_scalar(payload["metadata_json"]) + + _validate_ladder_metadata(metadata) + + puma = _int64_array(arrays["puma"], label="puma") + if len(puma) == 0: + raise ValueError("US PUMA ladder artifact has zero PUMAs.") + puma_state = puma // PUMA_GEOID_STATE_DIVISOR + if ((puma_state < 1) | (puma_state > 56)).any(): + bad = puma[(puma_state < 1) | (puma_state > 56)][:5].tolist() + raise ValueError( + "puma geoids must be state_fips*10**5+PUMA5CE with state in " + f"[1, 56]; invalid: {bad}." + ) + if len(np.unique(puma)) != len(puma): + raise ValueError("puma geoids must be unique.") + if (np.diff(puma) <= 0).any(): + raise ValueError("puma geoids must be sorted ascending.") + + population = np.asarray(arrays["puma_population"], dtype=np.float64) + if len(population) != len(puma): + raise ValueError("puma_population must align with puma.") + if not np.isfinite(population).all() or (population <= 0).any(): + raise ValueError("puma_population must be positive and finite for every PUMA.") + + puma_set = set(puma.tolist()) + cd_puma = _validated_overlap_puma(arrays["cd_overlap_puma"], puma_set, layer="cd") + cd_value = _int64_array(arrays["cd_overlap_cd"], label="cd_overlap_cd") + cd_pop = _overlap_population(arrays["cd_overlap_population"], cd_puma, layer="cd") + if ((cd_value < 100) | (cd_value > 9999)).any(): + bad = cd_value[(cd_value < 100) | (cd_value > 9999)][:5].tolist() + raise ValueError( + f"cd_overlap_cd must be state_fips*100+district; invalid: {bad}." + ) + _assert_layer_state_matches( + cd_puma, cd_value // 100, layer="congressional_district" + ) + _assert_conservation( + puma, population, cd_puma, cd_pop, layer="congressional_district" + ) + + county_puma = _validated_overlap_puma( + arrays["county_overlap_puma"], puma_set, layer="county" + ) + county_value = _int64_array( + arrays["county_overlap_county"], label="county_overlap_county" + ) + county_pop = _overlap_population( + arrays["county_overlap_population"], county_puma, layer="county" + ) + if ((county_value < 1001) | (county_value > 56999)).any(): + bad = county_value[(county_value < 1001) | (county_value > 56999)][:5].tolist() + raise ValueError(f"county_overlap_county must be 5-digit fips; invalid: {bad}.") + _assert_layer_state_matches(county_puma, county_value // 1000, layer="county") + _assert_conservation(puma, population, county_puma, county_pop, layer="county") + + tract_puma = _validated_overlap_puma( + arrays["tract_overlap_puma"], puma_set, layer="tract" + ) + tract_value = _int64_array( + arrays["tract_overlap_tract"], label="tract_overlap_tract" + ) + tract_pop = _overlap_population( + arrays["tract_overlap_population"], tract_puma, layer="tract" + ) + if ((tract_value < 10**9) | (tract_value >= 10**11)).any(): + bad = tract_value[(tract_value < 10**9) | (tract_value >= 10**11)][:5].tolist() + raise ValueError( + f"tract_overlap_tract must be 11-digit geoids; invalid: {bad}." + ) + _assert_layer_state_matches( + tract_puma, tract_value // COUNTY_FROM_TRACT_DIVISOR // 1000, layer="tract" + ) + if len(np.unique(tract_value)) != len(tract_value): + raise ValueError("tract_overlap_tract geoids must be unique (tracts nest).") + _assert_conservation(puma, population, tract_puma, tract_pop, layer="tract") + + return UsPumaLadder( + puma=puma, + puma_population=population, + cd_overlap_puma=cd_puma, + cd_overlap_cd=cd_value, + cd_overlap_population=cd_pop, + county_overlap_puma=county_puma, + county_overlap_county=county_value.astype(np.int32), + county_overlap_population=county_pop, + tract_overlap_puma=tract_puma, + tract_overlap_tract=tract_value, + tract_overlap_population=tract_pop, + metadata=metadata, + ) + + +def assign_us_puma_ladder( + household: pd.DataFrame, + ladder: UsPumaLadder, + *, + seed: int = 0, + assign_tract: bool = False, + expected_congressional_district_vintage: str | None = None, + state_fips_column: str = "state_fips", + puma_column: str = "puma", +) -> pd.DataFrame: + """Assign each household its PUMA anchor and the derived ladder columns. + + Records that already carry a valid ``puma`` (the ACS spine) keep it; records + that carry none (the ASEC spine) draw a PUMA within their ``state_fips`` + proportional to 2020 PUMA population. Every record then draws a congressional + district and a county within its PUMA proportional to the block-population + overlap, and — when ``assign_tract`` — a tract from which the county derives + structurally. Draws are consumed from one seeded generator in a fixed sorted + order (states, then PUMAs), so the result is reproducible from ``seed``. + + A household ``puma`` absent from the ladder, or a household ``state_fips`` + the ladder has no PUMAs for, is an error — never a silent partial join. + """ + + if state_fips_column not in household.columns: + raise ValueError( + f"household table must contain {state_fips_column!r} before PUMA-" + "ladder assignment." + ) + if expected_congressional_district_vintage is not None: + ladder_cd_vintage = ladder.layer_vintages["congressional_district"] + if ladder_cd_vintage != expected_congressional_district_vintage: + raise ValueError( + "US PUMA ladder congressional-district vintage " + f"{ladder_cd_vintage!r} does not match the vintage the caller " + f"expects ({expected_congressional_district_vintage!r})." + ) + + state_strings = _state_fips_strings(household[state_fips_column]) + state_ints = state_strings.astype(np.int64) + n = len(household) + rng = np.random.default_rng(seed) + + row_puma = np.full(n, -1, dtype=np.int64) + if puma_column in household.columns: + known_mask, known_puma = _parse_input_puma( + household[puma_column], ladder, state_ints + ) + row_puma[known_mask] = known_puma[known_mask] + unknown_mask = row_puma < 0 + if unknown_mask.any(): + _draw_puma_within_state(row_puma, unknown_mask, state_ints, ladder, rng) + + cd_puma, cd_value, cd_pop = ladder._overlap("congressional_district") + cd_values = _draw_layer_values( + row_puma, cd_puma, cd_value, cd_pop, rng, layer="congressional_district" + ) + + tract_values: np.ndarray | None = None + if assign_tract: + tract_puma, tract_value, tract_pop = ladder._overlap("tract") + tract_values = _draw_layer_values( + row_puma, tract_puma, tract_value, tract_pop, rng, layer="tract" + ) + county_values = tract_values // COUNTY_FROM_TRACT_DIVISOR + else: + county_puma, county_value, county_pop = ladder._overlap("county") + county_values = _draw_layer_values( + row_puma, county_puma, county_value, county_pop, rng, layer="county" + ) + + _assert_row_state_consistency( + state_ints, + puma_state=row_puma // PUMA_GEOID_STATE_DIVISOR, + cd_state=cd_values // 100, + county_state=county_values // 1000, + ) + + assigned = household.copy() + assigned["puma"] = np.array( + [f"{value:07d}" for value in row_puma.tolist()], dtype=object + ) + assigned[CONGRESSIONAL_DISTRICT_GEOID_COLUMN] = cd_values.astype(np.int64) + assigned["county_fips"] = np.array( + [f"{value:05d}" for value in county_values.tolist()], dtype=object + ) + if tract_values is not None: + assigned[US_PUMA_LADDER_TRACT_COLUMN] = np.array( + [f"{value:011d}" for value in tract_values.tolist()], dtype=object + ) + return assigned + + +def with_household_us_puma_ladder( + frame: Frame, + ladder: UsPumaLadder, + *, + seed: int = 0, + assign_tract: bool = False, + expected_congressional_district_vintage: str | None = None, +) -> Frame: + """Return ``frame`` with the household PUMA ladder assigned.""" + + tables = {entity: frame.table(entity).copy() for entity in frame.entities} + tables["household"] = assign_us_puma_ladder( + tables["household"], + ladder, + seed=seed, + assign_tract=assign_tract, + expected_congressional_district_vintage=( + expected_congressional_district_vintage + ), + ) + for link_name in frame.links: + tables[link_name] = frame.link(link_name).copy() + weights = {entity: frame.weights_for(entity) for entity in frame.weighted_entities} + return Frame( + tables, + frame.schema, + weights, + frame.strata, + mass_log=frame.mass_log, + ) + + +def us_puma_ladder_assignment_summary( + household: pd.DataFrame, + ladder: UsPumaLadder, + *, + weight_values: np.ndarray | None = None, + assign_tract: bool = False, +) -> dict[str, Any]: + """Summarize ladder assignment coverage for provenance logs.""" + + columns = list(US_PUMA_LADDER_COLUMNS) + if assign_tract: + columns.append(US_PUMA_LADDER_TRACT_COLUMN) + applied = all(column in household.columns for column in columns) + weights = ( + np.asarray(weight_values, dtype=np.float64) + if weight_values is not None + else np.ones(len(household), dtype=np.float64) + ) + if len(weights) != len(household): + raise ValueError( + f"weight_values length {len(weights)} does not match household " + f"rows {len(household)}." + ) + summary: dict[str, Any] = { + "applied": applied, + "household_rows": int(len(household)), + "ladder_pumas": int(len(ladder)), + "layer_vintages": ladder.layer_vintages, + "sampling_basis": str(ladder.metadata["sampling_basis"]), + } + if not applied: + return summary + total = float(weights.sum()) + for column in columns: + values = household[column].astype(str).to_numpy() + nonempty = values != "" + summary[f"assigned_{column}_values"] = int( + pd.Series(values[nonempty]).nunique() + ) + summary[f"{column}_nonempty_weighted_share"] = ( + float(weights[nonempty].sum() / total) if total > 0 else 0.0 + ) + county = household["county_fips"].astype(str).to_numpy() + nyc = np.isin(county, US_NYC_COUNTY_FIPS) + summary["nyc_weighted_household_share"] = ( + float(weights[nyc].sum() / total) if total > 0 else 0.0 + ) + return summary + + +def us_puma_ladder_gate( + household: pd.DataFrame, + weight_values: np.ndarray, + *, + assign_tract: bool = False, + state_fips_column: str = "state_fips", + nyc_share_bounds: tuple[float, float] = (0.005, 0.06), + nyc_share_of_new_york_bounds: tuple[float, float] = (0.20, 0.65), +) -> GateResult: + """Gate the exported PUMA ladder: structure, state consistency, and NYC mass. + + The NYC checks are the permanent form of the populace #34 regression: + recomputing ``in_nyc`` from the PUMA-derived county rung must never silently + collapse NYC to zero. NYC holds roughly 2.6% of national household weight and + roughly 42% of New York State's, so the default bounds fail on collapse (0%) + and on gross misassignment, not on calibration noise. + """ + + failures: list[str] = [] + details: dict[str, object] = {} + weights = np.asarray(weight_values, dtype=np.float64) + if len(weights) != len(household): + raise ValueError( + f"weight_values length {len(weights)} does not match household " + f"rows {len(household)}." + ) + + expected_columns = [*US_PUMA_LADDER_COLUMNS, state_fips_column] + if assign_tract: + expected_columns.append(US_PUMA_LADDER_TRACT_COLUMN) + missing_columns = [ + column for column in expected_columns if column not in household.columns + ] + if missing_columns: + failures.append( + f"household table is missing geography column(s): {missing_columns}" + ) + return GateResult( + name="us_puma_ladder", + passed=False, + failures=tuple(failures), + details=details, + ) + + state = _state_fips_strings(household[state_fips_column]) + puma = household["puma"].astype(str).to_numpy() + county = household["county_fips"].astype(str).to_numpy() + for label, values, width in (("puma", puma, 7), ("county_fips", county, 5)): + bad = ~(np.char.str_len(values.astype(np.str_)) == width) | ~np.char.isdigit( + values.astype(np.str_) + ) + if bad.any(): + failures.append( + f"{label}: {int(bad.sum())}/{len(values)} row(s) are not " + f"{width}-digit codes; examples {sorted(set(values[bad]))[:5]}" + ) + + prefix_checks = [ + ("puma", np.array([value[:2] for value in puma.tolist()]), state), + ("county_fips", np.array([value[:2] for value in county.tolist()]), state), + ] + cd_state = _household_cd_state(household[CONGRESSIONAL_DISTRICT_GEOID_COLUMN]) + prefix_checks.append((CONGRESSIONAL_DISTRICT_GEOID_COLUMN, cd_state, state)) + if assign_tract: + tract = household[US_PUMA_LADDER_TRACT_COLUMN].astype(str).to_numpy() + bad_tract = ~(np.char.str_len(tract.astype(np.str_)) == 11) | ~np.char.isdigit( + tract.astype(np.str_) + ) + if bad_tract.any(): + failures.append( + f"tract_geoid: {int(bad_tract.sum())}/{len(tract)} row(s) are " + f"not 11-digit codes; examples {sorted(set(tract[bad_tract]))[:5]}" + ) + prefix_checks.append( + ("tract_geoid", np.array([value[:5] for value in tract.tolist()]), county) + ) + for label, values, expected in prefix_checks: + mismatched = values != expected + if mismatched.any(): + failures.append( + f"{label}: {int(mismatched.sum())}/{len(values)} row(s) " + "disagree with the expected state/county prefix" + ) + + total = float(weights.sum()) + if total <= 0: + failures.append("household weights sum to zero; shares are undefined") + return GateResult( + name="us_puma_ladder", + passed=False, + failures=tuple(failures), + details=details, + ) + + nyc_mask = np.isin(county, US_NYC_COUNTY_FIPS) + nyc_share = float(weights[nyc_mask].sum() / total) + details["nyc_weighted_household_share"] = nyc_share + low, high = nyc_share_bounds + if not low <= nyc_share <= high: + failures.append( + f"NYC weighted household share {nyc_share:.4f} outside " + f"[{low}, {high}] (the #34 in_nyc collapse regression)" + ) + new_york = state == "36" + new_york_total = float(weights[new_york].sum()) + if new_york_total > 0: + nyc_of_ny = float(weights[nyc_mask & new_york].sum() / new_york_total) + details["nyc_weighted_share_of_new_york"] = nyc_of_ny + low, high = nyc_share_of_new_york_bounds + if not low <= nyc_of_ny <= high: + failures.append( + f"NYC weighted share of New York State {nyc_of_ny:.4f} " + f"outside [{low}, {high}]" + ) + else: + failures.append("New York State carries zero household weight") + + return GateResult( + name="us_puma_ladder", + passed=not failures, + failures=tuple(failures), + details=details, + ) + + +def _parse_input_puma( + values: Any, ladder: UsPumaLadder, state_ints: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """Split household ``puma`` into a known-mask and integer geoids. + + Null / non-positive entries are ASEC records to be drawn later. A positive + entry must be a PUMA the ladder covers, in the household's own state — an + absent PUMA is a vintage mismatch, never a silent partial join. + """ + + series = pd.Series(values).reset_index(drop=True) + numeric = pd.to_numeric(series, errors="coerce").to_numpy(dtype=np.float64) + known = np.isfinite(numeric) & (numeric > 0) + non_integer = known & (numeric != np.floor(numeric)) + if non_integer.any(): + bad = series.to_numpy()[non_integer][:5].tolist() + raise ValueError(f"puma contains non-integer geoid value(s): {bad}.") + puma_int = np.full(len(series), -1, dtype=np.int64) + puma_int[known] = numeric[known].astype(np.int64) + + ladder_set = set(ladder.puma.tolist()) + missing = sorted(set(puma_int[known].tolist()) - ladder_set) + if missing: + raise ValueError( + "US PUMA ladder has no anchor row for household puma geoid(s): " + f"{missing[:10]} ({len(missing)} total). The household PUMAs and the " + "ladder artifact must share a PUMA vintage." + ) + puma_state = puma_int[known] // PUMA_GEOID_STATE_DIVISOR + want_state = state_ints[known] + mismatched = puma_state != want_state + if mismatched.any(): + examples = [ + f"state {want:02d} -> puma {value:07d}" + for want, value in zip( + want_state[mismatched][:5].tolist(), + puma_int[known][mismatched][:5].tolist(), + strict=True, + ) + ] + raise ValueError(f"household puma state disagrees with state_fips: {examples}.") + return known, puma_int + + +def _draw_puma_within_state( + row_puma: np.ndarray, + unknown_mask: np.ndarray, + state_ints: np.ndarray, + ladder: UsPumaLadder, + rng: np.random.Generator, +) -> None: + """Draw a PUMA within state (proportional to PUMA population) in place.""" + + ladder_state = ladder.puma // PUMA_GEOID_STATE_DIVISOR + unknown_positions = np.nonzero(unknown_mask)[0] + unknown_states = state_ints[unknown_mask] + for state in np.unique(unknown_states): + rows = unknown_positions[unknown_states == state] + in_state = ladder_state == state + state_pumas = ladder.puma[in_state] + if len(state_pumas) == 0: + raise ValueError( + f"US PUMA ladder has no PUMAs for state_fips {int(state):02d}; " + "the ASEC household state and the ladder must share coverage." + ) + weights = ladder.puma_population[in_state].astype(np.float64) + row_puma[rows] = rng.choice( + state_pumas, size=len(rows), replace=True, p=weights / weights.sum() + ) + + +def _draw_layer_values( + row_puma: np.ndarray, + overlap_puma: np.ndarray, + overlap_value: np.ndarray, + overlap_population: np.ndarray, + rng: np.random.Generator, + *, + layer: str, +) -> np.ndarray: + """Draw one layer value per row proportional to its PUMA's overlap. + + ``overlap_*`` are sorted by PUMA, so each PUMA's slice is contiguous; + households are grouped by PUMA in sorted order for reproducibility. + """ + + n = len(row_puma) + out = np.full(n, -1, dtype=np.int64) + positions = pd.Series(np.arange(n, dtype=np.int64)) + for puma, group_positions in positions.groupby(pd.Series(row_puma), sort=True): + lo = int(np.searchsorted(overlap_puma, puma, side="left")) + hi = int(np.searchsorted(overlap_puma, puma, side="right")) + if hi <= lo: + raise ValueError( + f"US PUMA ladder has no {layer} overlap for PUMA {int(puma)}." + ) + values = overlap_value[lo:hi] + weights = overlap_population[lo:hi].astype(np.float64) + rows = group_positions.to_numpy() + out[rows] = rng.choice( + values, size=len(rows), replace=True, p=weights / weights.sum() + ) + if (out < 0).any(): + raise ValueError(f"internal error: some rows received no {layer} draw.") + return out + + +def _assert_row_state_consistency( + state_ints: np.ndarray, + *, + puma_state: np.ndarray, + cd_state: np.ndarray, + county_state: np.ndarray, +) -> None: + for label, values in ( + ("puma", puma_state), + (CONGRESSIONAL_DISTRICT_GEOID_COLUMN, cd_state), + ("county_fips", county_state), + ): + mismatched = values != state_ints + if mismatched.any(): + examples = [ + f"state {want:02d} -> {label} state {got:02d}" + for want, got in zip( + state_ints[mismatched][:5].tolist(), + values[mismatched][:5].tolist(), + strict=True, + ) + ] + raise ValueError( + f"assigned {label} disagrees with household state_fips " + f"(ladder is inconsistent): {examples}." + ) + + +def _metadata_from_scalar(value: np.ndarray) -> dict[str, Any]: + if value.shape != (): + raise ValueError("metadata_json must be a scalar JSON string.") + raw = value.item() + if isinstance(raw, bytes): + raw = raw.decode() + if not isinstance(raw, str): + raise ValueError("metadata_json must be a scalar JSON string.") + metadata = json.loads(raw) + if not isinstance(metadata, dict): + raise ValueError("metadata_json must decode to an object.") + return metadata + + +def _validate_ladder_metadata(metadata: Mapping[str, Any]) -> None: + if metadata.get("schema_version") != US_PUMA_LADDER_SCHEMA_VERSION: + raise ValueError( + "US PUMA ladder metadata schema_version must be " + f"{US_PUMA_LADDER_SCHEMA_VERSION}, got {metadata.get('schema_version')!r}." + ) + if metadata.get("kind") != US_PUMA_LADDER_KIND: + raise ValueError( + f"US PUMA ladder metadata kind must be {US_PUMA_LADDER_KIND!r}, " + f"got {metadata.get('kind')!r}." + ) + if not str(metadata.get("puma_vintage") or ""): + raise ValueError("US PUMA ladder metadata must record puma_vintage.") + if not str(metadata.get("sampling_basis") or ""): + raise ValueError("US PUMA ladder metadata must record sampling_basis.") + layers = metadata.get("layers") + if not isinstance(layers, Mapping): + raise ValueError("US PUMA ladder metadata must carry a layers mapping.") + missing = [layer for layer in US_PUMA_LADDER_DERIVED_LAYERS if layer not in layers] + if missing: + raise ValueError( + f"US PUMA ladder metadata layers are missing: {missing}. Every " + "derived layer records its own vintage (vintage_policy: error)." + ) + for layer in US_PUMA_LADDER_DERIVED_LAYERS: + spec = layers[layer] + if not isinstance(spec, Mapping) or not str(spec.get("vintage") or ""): + raise ValueError( + f"US PUMA ladder layer {layer!r} must record a non-empty vintage." + ) + if not str(spec.get("source") or ""): + raise ValueError( + f"US PUMA ladder layer {layer!r} must record a non-empty " + "source citation." + ) + + +def _validated_overlap_puma( + values: np.ndarray, puma_set: set[int], *, layer: str +) -> np.ndarray: + array = _int64_array(values, label=f"{layer}_overlap_puma") + if len(array) == 0: + raise ValueError(f"{layer} overlap table is empty.") + if (np.diff(array) < 0).any(): + raise ValueError(f"{layer}_overlap_puma must be sorted ascending.") + unknown = sorted(set(array.tolist()) - puma_set) + if unknown: + raise ValueError( + f"{layer} overlap references PUMA(s) absent from the anchor: {unknown[:5]}." + ) + return array + + +def _overlap_population( + values: np.ndarray, aligned: np.ndarray, *, layer: str +) -> np.ndarray: + array = np.asarray(values, dtype=np.float64) + if len(array) != len(aligned): + raise ValueError( + f"{layer}_overlap_population must align with {layer}_overlap_puma." + ) + if not np.isfinite(array).all() or (array <= 0).any(): + raise ValueError(f"{layer} overlap population must be positive and finite.") + return array + + +def _assert_layer_state_matches( + overlap_puma: np.ndarray, value_state: np.ndarray, *, layer: str +) -> None: + puma_state = overlap_puma // PUMA_GEOID_STATE_DIVISOR + mismatched = puma_state != value_state + if mismatched.any(): + examples = [ + f"puma {p} -> {layer} state {s}" + for p, s in zip( + overlap_puma[mismatched][:5].tolist(), + value_state[mismatched][:5].tolist(), + strict=True, + ) + ] + raise ValueError( + f"{layer} overlap state prefix must match its PUMA's state; " + f"mismatched: {examples}." + ) + + +def _assert_conservation( + puma: np.ndarray, + population: np.ndarray, + overlap_puma: np.ndarray, + overlap_population: np.ndarray, + *, + layer: str, +) -> None: + anchor = dict(zip(puma.tolist(), population.tolist(), strict=True)) + totals: dict[int, float] = {int(key): 0.0 for key in puma.tolist()} + for key, value in zip( + overlap_puma.tolist(), overlap_population.tolist(), strict=True + ): + totals[key] += value + for key, total in totals.items(): + if abs(total - anchor[key]) > 1e-6: + raise ValueError( + f"{layer} overlap for PUMA {key} sums to {total} but the anchor " + f"population is {anchor[key]} (must conserve exactly)." + ) + + +def _int64_array(values: np.ndarray, *, label: str) -> np.ndarray: + array = np.asarray(values) + if array.dtype.kind not in ("i", "u"): + raise ValueError(f"{label} must be an integer array, got dtype {array.dtype}.") + return array.astype(np.int64, copy=False) + + +def _state_fips_strings(values: Any) -> np.ndarray: + series = pd.Series(values) + numeric = pd.to_numeric(series, errors="coerce") + if numeric.isna().any(): + bad = series[numeric.isna()].head(5).tolist() + raise ValueError(f"state_fips contains non-numeric value(s): {bad}.") + numeric_values = numeric.to_numpy(dtype=np.float64) + if (numeric_values != np.floor(numeric_values)).any(): + raise ValueError("state_fips must contain integer values.") + integers = numeric_values.astype(np.int64) + if ((integers < 1) | (integers > 99)).any(): + raise ValueError( + "state_fips must contain two-digit state FIPS codes; invalid: " + f"{integers[(integers < 1) | (integers > 99)][:5].tolist()}." + ) + return np.array([f"{value:02d}" for value in integers.tolist()]) + + +def _household_cd_state(values: Any) -> np.ndarray: + numeric = pd.to_numeric(pd.Series(values), errors="coerce") + if numeric.isna().any(): + bad = pd.Series(values)[numeric.isna()].head(5).tolist() + raise ValueError( + f"{CONGRESSIONAL_DISTRICT_GEOID_COLUMN} contains non-numeric " + f"value(s): {bad}." + ) + integers = numeric.to_numpy(dtype=np.int64) + return np.array([f"{value // 100:02d}" for value in integers.tolist()]) + + +__all__ = [ + "PUMA_LADDER_ARTIFACT_SHA256_ATTR", + "PUMA_LADDER_VINTAGES_ATTR", + "US_PUMA_LADDER_COLUMNS", + "US_PUMA_LADDER_DERIVED_LAYERS", + "US_PUMA_LADDER_KIND", + "US_PUMA_LADDER_SCHEMA_VERSION", + "US_PUMA_LADDER_TRACT_COLUMN", + "UsPumaLadder", + "assign_us_puma_ladder", + "load_us_puma_ladder", + "us_puma_ladder_assignment_summary", + "us_puma_ladder_gate", + "with_household_us_puma_ladder", +] diff --git a/packages/populace-build/src/populace/build/us_runtime/puma_ladder_sources.py b/packages/populace-build/src/populace/build/us_runtime/puma_ladder_sources.py new file mode 100644 index 00000000..a35f324b --- /dev/null +++ b/packages/populace-build/src/populace/build/us_runtime/puma_ladder_sources.py @@ -0,0 +1,265 @@ +"""Parsers and assembly for the primary Census sources behind the US PUMA ladder. + +The PUMA ladder anchors every household at a 2020 Public Use Microdata Area — +the finest geography the ACS PUMS publishes — and derives every coarser layer +(congressional district, county, tract) from that anchor by a +population-weighted draw over the PUMA's overlap with each layer. One national +dataset, filterable at any grain (populace #275; no per-area files). + +Because 2020 tabulation blocks nest in 2020 census tracts, and 2020 tracts nest +in 2020 PUMAs, every overlap distribution the ladder needs comes from a single +block pass joined to one small tract-to-PUMA relationship file: + +- ``tract = block_geoid // 10**4`` and ``county = block_geoid // 10**10`` are + structural prefixes of the 15-digit block geoid. +- ``tract -> PUMA`` comes from the Census 2020 Census Tract to 2020 PUMA + relationship file (each 2020 tract belongs to exactly one 2020 PUMA). +- ``block -> 119th CD`` and ``block -> POP100`` reuse the block-ladder parsers + :func:`parse_national_cd_bef` and :func:`parse_pl_geo_blocks` unchanged. + +Summing block populations by ``(PUMA, CD)``, ``(PUMA, county)`` and +``(PUMA, tract)`` yields the three overlap tables; summing by PUMA yields the +anchor population used for the ASEC state -> PUMA draw. Every populated block +contributes to exactly one PUMA, one CD, one county and one tract, so each +overlap table conserves its PUMA's population exactly — a defect (a populated +block the CD file does not cover, or a tract the relationship file omits) fails +loudly here rather than shipping silent partial weights. + +Download orchestration lives in ``tools/build_us_puma_ladder_artifact.py``; +everything here is pure and unit-testable. +""" + +from __future__ import annotations + +import json +from collections.abc import Iterable, Mapping +from typing import Any + +import numpy as np + +#: PUMA geoids are ``state_fips * 10**5 + PUMA5CE`` (a 2-digit state prefix and +#: the 5-digit 2020 PUMA code), matching how block/tract geoids drop leading +#: zeros to an integer. Tract geoids are the 11-digit ``state+county+tract`` +#: integer; county fips the 5-digit ``state+county`` integer. +PUMA_GEOID_STATE_DIVISOR = 10**5 +TRACT_FROM_BLOCK_DIVISOR = 10**4 +COUNTY_FROM_BLOCK_DIVISOR = 10**10 +COUNTY_FROM_TRACT_DIVISOR = 10**6 + +_TRACT_TO_PUMA_HEADER = ("STATEFP", "COUNTYFP", "TRACTCE", "PUMA5CE") + + +def parse_tract_to_puma_relationship( + lines: Iterable[str], + *, + allowed_state_fips: frozenset[str] | None = None, +) -> dict[int, int]: + """Parse the 2020 tract-to-PUMA relationship file into tract geoid -> PUMA. + + The file is comma-delimited with header ``STATEFP,COUNTYFP,TRACTCE,PUMA5CE`` + (a UTF-8 BOM may prefix the first field). Rows for states outside + ``allowed_state_fips`` (Puerto Rico and the island territories, when the + caller passes the 50-states-plus-DC set) are skipped so the returned map + matches the block spine. Every populated tract maps to exactly one PUMA; + a tract that appears twice with conflicting PUMAs is a source defect. + """ + + iterator = iter(lines) + header = _required_header(iterator, source="tract-to-PUMA relationship") + cleaned = header.lstrip("") + if tuple(part.strip().upper() for part in cleaned.split(",")) != ( + _TRACT_TO_PUMA_HEADER + ): + raise ValueError( + "tract-to-PUMA relationship header must be " + f"'STATEFP,COUNTYFP,TRACTCE,PUMA5CE', got {header!r}." + ) + result: dict[int, int] = {} + for line_number, line in enumerate(iterator, start=2): + stripped = line.strip() + if not stripped: + continue + parts = stripped.split(",") + if len(parts) != 4: + raise ValueError( + f"tract-to-PUMA line {line_number} must have four fields, " + f"got {stripped!r}." + ) + state_fips, county_fips, tract_ce, puma_ce = (part.strip() for part in parts) + _require_digits(state_fips, width=2, source=f"tract-to-PUMA line {line_number}") + if allowed_state_fips is not None and state_fips not in allowed_state_fips: + continue + _require_digits( + county_fips, width=3, source=f"tract-to-PUMA line {line_number}" + ) + _require_digits(tract_ce, width=6, source=f"tract-to-PUMA line {line_number}") + _require_digits(puma_ce, width=5, source=f"tract-to-PUMA line {line_number}") + tract = int(f"{state_fips}{county_fips}{tract_ce}") + puma = int(f"{state_fips}{puma_ce}") + existing = result.get(tract) + if existing is not None and existing != puma: + raise ValueError( + f"tract-to-PUMA relationship assigns tract " + f"{state_fips}{county_fips}{tract_ce} to both PUMA " + f"{existing} and {puma}." + ) + result[tract] = puma + if not result: + raise ValueError("tract-to-PUMA relationship contained no rows.") + return result + + +def assemble_us_puma_ladder( + *, + block_population: Mapping[int, int], + cd_by_block: Mapping[int, int], + tract_to_puma: Mapping[int, int], + metadata: Mapping[str, Any], +) -> dict[str, np.ndarray]: + """Aggregate the block pass into the PUMA ladder artifact's NPZ payload. + + For every populated block the block's tract (a structural prefix) must map + to a PUMA and the block must carry a congressional district; either gap is + a source defect, not a skippable row, so the three overlap tables each + conserve their PUMA's population exactly. Returns arrays for the anchor + (``puma`` / ``puma_population``) and the three overlap tables, each sorted + by ``(puma, layer_value)`` for a stable, searchsorted-friendly artifact. + """ + + puma_population: dict[int, int] = {} + puma_cd_population: dict[tuple[int, int], int] = {} + puma_county_population: dict[tuple[int, int], int] = {} + puma_tract_population: dict[tuple[int, int], int] = {} + missing_puma: list[int] = [] + missing_cd: list[int] = [] + + for block, population in block_population.items(): + if population <= 0: + continue + tract = block // TRACT_FROM_BLOCK_DIVISOR + puma = tract_to_puma.get(tract) + if puma is None: + missing_puma.append(block) + continue + cd = cd_by_block.get(block) + if cd is None: + missing_cd.append(block) + continue + county = block // COUNTY_FROM_BLOCK_DIVISOR + puma_population[puma] = puma_population.get(puma, 0) + population + cd_key = (puma, cd) + puma_cd_population[cd_key] = puma_cd_population.get(cd_key, 0) + population + county_key = (puma, county) + puma_county_population[county_key] = ( + puma_county_population.get(county_key, 0) + population + ) + tract_key = (puma, tract) + puma_tract_population[tract_key] = ( + puma_tract_population.get(tract_key, 0) + population + ) + + if missing_puma: + examples = [f"{block:015d}" for block in sorted(missing_puma)[:5]] + raise ValueError( + f"{len(missing_puma)} populated block(s) have a tract absent from " + f"the tract-to-PUMA relationship file; examples: {examples}." + ) + if missing_cd: + examples = [f"{block:015d}" for block in sorted(missing_cd)[:5]] + raise ValueError( + f"{len(missing_cd)} populated block(s) have no congressional " + f"district in the CD BEF; examples: {examples}." + ) + if not puma_population: + raise ValueError("PUMA ladder assembly produced no populated PUMAs.") + + pumas = np.asarray(sorted(puma_population), dtype=np.int64) + anchor_population = np.asarray( + [puma_population[puma] for puma in pumas.tolist()], dtype=np.int64 + ) + cd_puma, cd_value, cd_pop = _overlap_arrays(puma_cd_population) + county_puma, county_value, county_pop = _overlap_arrays(puma_county_population) + tract_puma, tract_value, tract_pop = _overlap_arrays(puma_tract_population) + + _assert_conserves( + pumas, anchor_population, cd_puma, cd_pop, layer="congressional_district" + ) + _assert_conserves(pumas, anchor_population, county_puma, county_pop, layer="county") + _assert_conserves(pumas, anchor_population, tract_puma, tract_pop, layer="tract") + + return { + "puma": pumas, + "puma_population": anchor_population, + "cd_overlap_puma": cd_puma, + "cd_overlap_cd": cd_value, + "cd_overlap_population": cd_pop, + "county_overlap_puma": county_puma, + "county_overlap_county": county_value.astype(np.int32), + "county_overlap_population": county_pop, + "tract_overlap_puma": tract_puma, + "tract_overlap_tract": tract_value, + "tract_overlap_population": tract_pop, + "metadata_json": np.asarray(json.dumps(dict(metadata), sort_keys=True)), + } + + +def _overlap_arrays( + overlap: Mapping[tuple[int, int], int], +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Return ``(puma, value, population)`` arrays sorted by ``(puma, value)``.""" + + keys = sorted(overlap) + puma = np.asarray([key[0] for key in keys], dtype=np.int64) + value = np.asarray([key[1] for key in keys], dtype=np.int64) + population = np.asarray([overlap[key] for key in keys], dtype=np.int64) + return puma, value, population + + +def _assert_conserves( + pumas: np.ndarray, + anchor_population: np.ndarray, + overlap_puma: np.ndarray, + overlap_population: np.ndarray, + *, + layer: str, +) -> None: + """Every PUMA's overlap population must sum to its anchor population.""" + + totals = {int(puma): 0 for puma in pumas.tolist()} + for puma, population in zip( + overlap_puma.tolist(), overlap_population.tolist(), strict=True + ): + if puma not in totals: + raise ValueError( + f"{layer} overlap references PUMA {puma} absent from the anchor." + ) + totals[puma] += population + for puma, anchor in zip(pumas.tolist(), anchor_population.tolist(), strict=True): + if totals[puma] != anchor: + raise ValueError( + f"{layer} overlap for PUMA {puma} sums to {totals[puma]:,} but " + f"the anchor population is {anchor:,}." + ) + + +def _required_header(iterator: Iterable[str], *, source: str) -> str: + for line in iterator: + stripped = line.strip() + if stripped: + return stripped + raise ValueError(f"{source} is empty.") + + +def _require_digits(value: str, *, width: int, source: str) -> None: + if not (value.isdigit() and len(value) == width): + raise ValueError(f"{source}: expected a {width}-digit code, got {value!r}.") + + +__all__ = [ + "COUNTY_FROM_BLOCK_DIVISOR", + "COUNTY_FROM_TRACT_DIVISOR", + "PUMA_GEOID_STATE_DIVISOR", + "TRACT_FROM_BLOCK_DIVISOR", + "assemble_us_puma_ladder", + "parse_tract_to_puma_relationship", +] diff --git a/packages/populace-build/src/populace/build/us_runtime/us_puma_ladder.provenance.json b/packages/populace-build/src/populace/build/us_runtime/us_puma_ladder.provenance.json new file mode 100644 index 00000000..e6bb64c5 --- /dev/null +++ b/packages/populace-build/src/populace/build/us_runtime/us_puma_ladder.provenance.json @@ -0,0 +1,284 @@ +{ + "artifact": "us_puma_ladder_2020.npz", + "congressional_district_overlaps": 4043, + "congressional_districts": 436, + "counties": 3143, + "county_overlaps": 4620, + "layer_vintages": { + "congressional_district": "119th_congress", + "county": "2020_census", + "puma": "2020_puma", + "tract": "2020_census" + }, + "output_sha256": "39a2ab2abeab07a88362af7ab2940e0e1d50a297c919e4bbc6fb65bab51147d8", + "population": 331449281, + "pumas": 2462, + "source_files": { + "cd119_bef": { + "sha256": "1433feb5178dc7b4188ee30f5f7f715851f4400740b8fe1ce606a876c6294bd6", + "url": "https://www2.census.gov/programs-surveys/decennial/rdo/mapping-files/2025/119-congressional-district-befs/cd119.zip" + }, + "pl94171_ak": { + "sha256": "c90837c12ab2f5fe99132c558705e8d0b5174afdf2a0d28b71927f535fe63487", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Alaska/ak2020.pl.zip" + }, + "pl94171_al": { + "sha256": "bd3fbbf22615b166869f9119aa12017b9a7b82baf0f22a9170e92a88a397cd77", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Alabama/al2020.pl.zip" + }, + "pl94171_ar": { + "sha256": "940ad2916676e443dfdab98c3e145fcf62ef9dcd9f10c763ccdb195f4d5f132e", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Arkansas/ar2020.pl.zip" + }, + "pl94171_az": { + "sha256": "56997e8b1b5175148fd3e1157ab57eee7c1722221cafebfd5d12fe2e687b2aff", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Arizona/az2020.pl.zip" + }, + "pl94171_ca": { + "sha256": "e02c83c85a9a5d69ded4bf0d4c8ebe3f4bcd2d8e5911db765e9e42915ad1f753", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/California/ca2020.pl.zip" + }, + "pl94171_co": { + "sha256": "fce76008aedd10681bf0e57fa366ca788d280264edfe91e7d80a9b33ae86ab2c", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Colorado/co2020.pl.zip" + }, + "pl94171_ct": { + "sha256": "c4a2d42a79c52de10ea9397cca43db59c2df38e6d6e66b4396a0de8bd3f66780", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Connecticut/ct2020.pl.zip" + }, + "pl94171_dc": { + "sha256": "21b5189f45ce260c4607c9579f4ad52d11920133af827b14b183121b48d40ce4", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/District_of_Columbia/dc2020.pl.zip" + }, + "pl94171_de": { + "sha256": "2d7743cc442574dc785c728e9637c5d4648bef518beef41b5f03756eb56ea78a", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Delaware/de2020.pl.zip" + }, + "pl94171_fl": { + "sha256": "fdcc98180f634b69c40ed1d71a63fe00fdda4cda3e35fbaa6382336501f4b0cf", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Florida/fl2020.pl.zip" + }, + "pl94171_ga": { + "sha256": "26434588ceab435fa33a138ff999e4187c923a74732fcc488c3a0597b4adebcb", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Georgia/ga2020.pl.zip" + }, + "pl94171_hi": { + "sha256": "b7e1ba6350c5d75267b688504052af001ad227763995f6bb26ccdd92fc45261f", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Hawaii/hi2020.pl.zip" + }, + "pl94171_ia": { + "sha256": "c55e8276080b684431b8629269d4225c31080ea17257a737f9d99e3cf2ac6822", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Iowa/ia2020.pl.zip" + }, + "pl94171_id": { + "sha256": "89345b6d5d12ccdec8af4039b0a6f1f5ba4cd7ca86c8f73ad6250b5713e3df47", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Idaho/id2020.pl.zip" + }, + "pl94171_il": { + "sha256": "e73f82d6aea0b9b372e9f8572c6edbedfa72ae9fc2039da84faf0cdf0c478a32", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Illinois/il2020.pl.zip" + }, + "pl94171_in": { + "sha256": "5f4999237a9bf3b3337865494d0d62c0a1dc620cad7e48caebaa83811f95901a", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Indiana/in2020.pl.zip" + }, + "pl94171_ks": { + "sha256": "ac6b4ae5e2af4c8d639b46748d55927ee7821756f6073cc4f8238c1fb00498ee", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Kansas/ks2020.pl.zip" + }, + "pl94171_ky": { + "sha256": "252a04b760e4da69665d5654108be2b63d55d294c1c0c608bfc7e7fcf2360ac4", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Kentucky/ky2020.pl.zip" + }, + "pl94171_la": { + "sha256": "d8d9ea96195f2b848be70ef33dcb224d433d8d58025e11aaa0432c5211e5345c", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Louisiana/la2020.pl.zip" + }, + "pl94171_ma": { + "sha256": "9857362fbc1162d3943685134a77a91542a70eb667f3477559d07fd0a81e4889", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Massachusetts/ma2020.pl.zip" + }, + "pl94171_md": { + "sha256": "9af1e1f0185b6219a9432132ef78e47c2d74c474948d2ed6004c7061412393f3", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Maryland/md2020.pl.zip" + }, + "pl94171_me": { + "sha256": "f9bb6fe4cf13bb28b44e4d6a8d5d3e5209f9249ebc402db7ea5c24f8a673ac07", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Maine/me2020.pl.zip" + }, + "pl94171_mi": { + "sha256": "971bd53abeb1d905bb9b09bfe4dc1afe8514a916f24d285b289e0f66ec5cfb62", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Michigan/mi2020.pl.zip" + }, + "pl94171_mn": { + "sha256": "7d30974597e73893d25ffef0975c8e673afcc1166f9919d5ff2a466a2e0e188f", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Minnesota/mn2020.pl.zip" + }, + "pl94171_mo": { + "sha256": "900d8a3bcc59e2cc3c174730e2a8704bd653af1fc043c9fbd73196e1c0d7fba5", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Missouri/mo2020.pl.zip" + }, + "pl94171_ms": { + "sha256": "cb64e3465cee9623dc574605f4dfcada33eaea6833062d31f769811fd20a5f26", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Mississippi/ms2020.pl.zip" + }, + "pl94171_mt": { + "sha256": "e01288ce3ae37ae732eb61d40268c6ac3ec8f9772c8609c6e9ca538bdd0bc130", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Montana/mt2020.pl.zip" + }, + "pl94171_nc": { + "sha256": "04f2f250dad877d237d94f553f9cb2209594c2a021d15f9fdeedd4404db538db", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/North_Carolina/nc2020.pl.zip" + }, + "pl94171_nd": { + "sha256": "a37d18e05a73b311a8fd7ade36dd5b17cac16ff125318c9ec6d2edd7f6968b71", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/North_Dakota/nd2020.pl.zip" + }, + "pl94171_ne": { + "sha256": "29298d68a5905cf9cb14cd64e35750d14727a330c5940ccf5238a8ae6e110300", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Nebraska/ne2020.pl.zip" + }, + "pl94171_nh": { + "sha256": "31852c0cfa871ea81f933bb1df904f0d45997097e38ade83db87ba58255d4468", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/New_Hampshire/nh2020.pl.zip" + }, + "pl94171_nj": { + "sha256": "caf0da4a9e9d46ced3f8020a7a2bf5240cd8805ed0679716f0cb2b8fcc65e811", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/New_Jersey/nj2020.pl.zip" + }, + "pl94171_nm": { + "sha256": "77133ba4bb2fb2c28d39993bb4cff5f1c42823c9d13b75e51e8ecf986354c259", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/New_Mexico/nm2020.pl.zip" + }, + "pl94171_nv": { + "sha256": "ac7e792e0a1d6248ff1700159f4ec6a3a7b3330ebe10899bb4f47ab97c19952d", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Nevada/nv2020.pl.zip" + }, + "pl94171_ny": { + "sha256": "9d046014bee6b61afa808c5dfe841a6d036c77bff70197472fbdd291a40b8ce2", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/New_York/ny2020.pl.zip" + }, + "pl94171_oh": { + "sha256": "ceef338247ad2dccc267885f2879e167ab6730c7a332d3ceb5aef9d6e2338aa0", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Ohio/oh2020.pl.zip" + }, + "pl94171_ok": { + "sha256": "e95c99e6385f4ed3a5ea67962444a53ad7a731f26d1dd88df2efde504f101dff", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Oklahoma/ok2020.pl.zip" + }, + "pl94171_or": { + "sha256": "f8e7fa8c331446850b1bf545e6c4aa63733373f582b11edd6b3cf31016bc4509", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Oregon/or2020.pl.zip" + }, + "pl94171_pa": { + "sha256": "2d33a7dab29c8dd5692bbde203d253e06eebbc44fcbaa96b1caa958d454026ae", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Pennsylvania/pa2020.pl.zip" + }, + "pl94171_ri": { + "sha256": "a2a6de20d0bac0dfe91aef66d3c2dde0eb768863884fbc45d4f3df5c49f5dcd0", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Rhode_Island/ri2020.pl.zip" + }, + "pl94171_sc": { + "sha256": "e8996f6873405d6c048c42e8878231fc7f152564295498ff586a041648a46bb5", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/South_Carolina/sc2020.pl.zip" + }, + "pl94171_sd": { + "sha256": "30d4edc3ccc86ad76a2a388c763ecdaebe23f43a5dcb92658587c04ca5cfa1a5", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/South_Dakota/sd2020.pl.zip" + }, + "pl94171_tn": { + "sha256": "89f843dbc717409169fef44df2682548c308c06c0d825451bf2e78340234daf2", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Tennessee/tn2020.pl.zip" + }, + "pl94171_tx": { + "sha256": "6fee8ca5759eb4942db25a59b8cbd3439ad4188279fccff87dac70a89ceca6ff", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Texas/tx2020.pl.zip" + }, + "pl94171_ut": { + "sha256": "4e6de13cd7b5307bfdd327446dd7107e2fdb790903f59186bb8285881531d0f0", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Utah/ut2020.pl.zip" + }, + "pl94171_va": { + "sha256": "e33c7c1ac0aa84b5fe5ba20188d2e562c5c96b1c02f59e766e18d037b23df9e8", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Virginia/va2020.pl.zip" + }, + "pl94171_vt": { + "sha256": "6a18a3364dff5e0432b25f4b4972eaf22ac0699ea0273430509a9b6a763638e9", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Vermont/vt2020.pl.zip" + }, + "pl94171_wa": { + "sha256": "dde77f2db632a4af274a968c39a7ab58fac6cfd90718e4a72ec239b84980db5b", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Washington/wa2020.pl.zip" + }, + "pl94171_wi": { + "sha256": "e70f80ef4a803cdeb1e1767fc42e13f7d01ce9f98566c82ca4c2379aa9809f7f", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Wisconsin/wi2020.pl.zip" + }, + "pl94171_wv": { + "sha256": "8356c49f300e9c18658de56c06e551f5e8cb591c395fc3d1748077d46322889b", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/West_Virginia/wv2020.pl.zip" + }, + "pl94171_wy": { + "sha256": "98dd7d035942fea2758e0fcabe71efa0588a7ec7c9c33110bcde493fe8e8b942", + "url": "https://www2.census.gov/programs-surveys/decennial/2020/data/01-Redistricting_File--PL_94-171/Wyoming/wy2020.pl.zip" + }, + "tract_to_puma": { + "sha256": "5262f460b2c8dbe86b00549e916317b1791fdc595c9af72ab18f6f6a67040531", + "url": "https://www2.census.gov/geo/docs/maps-data/data/rel2020/2020_Census_Tract_to_2020_PUMA.txt" + } + }, + "states": [ + "01", + "02", + "04", + "05", + "06", + "08", + "09", + "10", + "11", + "12", + "13", + "15", + "16", + "17", + "18", + "19", + "20", + "21", + "22", + "23", + "24", + "25", + "26", + "27", + "28", + "29", + "30", + "31", + "32", + "33", + "34", + "35", + "36", + "37", + "38", + "39", + "40", + "41", + "42", + "44", + "45", + "46", + "47", + "48", + "49", + "50", + "51", + "53", + "54", + "55", + "56" + ], + "tract_overlaps": 83848 +} diff --git a/packages/populace-build/tests/test_us_puma_ladder.py b/packages/populace-build/tests/test_us_puma_ladder.py new file mode 100644 index 00000000..3bba90e3 --- /dev/null +++ b/packages/populace-build/tests/test_us_puma_ladder.py @@ -0,0 +1,467 @@ +"""US PUMA-anchored geography-ladder tests.""" + +import json +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from populace.build.us_runtime import ( + US_PUMA_LADDER_COLUMNS, + assign_us_puma_ladder, + load_us_puma_ladder, + us_puma_ladder_assignment_summary, + us_puma_ladder_gate, + with_household_us_puma_ladder, +) +from populace.frame import US_SCHEMA, Frame, WeightKind, Weights + + +def _ladder_metadata(**overrides: object) -> dict: + layers = { + "congressional_district": { + "vintage": "119th_congress", + "source": "Census 119th Congressional District BEF (NationalCD119.txt)", + }, + "county": { + "vintage": "2020_census", + "source": "Census 2020 Census Tract to 2020 PUMA relationship file", + }, + "tract": { + "vintage": "2020_census", + "source": "Census 2020 Census Tract to 2020 PUMA relationship file", + }, + } + metadata = { + "schema_version": 1, + "kind": "us_puma_ladder", + "puma_vintage": "2020_puma", + "sampling_basis": "population", + "layers": layers, + } + metadata.update(overrides) + return metadata + + +def _write_ladder(path: Path, **overrides: object) -> Path: + # PUMA 0100100 spans CDs 101/102 (900/100) and counties 01001/01003; + # PUMA 0100200 is CD 102 in county 01003; PUMA 0200100 is the state-02 + # at-large district. Every overlap table conserves its PUMA population. + arrays: dict[str, np.ndarray] = { + "puma": np.asarray([100100, 100200, 200100], dtype=np.int64), + "puma_population": np.asarray([1000.0, 500.0, 400.0]), + "cd_overlap_puma": np.asarray([100100, 100100, 100200, 200100], dtype=np.int64), + "cd_overlap_cd": np.asarray([101, 102, 102, 200], dtype=np.int64), + "cd_overlap_population": np.asarray([900.0, 100.0, 500.0, 400.0]), + "county_overlap_puma": np.asarray( + [100100, 100100, 100200, 200100], dtype=np.int64 + ), + "county_overlap_county": np.asarray([1001, 1003, 1003, 2013], dtype=np.int64), + "county_overlap_population": np.asarray([900.0, 100.0, 500.0, 400.0]), + "tract_overlap_puma": np.asarray( + [100100, 100100, 100200, 200100], dtype=np.int64 + ), + "tract_overlap_tract": np.asarray( + [1001000100, 1003000100, 1003000200, 2013000100], dtype=np.int64 + ), + "tract_overlap_population": np.asarray([900.0, 100.0, 500.0, 400.0]), + "metadata_json": np.asarray(json.dumps(_ladder_metadata())), + } + arrays.update({k: np.asarray(v) for k, v in overrides.items()}) + np.savez_compressed(path, **arrays) + return path + + +def _ladder(tmp_path: Path): + return load_us_puma_ladder(_write_ladder(tmp_path / "puma_ladder.npz")) + + +# --------------------------------------------------------------------------- # +# Loading and validation +# --------------------------------------------------------------------------- # + + +def test_load_round_trips_arrays_and_vintages(tmp_path) -> None: + ladder = _ladder(tmp_path) + + assert len(ladder) == 3 + assert ladder.layer_vintages == { + "puma": "2020_puma", + "congressional_district": "119th_congress", + "county": "2020_census", + "tract": "2020_census", + } + assert ladder.puma_population.dtype == np.float64 + assert ladder.county_overlap_county.dtype == np.int32 + + +def test_load_refuses_missing_array_key(tmp_path) -> None: + path = tmp_path / "puma_ladder.npz" + _write_ladder(path) + with np.load(path) as payload: + arrays = {k: payload[k] for k in payload.files if k != "cd_overlap_cd"} + np.savez_compressed(path, **arrays) + + with pytest.raises(ValueError, match="missing required key"): + load_us_puma_ladder(path) + + +def test_load_refuses_missing_layer_vintage(tmp_path) -> None: + metadata = _ladder_metadata() + del metadata["layers"]["tract"] + path = _write_ladder( + tmp_path / "puma_ladder.npz", metadata_json=np.asarray(json.dumps(metadata)) + ) + + with pytest.raises(ValueError, match="vintage_policy: error"): + load_us_puma_ladder(path) + + +def test_load_refuses_nonpositive_population(tmp_path) -> None: + path = _write_ladder( + tmp_path / "puma_ladder.npz", + puma_population=np.asarray([1000.0, 0.0, 400.0]), + ) + + with pytest.raises(ValueError, match="positive"): + load_us_puma_ladder(path) + + +def test_load_refuses_duplicate_or_unsorted_pumas(tmp_path) -> None: + path = _write_ladder( + tmp_path / "puma_ladder.npz", + puma=np.asarray([100200, 100100, 200100], dtype=np.int64), + ) + + with pytest.raises(ValueError, match="sorted ascending"): + load_us_puma_ladder(path) + + +def test_load_refuses_overlap_not_conserving_population(tmp_path) -> None: + path = _write_ladder( + tmp_path / "puma_ladder.npz", + cd_overlap_population=np.asarray([800.0, 100.0, 500.0, 400.0]), + ) + + with pytest.raises(ValueError, match="must conserve exactly"): + load_us_puma_ladder(path) + + +def test_load_refuses_overlap_puma_absent_from_anchor(tmp_path) -> None: + path = _write_ladder( + tmp_path / "puma_ladder.npz", + cd_overlap_puma=np.asarray([100100, 100100, 100200, 999999], dtype=np.int64), + ) + + with pytest.raises(ValueError, match="absent from the anchor"): + load_us_puma_ladder(path) + + +def test_load_refuses_cd_state_mismatch(tmp_path) -> None: + path = _write_ladder( + tmp_path / "puma_ladder.npz", + cd_overlap_cd=np.asarray([101, 102, 102, 600], dtype=np.int64), + ) + + with pytest.raises(ValueError, match="state prefix must match"): + load_us_puma_ladder(path) + + +# --------------------------------------------------------------------------- # +# Assignment — ACS spine (known PUMA) and ASEC spine (drawn PUMA) +# --------------------------------------------------------------------------- # + + +def test_acs_rows_keep_their_puma_and_derive_the_ladder(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame( + { + "household_id": [10, 20], + "state_fips": [1, 2], + "puma": [100100, 200100], + }, + index=[100, 200], + ) + + assigned = assign_us_puma_ladder(household, ladder, seed=0) + + assert assigned.index.tolist() == [100, 200] + for column in US_PUMA_LADDER_COLUMNS: + assert column in assigned.columns + # The state-02 at-large PUMA has exactly one CD and one county: fully + # determined regardless of seed. + assert assigned.loc[200, "puma"] == "0200100" + assert assigned.loc[200, "congressional_district_geoid"] == 200 + assert assigned.loc[200, "county_fips"] == "02013" + # The state-01 ACS household keeps PUMA 0100100 and draws within it. + assert assigned.loc[100, "puma"] == "0100100" + assert assigned.loc[100, "congressional_district_geoid"] in {101, 102} + assert assigned.loc[100, "county_fips"] in {"01001", "01003"} + + +def test_asec_rows_draw_a_puma_within_state(tmp_path) -> None: + ladder = _ladder(tmp_path) + # No puma column at all — the pure-ASEC spine. + household = pd.DataFrame({"household_id": [1, 2], "state_fips": [1, 2]}) + + assigned = assign_us_puma_ladder(household, ladder, seed=0) + + assert set(assigned.loc[assigned["state_fips"] == 1, "puma"]).issubset( + {"0100100", "0100200"} + ) + assert assigned.loc[assigned["state_fips"] == 2, "puma"].tolist() == ["0200100"] + + +def test_mixed_acs_and_asec_frame(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame( + { + "household_id": [1, 2, 3], + "state_fips": [1, 1, 2], + # Row 0 is ACS (known PUMA); rows 1-2 are ASEC (NaN → drawn). + "puma": [100200, np.nan, np.nan], + } + ) + + assigned = assign_us_puma_ladder(household, ladder, seed=1) + + assert assigned.loc[0, "puma"] == "0100200" + assert assigned.loc[0, "congressional_district_geoid"] == 102 + assert assigned.loc[2, "puma"] == "0200100" + assert assigned.loc[1, "puma"] in {"0100100", "0100200"} + + +def test_assignment_is_deterministic_and_population_weighted(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame( + { + "household_id": range(600), + "state_fips": [1] * 600, + "puma": [100100] * 600, + } + ) + + first = assign_us_puma_ladder(household, ladder, seed=7) + second = assign_us_puma_ladder(household, ladder, seed=7) + + assert first["congressional_district_geoid"].tolist() == ( + second["congressional_district_geoid"].tolist() + ) + # CD 101 carries 900/1000 of PUMA 0100100's population. + share_101 = (first["congressional_district_geoid"] == 101).mean() + assert 0.85 < share_101 < 0.95 + # County 01001 carries the same 900/1000. + assert 0.85 < (first["county_fips"] == "01001").mean() < 0.95 + + +def test_state_to_puma_draw_matches_population_weights(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame({"household_id": range(6000), "state_fips": [1] * 6000}) + + assigned = assign_us_puma_ladder(household, ladder, seed=3) + + # PUMA 0100100 holds 1000/1500 of state 01's population. + share = (assigned["puma"] == "0100100").mean() + assert abs(share - (1000 / 1500)) < 0.03 + + +def test_assign_tract_derives_a_consistent_county(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame( + {"household_id": range(200), "state_fips": [1] * 200, "puma": [100100] * 200} + ) + + assigned = assign_us_puma_ladder(household, ladder, seed=0, assign_tract=True) + + assert "tract_geoid" in assigned.columns + tracts = assigned["tract_geoid"] + assert set(tracts).issubset({"01001000100", "01003000100"}) + # County derives structurally from the drawn tract — never an independent + # draw that could disagree. + assert (assigned["county_fips"] == tracts.str[:5]).all() + + +def test_congressional_district_is_stable_across_the_tract_flag(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame( + {"household_id": range(200), "state_fips": [1] * 200, "puma": [100100] * 200} + ) + + without_tract = assign_us_puma_ladder(household, ladder, seed=5) + with_tract = assign_us_puma_ladder(household, ladder, seed=5, assign_tract=True) + + assert without_tract["congressional_district_geoid"].tolist() == ( + with_tract["congressional_district_geoid"].tolist() + ) + + +def test_assignment_requires_state_fips(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame({"household_id": [1], "puma": [100100]}) + + with pytest.raises(ValueError, match="must contain 'state_fips'"): + assign_us_puma_ladder(household, ladder) + + +def test_assignment_refuses_puma_absent_from_ladder(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame({"household_id": [1], "state_fips": [1], "puma": [109999]}) + + with pytest.raises(ValueError, match="share a PUMA vintage"): + assign_us_puma_ladder(household, ladder) + + +def test_assignment_refuses_puma_state_disagreement(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame({"household_id": [1], "state_fips": [2], "puma": [100100]}) + + with pytest.raises(ValueError, match="disagrees with state_fips"): + assign_us_puma_ladder(household, ladder) + + +def test_assignment_refuses_state_without_pumas(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame({"household_id": [1], "state_fips": [5]}) + + with pytest.raises(ValueError, match="no PUMAs for state_fips 05"): + assign_us_puma_ladder(household, ladder) + + +def test_assignment_refuses_mismatched_cd_vintage(tmp_path) -> None: + ladder = _ladder(tmp_path) + household = pd.DataFrame({"household_id": [1], "state_fips": [1], "puma": [100100]}) + + with pytest.raises(ValueError, match="does not match the vintage"): + assign_us_puma_ladder( + household, + ladder, + expected_congressional_district_vintage="120th_congress", + ) + + +# --------------------------------------------------------------------------- # +# Frame integration, summary, and gate +# --------------------------------------------------------------------------- # + + +def test_with_household_preserves_frame_mass(tmp_path) -> None: + ladder = _ladder(tmp_path) + frame = _minimal_us_frame() + + assigned = with_household_us_puma_ladder( + frame, + ladder, + seed=0, + expected_congressional_district_vintage="119th_congress", + ) + + household = assigned.table("household") + assert household["puma"].tolist() == ["0100200", "0200100"] + assert assigned.weights_for("household").values.tolist() == [100.0, 300.0] + assert assigned.strata.tolist() == frame.strata.tolist() + + summary = us_puma_ladder_assignment_summary( + household, + ladder, + weight_values=assigned.weights_for("household").values, + ) + assert summary["applied"] is True + assert summary["ladder_pumas"] == 3 + assert summary["assigned_puma_values"] == 2 + assert summary["layer_vintages"]["congressional_district"] == "119th_congress" + + +def test_gate_passes_on_a_consistent_ladder() -> None: + household, weights = _gated_household() + + result = us_puma_ladder_gate(household, weights) + + assert result.passed, result.failures + assert 0.005 <= result.details["nyc_weighted_household_share"] <= 0.06 + + +def test_gate_fails_when_nyc_collapses_to_zero() -> None: + household, weights = _gated_household() + household.loc[household["county_fips"].isin(["36061"]), ["county_fips", "puma"]] = [ + "36001", + "3600100", + ] + + result = us_puma_ladder_gate(household, weights) + + assert not result.passed + assert any("in_nyc collapse" in failure for failure in result.failures) + + +def test_gate_fails_on_state_prefix_inconsistency() -> None: + household, weights = _gated_household() + household.loc[household.index[0], "county_fips"] = "06037" + + result = us_puma_ladder_gate(household, weights) + + assert not result.passed + assert any("disagree with the expected" in f for f in result.failures) + + +def test_gate_fails_when_columns_are_missing() -> None: + household, weights = _gated_household() + household = household.drop(columns=["county_fips"]) + + result = us_puma_ladder_gate(household, weights) + + assert not result.passed + assert any("missing geography column" in f for f in result.failures) + + +def _gated_household() -> tuple[pd.DataFrame, np.ndarray]: + """A weighted household table whose NYC mass sits inside the gate bounds. + + NYC is 2.6% of national weight and 41.9% of New York State's. + """ + + household = pd.DataFrame( + { + "household_id": [1, 2, 3, 4], + "state_fips": [36, 36, 6, 48], + "puma": ["3603801", "3600100", "0603701", "4800100"], + "congressional_district_geoid": [3612, 3620, 653, 4802], + "county_fips": ["36061", "36001", "06037", "48001"], + } + ) + weights = np.asarray([2.6, 3.6, 73.8, 20.0]) + return household, weights + + +def _minimal_us_frame() -> Frame: + tables = { + "person": pd.DataFrame( + { + "person_id": np.asarray([1, 2], dtype="int64"), + "person_household_id": np.asarray([1, 2], dtype="int64"), + "person_tax_unit_id": np.asarray([10, 20], dtype="int64"), + "person_spm_unit_id": np.asarray([100, 200], dtype="int64"), + "person_family_id": np.asarray([1000, 2000], dtype="int64"), + "person_marital_unit_id": np.asarray([10000, 20000], dtype="int64"), + } + ), + "household": pd.DataFrame( + { + "household_id": np.asarray([1, 2], dtype="int64"), + "state_fips": np.asarray([1, 2], dtype="int64"), + "puma": np.asarray([100200, 200100], dtype="int64"), + } + ), + "tax_unit": pd.DataFrame({"tax_unit_id": np.asarray([10, 20])}), + "spm_unit": pd.DataFrame({"spm_unit_id": np.asarray([100, 200])}), + "family": pd.DataFrame({"family_id": np.asarray([1000, 2000])}), + "marital_unit": pd.DataFrame({"marital_unit_id": np.asarray([10000, 20000])}), + } + weights = { + "household": Weights( + values=np.asarray([100.0, 300.0]), + kind=WeightKind.DESIGN, + ) + } + strata = pd.Series(["acs_2023", "acs_2023"], name="stratum") + return Frame(tables, US_SCHEMA, weights, strata) diff --git a/packages/populace-build/tests/test_us_puma_ladder_sources.py b/packages/populace-build/tests/test_us_puma_ladder_sources.py new file mode 100644 index 00000000..f464142d --- /dev/null +++ b/packages/populace-build/tests/test_us_puma_ladder_sources.py @@ -0,0 +1,191 @@ +"""Tests for the pure PUMA-ladder source parsers and assembler.""" + +import json + +import numpy as np +import pytest + +from populace.build.us_runtime.puma_ladder_sources import ( + assemble_us_puma_ladder, + parse_tract_to_puma_relationship, +) + +# Four populated blocks: PUMA 0100100 spans counties 01001/01003 and CDs +# 101/102; PUMA 0100200 sits in county 01003; PUMA 0200100 is the state-02 +# at-large district. Block ints drop the leading state zero, exactly as the +# P.L. 94-171 and CD BEF parsers return them. +_BLOCK_POPULATION = { + 10010001001000: 900, # tract 01001000100, county 01001 + 10030001001000: 100, # tract 01003000100, county 01003 + 10030002001000: 500, # tract 01003000200, county 01003 + 20130001001000: 400, # tract 02013000100, county 02013 +} +_CD_BY_BLOCK = { + 10010001001000: 101, + 10030001001000: 102, + 10030002001000: 102, + 20130001001000: 200, +} +_TRACT_TO_PUMA = { + 1001000100: 100100, + 1003000100: 100100, + 1003000200: 100200, + 2013000100: 200100, +} +_METADATA = { + "schema_version": 1, + "kind": "us_puma_ladder", + "puma_vintage": "2020_puma", + "sampling_basis": "population", + "layers": { + "congressional_district": {"vintage": "119th_congress", "source": "cd119"}, + "county": {"vintage": "2020_census", "source": "tract-to-puma"}, + "tract": {"vintage": "2020_census", "source": "tract-to-puma"}, + }, +} + + +def _relationship_lines(extra: list[str] | None = None) -> list[str]: + lines = [ + "STATEFP,COUNTYFP,TRACTCE,PUMA5CE", + "01,001,000100,00100", + "01,003,000100,00100", + "01,003,000200,00200", + "02,013,000100,00100", + "72,001,000100,00100", # Puerto Rico — filtered out. + ] + return lines + (extra or []) + + +def test_parse_tract_to_puma_filters_territories_and_builds_geoids() -> None: + mapping = parse_tract_to_puma_relationship( + _relationship_lines(), allowed_state_fips=frozenset({"01", "02"}) + ) + + assert mapping == _TRACT_TO_PUMA + # The Puerto Rico tract is excluded by the state filter. + assert 72001000100 not in mapping + + +def test_parse_tract_to_puma_keeps_all_states_without_a_filter() -> None: + mapping = parse_tract_to_puma_relationship(_relationship_lines()) + + assert 72001000100 in mapping + assert mapping[72001000100] == 7200100 + + +def test_parse_tract_to_puma_rejects_a_wrong_header() -> None: + with pytest.raises(ValueError, match="header must be"): + parse_tract_to_puma_relationship(["STATE,COUNTY,TRACT,PUMA", "01,001,1,1"]) + + +def test_parse_tract_to_puma_rejects_a_malformed_row() -> None: + with pytest.raises(ValueError, match="four fields"): + parse_tract_to_puma_relationship( + ["STATEFP,COUNTYFP,TRACTCE,PUMA5CE", "01,001,000100"] + ) + + +def test_parse_tract_to_puma_rejects_a_bad_width() -> None: + with pytest.raises(ValueError, match="5-digit"): + parse_tract_to_puma_relationship( + ["STATEFP,COUNTYFP,TRACTCE,PUMA5CE", "01,001,000100,100"] + ) + + +def test_parse_tract_to_puma_rejects_conflicting_pumas() -> None: + with pytest.raises(ValueError, match="both PUMA"): + parse_tract_to_puma_relationship( + [ + "STATEFP,COUNTYFP,TRACTCE,PUMA5CE", + "01,001,000100,00100", + "01,001,000100,00200", + ] + ) + + +def test_assemble_builds_conserving_overlap_tables() -> None: + payload = assemble_us_puma_ladder( + block_population=_BLOCK_POPULATION, + cd_by_block=_CD_BY_BLOCK, + tract_to_puma=_TRACT_TO_PUMA, + metadata=_METADATA, + ) + + assert payload["puma"].tolist() == [100100, 100200, 200100] + assert payload["puma_population"].tolist() == [1000, 500, 400] + + # CD overlap sorted by (puma, cd), conserving each PUMA's population. + assert list( + zip( + payload["cd_overlap_puma"].tolist(), + payload["cd_overlap_cd"].tolist(), + payload["cd_overlap_population"].tolist(), + strict=True, + ) + ) == [ + (100100, 101, 900), + (100100, 102, 100), + (100200, 102, 500), + (200100, 200, 400), + ] + # County overlap: PUMA 0100100 straddles counties 01001 and 01003. + assert list( + zip( + payload["county_overlap_puma"].tolist(), + payload["county_overlap_county"].tolist(), + payload["county_overlap_population"].tolist(), + strict=True, + ) + ) == [ + (100100, 1001, 900), + (100100, 1003, 100), + (100200, 1003, 500), + (200100, 2013, 400), + ] + assert payload["tract_overlap_tract"].tolist() == [ + 1001000100, + 1003000100, + 1003000200, + 2013000100, + ] + metadata = json.loads(str(payload["metadata_json"])) + assert metadata["kind"] == "us_puma_ladder" + + +def test_assemble_refuses_a_block_with_no_puma() -> None: + with pytest.raises(ValueError, match="tract absent from the tract-to-PUMA"): + assemble_us_puma_ladder( + block_population={**_BLOCK_POPULATION, 30070001001000: 50}, + cd_by_block={**_CD_BY_BLOCK, 30070001001000: 301}, + tract_to_puma=_TRACT_TO_PUMA, # no tract 3007000100 + metadata=_METADATA, + ) + + +def test_assemble_refuses_a_block_with_no_congressional_district() -> None: + with pytest.raises(ValueError, match="no congressional district"): + assemble_us_puma_ladder( + block_population=_BLOCK_POPULATION, + cd_by_block={ + key: value + for key, value in _CD_BY_BLOCK.items() + if key != 20130001001000 + }, + tract_to_puma=_TRACT_TO_PUMA, + metadata=_METADATA, + ) + + +def test_assemble_ignores_zero_population_blocks() -> None: + payload = assemble_us_puma_ladder( + block_population={**_BLOCK_POPULATION, 10010001009000: 0}, + cd_by_block=_CD_BY_BLOCK, + tract_to_puma=_TRACT_TO_PUMA, + metadata=_METADATA, + ) + + # The zero-population block never needs a CD or tract lookup and does not + # change any PUMA's conserved total. + assert payload["puma_population"].tolist() == [1000, 500, 400] + assert np.asarray(payload["cd_overlap_population"]).sum() == 1900 diff --git a/tools/build_us_puma_ladder_artifact.py b/tools/build_us_puma_ladder_artifact.py new file mode 100644 index 00000000..2d86eb72 --- /dev/null +++ b/tools/build_us_puma_ladder_artifact.py @@ -0,0 +1,257 @@ +"""Build the US PUMA-ladder artifact from primary Census sources. + +Downloads (with a local cache shared with the block-ladder builder) the 2020 +Census Tract to 2020 PUMA relationship file, the 119th Congressional District +block equivalency file, and the 2020 P.L. 94-171 geographic headers (block +populations); joins them at 2020-tabulation-block grain and aggregates to the +PUMA anchor plus three ``(PUMA, layer_value, population)`` overlap tables; and +writes one national NPZ artifact whose embedded metadata records a vintage and +source per derived layer (``vintage_policy: error`` — the loader refuses an +artifact missing any of them, or one whose overlaps do not conserve each PUMA's +population). No per-area files, the standing rule. + +Because 2020 blocks nest in 2020 tracts nest in 2020 PUMAs, the only source +this builder adds over the block ladder is the small tract-to-PUMA relationship +file; block populations and 119th-CD assignments reuse the block-ladder +parsers unchanged. + +The artifact is self-checked by loading it back through +``populace.build.us_runtime.load_us_puma_ladder`` before the summary is +written, so a published ladder is by construction a loadable ladder. Every +source's URL + SHA-256, the output SHA-256, and the national conservation +totals are recorded in the summary JSON — the pinned-source manifest. + +Example: + uv run python tools/build_us_puma_ladder_artifact.py \ + --out build/us/us_puma_ladder_2020.npz \ + --cache-dir ~/.cache/populace-us-geography + + # Smoke run over two small states: + uv run python tools/build_us_puma_ladder_artifact.py \ + --out /tmp/puma_ladder_smoke.npz --states 10,50 +""" + +from __future__ import annotations + +import argparse +import hashlib +import io +import json +import sys +import urllib.request +import zipfile +from pathlib import Path + +import numpy as np + +from populace.build.us_runtime import load_us_puma_ladder +from populace.build.us_runtime.block_ladder_sources import ( + US_STATES, + parse_national_cd_bef, + parse_pl_geo_blocks, +) +from populace.build.us_runtime.puma_ladder_sources import ( + assemble_us_puma_ladder, + parse_tract_to_puma_relationship, +) + +TRACT_TO_PUMA_URL = ( + "https://www2.census.gov/geo/docs/maps-data/data/rel2020/" + "2020_Census_Tract_to_2020_PUMA.txt" +) +CD119_BEF_URL = ( + "https://www2.census.gov/programs-surveys/decennial/rdo/mapping-files/" + "2025/119-congressional-district-befs/cd119.zip" +) +CD119_NATIONAL_MEMBER = "NationalCD119.txt" +PL94171_URL_TEMPLATE = ( + "https://www2.census.gov/programs-surveys/decennial/2020/data/" + "01-Redistricting_File--PL_94-171/{dirname}/{usps_lower}2020.pl.zip" +) + +LAYER_VINTAGES = { + "congressional_district": { + "vintage": "119th_congress", + "source": "Census 119th Congressional District BEF (NationalCD119.txt)", + "url": CD119_BEF_URL, + }, + "county": { + "vintage": "2020_census", + "source": ( + "Census 2020 Census Tract to 2020 PUMA relationship file " + "(county is the tract geoid's structural prefix)" + ), + "url": TRACT_TO_PUMA_URL, + }, + "tract": { + "vintage": "2020_census", + "source": ( + "Census 2020 Census Tract to 2020 PUMA relationship file, weighted " + "by 2020 P.L. 94-171 block populations" + ), + "url": TRACT_TO_PUMA_URL, + }, +} + + +def _parse_args(argv: list[str] | None = None) -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Build the national US PUMA-ladder NPZ artifact." + ) + parser.add_argument("--out", required=True, type=Path) + parser.add_argument( + "--cache-dir", + type=Path, + default=Path.home() / ".cache" / "populace-us-geography", + help="Download cache; re-runs reuse verified files.", + ) + parser.add_argument( + "--states", + help=( + "Optional comma-separated state FIPS subset (smoke runs only; a " + "published ladder covers all 51)." + ), + ) + parser.add_argument( + "--summary-json", + type=Path, + help="Path for the build summary. Defaults beside --out.", + ) + return parser.parse_args(argv) + + +def _log(message: str) -> None: + print(message, file=sys.stderr, flush=True) + + +def _sha256(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as stream: + for chunk in iter(lambda: stream.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def _download(url: str, cache_dir: Path) -> Path: + cache_dir.mkdir(parents=True, exist_ok=True) + destination = cache_dir / url.rsplit("/", 1)[-1] + if destination.exists() and destination.stat().st_size > 0: + return destination + _log(f" downloading {url}") + request = urllib.request.Request( + url, headers={"User-Agent": "populace-build (puma ladder artifact)"} + ) + with urllib.request.urlopen(request) as response: + payload = response.read() + if not payload: + raise RuntimeError(f"Empty download from {url}") + partial = destination.with_suffix(destination.suffix + ".partial") + partial.write_bytes(payload) + partial.replace(destination) + return destination + + +def _zip_member_lines(archive_path: Path, member: str) -> list[str]: + with zipfile.ZipFile(archive_path) as archive: + with archive.open(member) as stream: + return io.TextIOWrapper(stream, encoding="latin-1").readlines() + + +def _text_lines(path: Path) -> list[str]: + with path.open(encoding="utf-8-sig") as stream: + return stream.readlines() + + +def main(argv: list[str] | None = None) -> None: + args = _parse_args(argv) + selected = ( + {value.strip() for value in args.states.split(",")} if args.states else None + ) + states = [entry for entry in US_STATES if selected is None or entry[0] in selected] + if selected is not None: + unknown = selected - {entry[0] for entry in US_STATES} + if unknown: + raise SystemExit(f"Unknown state FIPS in --states: {sorted(unknown)}") + allowed_state_fips = frozenset(fips for fips, _, _ in states) + _log(f"Building US PUMA ladder for {len(states)} state(s)") + + source_files: dict[str, dict[str, str]] = {} + + tract_puma_txt = _download(TRACT_TO_PUMA_URL, args.cache_dir) + source_files["tract_to_puma"] = { + "url": TRACT_TO_PUMA_URL, + "sha256": _sha256(tract_puma_txt), + } + _log(" parsing tract-to-PUMA relationship") + tract_to_puma = parse_tract_to_puma_relationship( + _text_lines(tract_puma_txt), allowed_state_fips=allowed_state_fips + ) + _log(f" tract-to-PUMA: {len(tract_to_puma):,} tracts") + + cd_zip = _download(CD119_BEF_URL, args.cache_dir) + source_files["cd119_bef"] = {"url": CD119_BEF_URL, "sha256": _sha256(cd_zip)} + _log(" parsing national CD119 BEF") + cd_by_block = parse_national_cd_bef( + _zip_member_lines(cd_zip, CD119_NATIONAL_MEMBER) + ) + + block_population: dict[int, int] = {} + for fips, usps, dirname in states: + _log(f" state {fips} {usps}") + pl_url = PL94171_URL_TEMPLATE.format(dirname=dirname, usps_lower=usps.lower()) + pl_zip = _download(pl_url, args.cache_dir) + source_files[f"pl94171_{usps.lower()}"] = { + "url": pl_url, + "sha256": _sha256(pl_zip), + } + geo_member = f"{usps.lower()}geo2020.pl" + state_blocks = parse_pl_geo_blocks( + _zip_member_lines(pl_zip, geo_member), state_fips=fips + ) + block_population.update(state_blocks) + + metadata = { + "schema_version": 1, + "kind": "us_puma_ladder", + "puma_vintage": "2020_puma", + "sampling_basis": "population", + "layers": LAYER_VINTAGES, + "states": [fips for fips, _, _ in states], + "source_files": source_files, + } + _log(f" aggregating {len(block_population):,} populated blocks") + payload = assemble_us_puma_ladder( + block_population=block_population, + cd_by_block=cd_by_block, + tract_to_puma=tract_to_puma, + metadata=metadata, + ) + args.out.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed(args.out, **payload) + + ladder = load_us_puma_ladder(args.out) # self-check: must load cleanly + summary = { + "artifact": args.out.name, + "output_sha256": _sha256(args.out), + "pumas": int(len(ladder)), + "population": int(ladder.puma_population.sum()), + "congressional_district_overlaps": int(len(ladder.cd_overlap_puma)), + "county_overlaps": int(len(ladder.county_overlap_puma)), + "tract_overlaps": int(len(ladder.tract_overlap_puma)), + "congressional_districts": int(np.unique(ladder.cd_overlap_cd).size), + "counties": int(np.unique(ladder.county_overlap_county).size), + "layer_vintages": ladder.layer_vintages, + "states": [fips for fips, _, _ in states], + "source_files": source_files, + } + summary_path = ( + args.summary_json + if args.summary_json is not None + else args.out.with_suffix(".summary.json") + ) + summary_path.write_text(json.dumps(summary, indent=2, sort_keys=True) + "\n") + print(json.dumps(summary, indent=2, sort_keys=True)) + + +if __name__ == "__main__": + main()