From 35f96fe5b9bd2b78848fe7b2c74b74a2925a7f1d Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Wed, 7 Oct 2026 08:17:27 -0400 Subject: [PATCH 1/6] Add donor transport operators and weighted quantile mapping Authenticate local populace-US donor inputs before decoding concepts; recode household roster pointers, preserve donor weights and support provenance, and derive stable take-up seeds from exact scoped source identities. Add one-step currency conversion and signed weighted quantile maps with declared interpolation, component isolation, and pro-rata entity aggregation. Cover the invariants with synthetic round trips, independent differentials, Hypothesis properties, and explicit float64 precision refusals. Co-Authored-By: GPT-6.1 Sol --- changelog.d/nz-transport-operators.added.md | 3 + docs/concept-schema-transport-adr.md | 14 +- .../src/microcosm/frame/transport.py | 657 ++++++++++++++++++ .../shared/test_transport_donor.py | 558 +++++++++++++++ .../shared/test_transport_quantile.py | 484 +++++++++++++ .../fixtures/transport/synthetic_donor.json | 236 +++++++ 6 files changed, 1947 insertions(+), 5 deletions(-) create mode 100644 changelog.d/nz-transport-operators.added.md create mode 100644 packages/microcosm-frame/src/microcosm/frame/transport.py create mode 100644 packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py create mode 100644 packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py create mode 100644 packages/microcosm-frame/tests/fixtures/transport/synthetic_donor.json diff --git a/changelog.d/nz-transport-operators.added.md b/changelog.d/nz-transport-operators.added.md new file mode 100644 index 0000000000..75ad7b707f --- /dev/null +++ b/changelog.d/nz-transport-operators.added.md @@ -0,0 +1,3 @@ +Add engine-free transport operators for authenticated local US donor concepts, +stable source-ID seeds, one-step currency conversion, and weighted signed +quantile maps with component and entity-aggregate support. diff --git a/docs/concept-schema-transport-adr.md b/docs/concept-schema-transport-adr.md index 3a34d1b265..13bf329550 100644 --- a/docs/concept-schema-transport-adr.md +++ b/docs/concept-schema-transport-adr.md @@ -1,9 +1,10 @@ # Engine-neutral concepts and transport (ADR) **Status:** accepted for the concept schema, the per-engine mappings and their -coverage reports, which this change adds. The transport pipeline described -below is the design those pieces serve; none of it is built yet, and no build -has migrated. **Date:** 2026-09-28. **Source:** Max's 27 September rulings on +coverage reports. The transport operators now include a pinned local US donor +reader, source-ID seeds, a currency bridge and weighted quantile mapping; +the composed transport pipeline below remains a design, and no build has +migrated. **Date:** 2026-09-28. **Source:** Max's 27 September rulings on international populations and law-anchored concepts. **Amended 2026-10-07:** schema version 2 adds `liquid_financial_assets`, the first stock concept, for New Zealand's Accommodation Supplement cash-asset test. Only the Axiom New @@ -284,8 +285,11 @@ These hold for every input and are tested (Hypothesis properties unless noted): encoding it through the adapter's mapping, one build at a time. - The spec engine's `imputation.yaml` `concepts:` block is unrelated. It groups predictor columns for imputation, and its names are not concept ids. -- Transport still needs four pieces beyond the schema: the donor bank, the - income transformation (step 3), a unit-construction step that builds engine +- `microcosm.frame.transport` supplies the donor bank and the income + transformation operators (step 3). Its bank declares US donor provenance + and source-currency amounts; it does not declare destination residents. + The consuming builder supplies the destination frame declaration after + transport. Transport still needs a unit-construction step that builds engine units from the concept pointers (until it exists, 13 of New Zealand's 80 and 46 of Belgium's 105 module bindings, and the take-up and housing bindings on US and UK group entities, stay deferred), and a country pack with geography, diff --git a/packages/microcosm-frame/src/microcosm/frame/transport.py b/packages/microcosm-frame/src/microcosm/frame/transport.py new file mode 100644 index 0000000000..e37888a7cc --- /dev/null +++ b/packages/microcosm-frame/src/microcosm/frame/transport.py @@ -0,0 +1,657 @@ +"""Country-neutral operators for transporting a declared donor population.""" + +import json +from collections.abc import Mapping, Sequence +from dataclasses import dataclass +from hashlib import blake2b, file_digest +from math import fsum +from pathlib import Path +from types import MappingProxyType + +import numpy as np +import pandas as pd +from numpy.typing import ArrayLike + +from microcosm.frame.adapters.policyengine_us_concepts import ( + POLICYENGINE_US_CONCEPT_MAPPING, +) +from microcosm.frame.concepts import ( + HOUSEHOLD_ID_COLUMN, + PERSON_HOUSEHOLD_ID_COLUMN, + PERSON_ID_COLUMN, + ContentBasis, + Unit, + concept_for_column, + validate_concept_tables, +) +from microcosm.frame.materialize import read_frame_table +from microcosm.frame.scaling import apply_scale + +__all__ = [ + "DonorBank", + "currency_bridge", + "derive_transport_seed", + "quantile_map", + "read_populace_us_donor", +] + + +@dataclass(frozen=True) +class _QuantileBand: + lower: float + upper: float | None + share: float + + +def _quantile_bands(rows: Sequence[Mapping[str, object]]) -> tuple[_QuantileBand, ...]: + if not rows: + raise ValueError("Target bands must be nonempty.") + bands = [] + for row in rows: + if not isinstance(row, Mapping) or set(row) != {"lower", "upper", "share"}: + raise ValueError( + "Each target band requires exactly lower, upper and share." + ) + lower = float(row["lower"]) + upper = None if row["upper"] is None else float(row["upper"]) + share = float(row["share"]) + if not np.isfinite(lower) or lower < 0: + raise ValueError("Band lower bounds must be finite and nonnegative.") + if upper is not None and (not np.isfinite(upper) or upper <= lower): + raise ValueError( + "Band upper bounds must be finite and exceed lower bounds." + ) + if not np.isfinite(share) or share <= 0: + raise ValueError("Band shares must be finite and positive.") + if bands and bands[-1].upper != lower: + raise ValueError("Target bands must be ordered and contiguous.") + bands.append(_QuantileBand(lower, upper, share)) + if not np.isclose(fsum(band.share for band in bands), 1.0, rtol=0, atol=1e-12): + raise ValueError("Target band shares must sum to one.") + return tuple(bands) + + +def _quantile_interpolation( + interpolation: str | Mapping[str, object], +) -> tuple[str, float | None]: + if interpolation == "uniform": + return "uniform", None + if not isinstance(interpolation, Mapping) or set(interpolation) != { + "method", + "alpha", + }: + raise ValueError( + "Interpolation must be uniform or {method: pareto, alpha: ...}." + ) + if interpolation["method"] != "pareto": + raise ValueError("The declared interpolation method must be pareto.") + alpha = float(interpolation["alpha"]) + if not np.isfinite(alpha) or alpha <= 0: + raise ValueError("Pareto alpha must be finite and positive.") + return "pareto", alpha + + +def _quantile_label_codes(labels: Sequence[object], size: int, name: str) -> np.ndarray: + if isinstance(labels, (str, bytes)): + raise ValueError(f"{name} must be a sequence of labels.") + labels = pd.Series(list(labels), dtype=object) + if len(labels) != size or labels.isna().any(): + raise ValueError(f"{name} must have one nonmissing label per value.") + try: + codes, _ = pd.factorize(labels, sort=False) + except TypeError as error: + raise ValueError(f"{name} labels must be hashable.") from error + return codes + + +def _quantile_component( + magnitudes: np.ndarray, + weights: np.ndarray, + bands: tuple[_QuantileBand, ...], + method: str, + alpha: float | None, +) -> np.ndarray: + # Pool equal source values: splitting tied records must not change their rank. + _, inverse = np.unique(magnitudes, return_inverse=True) + scaled_weights = weights / weights.max() + if (scaled_weights <= 0).any(): + raise ValueError("Weight ratios are too large for finite weighted ranks.") + masses = np.bincount(inverse, weights=scaled_weights) + total = fsum(masses) + ranks = (np.cumsum(masses) - masses / 2) / total + ranks = np.clip(ranks, np.nextafter(0.0, 1.0), np.nextafter(1.0, 0.0)) + shares = np.array([band.share for band in bands]) + boundaries = np.cumsum(shares / fsum(shares)) + boundaries[-1] = 1.0 + positions = np.searchsorted(boundaries, ranks, side="right") + mapped = np.empty_like(ranks) + for position, band in enumerate(bands): + selected = positions == position + start = 0.0 if position == 0 else boundaries[position - 1] + fraction = (ranks[selected] - start) / (boundaries[position] - start) + if method == "pareto" and position == len(bands) - 1: + assert alpha is not None + # A finite upper bound declares a truncated Pareto distribution. + tail = 0.0 if band.upper is None else (band.lower / band.upper) ** alpha + with np.errstate(over="ignore", invalid="ignore", divide="ignore"): + mapped[selected] = band.lower * np.exp( + -np.log1p(-fraction * (1 - tail)) / alpha + ) + else: + assert band.upper is not None + mapped[selected] = band.lower + fraction * (band.upper - band.lower) + if not np.isfinite(mapped).all() or (mapped <= 0).any(): + raise ValueError("Mapped nonzero values must have finite positive magnitudes.") + if (np.diff(mapped) <= 0).any(): + raise ValueError( + "Mapped distinct ranks must remain distinct at float64 precision." + ) + return mapped[inverse] + + +def _quantile_map_components( + values: np.ndarray, + weights: np.ndarray, + components: np.ndarray, + sign_bands: Mapping[str, tuple[_QuantileBand, ...]], + method: str, + alpha: float | None, +) -> np.ndarray: + result = values.copy() + order = np.argsort(components, kind="stable") + blocks = np.split(order, np.flatnonzero(np.diff(components[order])) + 1) + for block in blocks: + for sign, name in ((1, "positive"), (-1, "negative")): + selected = block[values[block] * sign > 0] + if len(selected): + result[selected] = sign * _quantile_component( + np.abs(values[selected]), + weights[selected], + sign_bands[name], + method, + alpha, + ) + return result + + +def quantile_map( + values: ArrayLike, + weights: ArrayLike, + target_bands: Sequence[Mapping[str, object]] + | Mapping[str, Sequence[Mapping[str, object]]], + *, + interpolation: str | Mapping[str, object], + group: Sequence[object] | Mapping[str, Sequence[object]] | None = None, +) -> np.ndarray: + """Map weighted midpoint ranks to declared bands, returning a new float64 array. + + ``target_bands`` is a sequence of ``{lower, upper, share}`` rows over + positive *magnitudes*, or a mapping with exactly ``positive`` and + ``negative`` sequences. Bands are contiguous and increasing; shares sum + to one separately for each sign, conditional on being nonzero. Zero + values retain their value and do not enter either rank distribution. + Inputs must be finite one-dimensional arrays with strictly positive + weights. Weights are never changed. + + ``interpolation="uniform"`` requires finite upper bounds. A declaration + ``{"method": "pareto", "alpha": ...}`` applies Pareto interpolation to + the top band only, with a positive lower bound. Its upper bound may be + ``None`` (unbounded) or finite (truncated); lower bands remain uniform. + + A sequence passed as ``group`` gives one component label per value. + Ranks are computed independently within each component and sign. Equal + magnitudes share their pooled weighted midpoint rank. Consequently the + error in a band's weighted share is bounded by the largest tied block's + fraction of the component/sign weight: a tied block is indivisible. + + Alternatively ``group={"entity_ids": [...], "components": [...]}`` + declares member-to-entity memberships. ``components`` is optional. + Members of an entity must repeat its weight and component label. Their + values are summed, entity totals are mapped using each entity weight + once, and members are rescaled pro rata by a positive factor. A floating + rounding residual is assigned to the largest nonzero member so that + tied entity totals remain exactly tied on repeated mapping. This also + preserves member signs for mixed-sign assets. A zero-total entity is + unchanged, including cancelling nonzero members. Component labels and + entity ids must be nonmissing hashable scalars. Mapping fails if float64 + cannot retain distinct ranks, signs, pro rata precision, or the mapped + aggregate totals (including severe mixed-sign cancellation). + """ + values = np.asarray(values, dtype=np.float64) + weights = np.asarray(weights, dtype=np.float64) + if values.ndim != 1 or weights.shape != values.shape: + raise ValueError( + "Values and weights must have matching one-dimensional shapes." + ) + if not np.isfinite(values).all(): + raise ValueError("Values must be finite.") + if not np.isfinite(weights).all() or (weights <= 0).any(): + raise ValueError("Weights must be finite and positive.") + method, alpha = _quantile_interpolation(interpolation) + if isinstance(target_bands, Mapping): + if set(target_bands) != {"positive", "negative"}: + raise ValueError( + "Signed target bands require exactly positive and negative." + ) + sign_bands = { + sign: _quantile_bands(rows) for sign, rows in target_bands.items() + } + else: + bands = _quantile_bands(target_bands) + sign_bands = {"positive": bands, "negative": bands} + for bands in sign_bands.values(): + if method == "uniform" and bands[-1].upper is None: + raise ValueError("Uniform interpolation requires finite upper bounds.") + if method == "pareto" and bands[-1].lower <= 0: + raise ValueError("The Pareto top band requires a positive lower bound.") + if not isinstance(group, Mapping): + components = ( + np.zeros(len(values), dtype=np.int64) + if group is None + else _quantile_label_codes(group, len(values), "Components") + ) + return _quantile_map_components( + values, weights, components, sign_bands, method, alpha + ) + if "entity_ids" not in group or not set(group) <= {"entity_ids", "components"}: + raise ValueError("Aggregate groups require entity_ids and optional components.") + entities = _quantile_label_codes(group["entity_ids"], len(values), "Entity ids") + components = ( + np.zeros(len(values), dtype=np.int64) + if "components" not in group + else _quantile_label_codes(group["components"], len(values), "Components") + ) + count = int(entities.max()) + 1 if len(entities) else 0 + totals = np.empty(count) + entity_weights = np.empty(count) + entity_components = np.empty(count, dtype=np.int64) + order = np.argsort(entities, kind="stable") + members = ( + np.split(order, np.flatnonzero(np.diff(entities[order])) + 1) if count else [] + ) + for entity, selected in enumerate(members): + if (weights[selected] != weights[selected[0]]).any(): + raise ValueError("Members must repeat one entity weight.") + if (components[selected] != components[selected[0]]).any(): + raise ValueError("Members must repeat one entity component label.") + try: + totals[entity] = fsum(values[selected]) + except OverflowError as error: + raise ValueError("Entity totals must be finite.") from error + entity_weights[entity] = weights[selected[0]] + entity_components[entity] = components[selected[0]] + mapped = _quantile_map_components( + totals, entity_weights, entity_components, sign_bands, method, alpha + ) + result = values.copy() + with np.errstate(over="ignore", invalid="ignore", divide="ignore"): + for entity, selected in enumerate(members): + if totals[entity] != 0: + result[selected] = (values[selected] / totals[entity]) * mapped[entity] + if not np.isfinite(result[selected]).all(): + raise ValueError( + "Pro rata mapping must preserve finite member values." + ) + nonzero = selected[values[selected] != 0] + # Put the sum-rounding residual on the largest member so an + # otherwise tiny holding is never materially changed. Severe + # cancellation may make the exact total unrepresentable. + anchor = nonzero[np.argmax(np.abs(result[nonzero]))] + original_anchor = result[anchor] + for attempt in range(3): + current_total = fsum(result[selected]) + if current_total == mapped[entity]: + break + if attempt == 0: + result[anchor] += mapped[entity] - current_total + else: + direction = ( + np.inf if current_total < mapped[entity] else -np.inf + ) + result[anchor] = np.nextafter(result[anchor], direction) + if fsum(result[selected]) != mapped[entity]: + raise ValueError( + "Pro rata mapping cannot retain exact totals at float64 precision." + ) + correction_bound = 4 * len(selected) * abs(np.spacing(original_anchor)) + if abs(result[anchor] - original_anchor) > correction_bound: + raise ValueError( + "Pro rata correction exceeds float64 rounding precision." + ) + if not np.isfinite(result).all() or (np.sign(result) != np.sign(values)).any(): + raise ValueError( + "Pro rata mapping must preserve finite member values and signs." + ) + return result + + +@dataclass(frozen=True) +class DonorBank: + """US donor concepts and aligned provenance, never destination microdata. + + ``weights`` and ``support_strata`` follow the household table's row order; + ``source_person_ids`` follows the person table. Weights are the untouched + donor design weights. Explicit ``person_source_id`` values retain their + exact integer/string identities, including repetitions on support clones. + Otherwise full source lineage is encoded as a JSON string of the support + stratum and typed year/household/person IDs, avoiding household-local ID + collisions; files without lineage use their pinned donor ``person_id``. + ``currency`` describes source + amounts; the destination producer supplies its own frame declaration after + transport. No destination country or year is inferred by this reader. + """ + + tables: Mapping[str, pd.DataFrame] + weights: np.ndarray + support_strata: np.ndarray + source_person_ids: np.ndarray + content_basis: ContentBasis = ContentBasis.TRANSPORT + donor_country: str = "us" + currency: str = "USD" + + +def _line_numbers(values: pd.Series, *, nullable: bool = False) -> np.ndarray: + numbers = pd.to_numeric(values, errors="raise").to_numpy( + dtype=np.float64, na_value=np.nan + ) + if nullable: + numbers = np.where(np.isnan(numbers), 0, numbers) + if not np.isfinite(numbers).all() or np.any(numbers != np.floor(numbers)): + raise ValueError(f"Donor {values.name} must contain integer line numbers.") + if np.any(numbers < (0 if nullable else 1)) or np.any(numbers >= 2**63): + raise ValueError(f"Donor {values.name} contains invalid line numbers.") + return numbers.astype(np.int64) + + +def _donor_extras(person: pd.DataFrame, household: pd.DataFrame) -> dict: + """Translate donor-only columns before crossing the concept decode seam.""" + required = ( + "A_LINENO", + "A_SPOUSE", + "PEPAR1", + "PEPAR2", + "A_EXPRRP", + "bank_account_assets", + "stock_assets", + "bond_assets", + ) + missing = set(required) - set(person.columns) + if missing: + raise ValueError(f"Donor person table lacks columns {sorted(missing)}.") + if "household_weight" not in household: + raise ValueError("Donor household table lacks household_weight.") + # Validate the structural IDs without turning large integers into floats. + validate_concept_tables( + { + "person": person[[PERSON_ID_COLUMN, PERSON_HOUSEHOLD_ID_COLUMN]], + "household": household[[HOUSEHOLD_ID_COLUMN]], + } + ) + line = _line_numbers(person["A_LINENO"]) + household_ids = person[PERSON_HOUSEHOLD_ID_COLUMN].to_numpy() + roster = pd.MultiIndex.from_arrays([household_ids, line]) + if not roster.is_unique: + raise ValueError("Donor line numbers must be unique within a household.") + if (person[PERSON_ID_COLUMN] > np.iinfo(np.int64).max).any(): + raise ValueError("Donor person IDs must fit the concept pointer's int64 dtype.") + person_ids = person[PERSON_ID_COLUMN].to_numpy(dtype=np.int64) + pointers = {} + for source, destination in ( + ("A_SPOUSE", "partner_person_id"), + ("PEPAR1", "parent_1_person_id"), + ("PEPAR2", "parent_2_person_id"), + ): + pointed_line = _line_numbers(person[source], nullable=True) + present = pointed_line != 0 + rows = roster.get_indexer( + pd.MultiIndex.from_arrays([household_ids, pointed_line]) + ) + if np.any(present & (rows < 0)): + raise ValueError(f"Donor {source} points to an unknown household line.") + pointer = pd.array([pd.NA] * len(person), dtype="Int64") + pointer[present] = person_ids[rows[present]] + pointers[destination] = pointer + # A sole co-resident parent is always parent 1 in the concept contract. + move = ( + pointers["parent_1_person_id"].isna() & ~pointers["parent_2_person_id"].isna() + ) + pointers["parent_1_person_id"][move] = pointers["parent_2_person_id"][move] + pointers["parent_2_person_id"][move] = pd.NA + relationship = _line_numbers(person["A_EXPRRP"]) + reference = np.isin(relationship, (1, 2)) + reference_households = pd.Index(household_ids[reference]) + if not reference_households.is_unique: + raise ValueError("Donor household has multiple reference persons.") + reference_rows = reference_households.get_indexer(household[HOUSEHOLD_ID_COLUMN]) + if np.any(reference_rows < 0): + raise ValueError("Donor household has no reference person.") + reference_ids = person_ids[reference][reference_rows] + # The CPS relationship recode identifies the reference person's unmarried + # partner as 13 (the ACS-to-CPS crosswalk uses the same code). A_SPOUSE + # alone does not capture these cohabiting couples. + unmarried_rows = np.flatnonzero(relationship == 13) + if not pd.Index(household_ids[unmarried_rows]).is_unique: + raise ValueError("Donor household has multiple unmarried partners.") + reference_positions = np.flatnonzero(reference) + for partner_row in unmarried_rows: + reference_row = reference_positions[ + reference_households.get_indexer([household_ids[partner_row]])[0] + ] + partners = pointers["partner_person_id"] + for own_row, other_row in ( + (partner_row, reference_row), + (reference_row, partner_row), + ): + if ( + not pd.isna(partners[own_row]) + and partners[own_row] != person_ids[other_row] + ): + raise ValueError( + "Donor unmarried-partner roster conflicts with spouse pointers." + ) + partners[partner_row] = person_ids[reference_row] + partners[reference_row] = person_ids[partner_row] + # The donor's three leaves are not invertible through an encoding mapping: + # their sum is the reviewed donor-side stock, without inventing a split. + leaves = person[["bank_account_assets", "stock_assets", "bond_assets"]].to_numpy( + dtype=np.float64 + ) + if not np.isfinite(leaves).all() or np.any(leaves < 0): + raise ValueError("Donor liquid-asset leaves must be finite and nonnegative.") + assets = leaves.sum(axis=1) + weights = household["household_weight"].to_numpy(dtype=np.float64, copy=True) + if not np.isfinite(weights).all() or np.any(weights < 0) or not np.any(weights > 0): + raise ValueError( + "Donor weights must be finite, nonnegative, with positive mass." + ) + # Older Build P files have one undifferentiated US support stratum. Never + # infer a survey channel from income or engine inputs when labels are absent. + channel_column = "household_support_channel" + if channel_column in household: + channels = household[channel_column].to_numpy(dtype=object, copy=True) + if any(not isinstance(value, str) or not value for value in channels): + raise ValueError( + "Donor household support channels must be nonempty strings." + ) + strata = np.asarray([f"us:{value}" for value in channels], dtype=object) + else: + channels = None + strata = np.full(len(household), "us", dtype=object) + if "person_support_channel" in person: + if channels is None: + raise ValueError( + "Person support channels require household support channels." + ) + hh_rows = pd.Index(household[HOUSEHOLD_ID_COLUMN]).get_indexer(household_ids) + if not np.array_equal( + person["person_support_channel"].to_numpy(), channels[hh_rows] + ): + raise ValueError("Donor person and household support channels disagree.") + lineage = ("source_year", "source_household_id", "source_person_id") + if "person_source_id" in person: + source_ids = person["person_source_id"].to_numpy(copy=True) + elif any(name in person for name in lineage): + if not all(name in person for name in lineage): + raise ValueError( + "Donor source lineage requires year, household and person IDs." + ) + parts = [_source_id_payloads(person[name]) for name in lineage] + hh_rows = pd.Index(household[HOUSEHOLD_ID_COLUMN]).get_indexer(household_ids) + source_ids = np.asarray( + [ + json.dumps( + [strata[hh_rows[index]], *items], + ensure_ascii=False, + separators=(",", ":"), + ) + for index, items in enumerate(zip(*parts, strict=True)) + ], + dtype=object, + ) + else: + source_ids = person[PERSON_ID_COLUMN].to_numpy(copy=True) + # Check exact identities with the same contract used by seed derivation. + _source_id_payloads(source_ids) + return { + "pointers": pointers, + "reference_person_id": pd.array(reference_ids, dtype="Int64"), + "liquid_financial_assets": assets, + "weights": weights, + "support_strata": strata, + "source_person_ids": source_ids, + } + + +def read_populace_us_donor(path: str | Path, *, sha256: str, size: int) -> DonorBank: + """Authenticate a local donor H5, then decode its person/household inputs. + + Size and SHA-256 are mandatory and verified before any HDF parsing. There + is no download or alternate donor fallback. CPS pointers are resolved by + (household id, line number), and invalid/asymmetric relations fail concept + validation. Donor-only enrichment is staged before ``decode``; all work + after that seam reads concept names. The result declares transport content + with a US donor and unchanged source-currency amounts and design weights. + """ + path = Path(path) + if isinstance(size, bool) or not isinstance(size, int) or size < 0: + raise ValueError("Donor size must be a nonnegative integer.") + if ( + not isinstance(sha256, str) + or len(sha256) != 64 + or any(c not in "0123456789abcdef" for c in sha256) + ): + raise ValueError("Donor SHA-256 must be 64 lowercase hexadecimal characters.") + with path.open("rb") as source: + if path.stat().st_size != size: + raise ValueError("Donor size mismatch.") + if file_digest(source, "sha256").hexdigest() != sha256: + raise ValueError("Donor SHA-256 mismatch.") + with pd.HDFStore(path, mode="r") as store: + raw = { + entity: read_frame_table(store, entity) + for entity in ("person", "household") + } + staged = _donor_extras(raw["person"], raw["household"]) + # Use the validated roster as the input mapping's householder flag, so + # stale engine flags cannot override the survey reference-person pointer. + raw["person"]["is_household_head"] = raw["person"][PERSON_ID_COLUMN].isin( + staged["reference_person_id"] + ) + tables = POLICYENGINE_US_CONCEPT_MAPPING.decode(raw) + del raw + for name, values in staged["pointers"].items(): + tables["person"][name] = values + tables["person"]["liquid_financial_assets"] = staged["liquid_financial_assets"] + tables["household"]["reference_person_id"] = staged["reference_person_id"] + violations = validate_concept_tables(tables) + if violations: + details = "; ".join(f"{v.entity}.{v.column}: {v.message}" for v in violations) + raise ValueError(f"Invalid donor concepts: {details}") + return DonorBank( + MappingProxyType(tables), + staged["weights"], + staged["support_strata"], + staged["source_person_ids"], + ) + + +def _source_id_payloads(source_ids) -> list[tuple[str, str]]: + values = np.asarray(source_ids, dtype=object) + if values.ndim != 1: + raise ValueError("Source person IDs must be one-dimensional.") + payloads = [] + for value in values: + if isinstance(value, str) and value: + payloads.append(("str", value)) + elif isinstance(value, (int, np.integer)) and not isinstance( + value, (bool, np.bool_) + ): + payloads.append(("int", str(value))) + else: + raise ValueError( + "Source person IDs must be nonempty strings or exact integers." + ) + return payloads + + +def derive_transport_seed(source_ids, stream: str) -> np.ndarray: + """Return order-independent float64 seeds in [0, 1) from exact source IDs. + + A domain-separated BLAKE2b hash binds the declared stream and typed source + identity. Repeated source IDs (support clones) share a seed. Keeping the + top 53 hash bits makes the result exactly representable, including the + largest possible draw, without a rounded value of 1. + """ + if not isinstance(stream, str) or not stream: + raise ValueError("A transport seed stream must be a nonempty string.") + payloads = _source_id_payloads(source_ids) + seeds = np.empty(len(payloads), dtype=np.float64) + for index, (kind, identity) in enumerate(payloads): + payload = json.dumps( + ["microcosm.transport.seed.v1", stream, kind, identity], + ensure_ascii=False, + separators=(",", ":"), + ).encode("utf-8") + bits = int.from_bytes(blake2b(payload, digest_size=8).digest(), "big") + seeds[index] = (bits >> 11) / 2**53 + return seeds + + +def currency_bridge( + tables: Mapping[str, pd.DataFrame], + columns: Mapping[str, Sequence[str]], + rate: float, +) -> dict[str, pd.DataFrame]: + """Copy tables and multiply each declared money column once, rounded to 2 dp. + + ``columns`` maps entities to their concept column names. Every selected + column must be a currency amount; duplicates and absent columns fail. + ``rate`` is supplied by the caller's specification, validated by the shared + ``apply_scale`` operator. Rounding uses NumPy's round-to-nearest-even rule. + Unselected columns and source tables are unchanged. + """ + if isinstance(rate, (bool, np.bool_)) or not np.isscalar(rate): + raise ValueError("A currency bridge requires one scalar rate.") + rate = float(rate) + apply_scale(np.empty(0), rate) + out = {entity: table.copy() for entity, table in tables.items()} + for entity, names in columns.items(): + if isinstance(names, str) or len(set(names)) != len(names): + raise ValueError("Currency bridge columns must be distinct column names.") + if entity not in out: + raise ValueError(f"Unknown currency bridge entity {entity!r}.") + for name in names: + item = concept_for_column(entity, name) + if item is None or item.unit is not Unit.BASE_CURRENCY: + raise ValueError( + f"Currency bridge requires an amount concept: {entity}.{name}." + ) + if name not in out[entity]: + raise ValueError(f"Currency bridge column is absent: {entity}.{name}.") + scaled = apply_scale(out[entity][name].to_numpy(dtype=np.float64), rate) + with np.errstate(over="ignore", invalid="ignore"): + rounded = np.round(scaled, 2) + if not np.isfinite(rounded).all(): + raise ValueError("Rounded currency amounts must be finite.") + out[entity][name] = rounded + return out diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py new file mode 100644 index 0000000000..7b68fab16c --- /dev/null +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py @@ -0,0 +1,558 @@ +"""Authenticated synthetic donor, exact roster, seed and currency contracts.""" + +import json +from hashlib import sha256 +from pathlib import Path +from types import SimpleNamespace +from unittest.mock import Mock + +import numpy as np +import pandas as pd +import pytest +from hypothesis import HealthCheck, given, settings +from hypothesis import strategies as st + +from microcosm.frame import SignedScale +from microcosm.frame.adapters.policyengine_us_concepts import ( + POLICYENGINE_US_CONCEPT_MAPPING, +) +from microcosm.frame.concept_mapping import ConceptMapping +from microcosm.frame.concepts import ContentBasis, validate_concept_tables +from microcosm.frame.transport import ( + currency_bridge, + derive_transport_seed, + quantile_map, + read_populace_us_donor, +) +from test_support.paths import paths_for + +FIXTURE = paths_for("microcosm-frame").tests / "fixtures/transport/synthetic_donor.json" +PROPERTY = settings(max_examples=75, deadline=None, database=None) +READER_PROPERTY = settings( + max_examples=40, + deadline=None, + database=None, + suppress_health_check=[HealthCheck.function_scoped_fixture], +) + + +def _fixture() -> dict: + return json.loads(FIXTURE.read_text()) + + +def _tables() -> dict[str, pd.DataFrame]: + return {entity: pd.DataFrame(rows) for entity, rows in _fixture()["tables"].items()} + + +def _pin(path: Path) -> dict: + return { + "size": path.stat().st_size, + "sha256": sha256(path.read_bytes()).hexdigest(), + } + + +def _write_h5(path: Path, tables: dict[str, pd.DataFrame]) -> dict: + with pd.HDFStore(path, mode="w") as store: + for entity, table in tables.items(): + store.put(entity, table, format="table") + return _pin(path) + + +def _fake_donor(monkeypatch, tmp_path, tables): + """Exercise public authentication/reader without repeated fixture H5 writes.""" + path = tmp_path / "synthetic-roster.h5" + path.write_bytes(b"synthetic donor bytes for mocked HDF decoding") + + class SyntheticStore: + def __init__(self, path, *, mode): + assert mode == "r" + + def __enter__(self): + return self + + def __exit__(self, *_): + return False + + def __getitem__(self, key): + return tables[key].copy(deep=True) + + def get_storer(self, _): + return SimpleNamespace(attrs=SimpleNamespace()) + + monkeypatch.setattr(pd, "HDFStore", SyntheticStore) + return read_populace_us_donor(path, **_pin(path)) + + +def test_synthetic_donor_matches_decode_and_exact_fixture(tmp_path, monkeypatch): + raw = _tables() + path = tmp_path / "synthetic.h5" + pin = _write_h5(path, raw) + original_decode = ConceptMapping.decode + seen = [] + + def spy_decode(mapping, tables): + seen.append({entity: table.copy(deep=True) for entity, table in tables.items()}) + return original_decode(mapping, tables) + + monkeypatch.setattr(ConceptMapping, "decode", spy_decode) + donor = read_populace_us_donor(path, **pin) + assert len(seen) == 1 + expected_decode = original_decode(POLICYENGINE_US_CONCEPT_MAPPING, seen[0]) + for entity, table in expected_decode.items(): + # Roster pointers use nullable Int64 even when the decoded reference + # pointer contains no nulls and pandas chooses a NumPy integer dtype. + pd.testing.assert_frame_equal( + donor.tables[entity][table.columns], table, check_dtype=False + ) + expected = _fixture()["expected"] + person = donor.tables["person"] + for column in ("partner_person_id", "parent_1_person_id", "parent_2_person_id"): + pd.testing.assert_series_equal( + person[column], + pd.Series(pd.array(expected[column], dtype="Int64"), name=column), + ) + np.testing.assert_array_equal( + donor.tables["household"]["reference_person_id"], + expected["reference_person_id"], + ) + np.testing.assert_array_equal( + person["liquid_financial_assets"], expected["liquid_financial_assets"] + ) + np.testing.assert_array_equal( + person["employment_income"], expected["employment_income"] + ) + np.testing.assert_array_equal(donor.weights, raw["household"]["household_weight"]) + np.testing.assert_array_equal( + donor.source_person_ids, raw["person"]["source_person_id"] + ) + assert donor.weights.dtype == np.float64 + assert donor.support_strata.tolist() == ["us:asec", "us:puf_tax_detail"] + assert donor.content_basis is ContentBasis.TRANSPORT + assert donor.donor_country == "us" + assert donor.currency == "USD" + assert validate_concept_tables(donor.tables) == () + assert "A_LINENO" not in person and "bank_account_assets" not in person + assert "household_weight" not in donor.tables["household"] + + +def test_synthetic_donor_transport_and_mapping_round_trip(tmp_path): + document = _fixture() + raw = _tables() + path = tmp_path / "roundtrip.h5" + donor = read_populace_us_donor(path, **_write_h5(path, raw)) + person = donor.tables["person"] + weights = ( + person["person_household_id"] + .map(pd.Series(donor.weights, index=donor.tables["household"]["household_id"])) + .to_numpy() + ) + before = {entity: table.copy(deep=True) for entity, table in donor.tables.items()} + original_weights = donor.weights.copy() + bands = document["target_bands"] + employment = quantile_map( + person["employment_income"], weights, bands, interpolation="uniform" + ) + assets = quantile_map( + person["liquid_financial_assets"], + weights, + bands, + interpolation="uniform", + group={"entity_ids": person["person_household_id"]}, + ) + np.testing.assert_allclose( + employment, document["expected"]["mapped_employment_income"] + ) + np.testing.assert_allclose( + assets, document["expected"]["mapped_liquid_financial_assets"] + ) + np.testing.assert_allclose( + quantile_map(employment, weights, bands, interpolation="uniform"), employment + ) + np.testing.assert_allclose( + quantile_map( + assets, + weights, + bands, + interpolation="uniform", + group={"entity_ids": person["person_household_id"]}, + ), + assets, + ) + transported = { + entity: table.copy(deep=True) for entity, table in donor.tables.items() + } + transported["person"]["employment_income"] = employment + encoded = POLICYENGINE_US_CONCEPT_MAPPING.encode( + transported, shares={"us.taxable_interest_share": 0.6} + ) + decoded = POLICYENGINE_US_CONCEPT_MAPPING.decode(encoded.tables) + for entity, table in decoded.items(): + pd.testing.assert_frame_equal( + table, transported[entity][table.columns], check_dtype=False + ) + for entity, table in before.items(): + pd.testing.assert_frame_equal(donor.tables[entity], table) + np.testing.assert_array_equal(donor.weights, original_weights) + + +@pytest.mark.parametrize("mismatch", ["size", "sha256"]) +def test_reader_refuses_pin_mismatch_before_hdf_parse(tmp_path, monkeypatch, mismatch): + path = tmp_path / "not-an-hdf.h5" + path.write_bytes(b"not an HDF file") + pin = _pin(path) + pin[mismatch] = pin[mismatch] + 1 if mismatch == "size" else "0" * 64 + parse = Mock(side_effect=AssertionError("HDF parsing preceded authentication")) + monkeypatch.setattr(pd, "HDFStore", parse) + with pytest.raises(ValueError, match="mismatch"): + read_populace_us_donor(path, **pin) + parse.assert_not_called() + + +@READER_PROPERTY +@given( + payload=st.binary(min_size=1, max_size=128), + excess_size=st.integers(min_value=1, max_value=1_000_000), +) +def test_reader_pin_mismatch_property_precedes_hdf_parse( + tmp_path, monkeypatch, payload, excess_size +): + path = tmp_path / "generated-donor.h5" + path.write_bytes(payload) + actual = _pin(path) + altered_digest = ("0" if actual["sha256"][0] != "0" else "1") + actual["sha256"][1:] + parse = Mock(side_effect=AssertionError("HDF parsing preceded authentication")) + monkeypatch.setattr(pd, "HDFStore", parse) + for pin in ( + {**actual, "size": actual["size"] + excess_size}, + {**actual, "sha256": altered_digest}, + ): + with pytest.raises(ValueError, match="mismatch"): + read_populace_us_donor(path, **pin) + parse.assert_not_called() + + +@pytest.mark.parametrize( + ("change", "message"), + [ + ("dangling_partner", "unknown household line"), + ("dangling_parent", "unknown household line"), + ("asymmetric_partner", "Partners must point at each other"), + ("duplicate_line", "unique within a household"), + ("multiple_reference", "multiple reference persons"), + ("missing_reference", "no reference person"), + ("channel_mismatch", "support channels disagree"), + ], +) +def test_reader_refuses_invalid_rosters(tmp_path, monkeypatch, change, message): + tables = _tables() + person = tables["person"] + mutations = { + "dangling_partner": (0, "A_SPOUSE", 99), + "dangling_parent": (2, "PEPAR1", 99), + "asymmetric_partner": (1, "A_SPOUSE", 0), + "duplicate_line": (1, "A_LINENO", 1), + "multiple_reference": (1, "A_EXPRRP", 1), + "missing_reference": (0, "A_EXPRRP", 5), + "channel_mismatch": (0, "person_support_channel", "another-channel"), + } + row, column, value = mutations[change] + person.loc[row, column] = value + with pytest.raises(ValueError, match=message): + _fake_donor(monkeypatch, tmp_path, tables) + + +@READER_PROPERTY +@given( + offset=st.integers(min_value=2**53 + 1, max_value=2**63 - 7), + order=st.permutations(tuple(range(6))), + reverse_households=st.booleans(), +) +def test_pointer_recode_is_exact_bijective_and_partner_symmetric( + tmp_path, monkeypatch, offset, order, reverse_households +): + tables = _tables() + tables["person"]["person_id"] = np.arange(offset, offset + 6, dtype=np.int64) + tables["person"] = tables["person"].iloc[list(order)].reset_index(drop=True) + if reverse_households: + tables["household"] = tables["household"].iloc[::-1].reset_index(drop=True) + donor = _fake_donor(monkeypatch, tmp_path, tables) + raw = tables["person"] + lookup = { + (row.person_household_id, row.A_LINENO): row.person_id + for row in raw.itertuples() + } + assert len(lookup) == len(raw) + concept = donor.tables["person"] + partners = concept.set_index("person_id")["partner_person_id"] + for row, recovered in zip(raw.itertuples(), concept.itertuples(), strict=True): + assert recovered.person_id == lookup[(row.person_household_id, row.A_LINENO)] + if row.A_SPOUSE: + partner = lookup[(row.person_household_id, row.A_SPOUSE)] + assert recovered.partner_person_id == partner + assert partners.loc[partner] == row.person_id + else: + assert pd.isna(recovered.partner_person_id) + parents = [ + lookup[(row.person_household_id, line)] + for line in (row.PEPAR1, row.PEPAR2) + if line + ] + for position, column in enumerate(("parent_1_person_id", "parent_2_person_id")): + observed = getattr(recovered, column) + if position < len(parents): + assert observed == parents[position] + else: + assert pd.isna(observed) + for row in donor.tables["household"].itertuples(): + assert row.reference_person_id == lookup[(row.household_id, 1)] + np.testing.assert_array_equal(donor.source_person_ids, raw["source_person_id"]) + np.testing.assert_array_equal( + donor.weights, tables["household"]["household_weight"] + ) + + +def test_reader_preserves_undifferentiated_us_support_and_exact_id_fallback( + tmp_path, monkeypatch +): + tables = _tables() + tables["person"] = tables["person"].drop( + columns=[ + "person_source_id", + "source_person_id", + "source_household_id", + "source_year", + "person_support_channel", + ] + ) + tables["household"] = tables["household"].drop(columns="household_support_channel") + donor = _fake_donor(monkeypatch, tmp_path, tables) + assert donor.support_strata.tolist() == ["us", "us"] + np.testing.assert_array_equal( + donor.source_person_ids, tables["person"]["person_id"] + ) + + +def test_reader_keeps_repeated_exact_source_ids_on_support_clones( + tmp_path, monkeypatch +): + tables = _tables() + tables["person"].loc[3, "person_source_id"] = tables["person"].loc[ + 0, "person_source_id" + ] + donor = _fake_donor(monkeypatch, tmp_path, tables) + assert donor.source_person_ids[0] == donor.source_person_ids[3] + seeds = derive_transport_seed(donor.source_person_ids, "synthetic-clone") + assert seeds[0] == seeds[3] + + +@pytest.mark.parametrize("scope", ["household", "year"]) +def test_reader_scopes_repeated_person_id_to_exact_source_lineage( + tmp_path, monkeypatch, scope +): + tables = _tables() + person = tables["person"].drop(columns="person_source_id") + tables["person"] = person + tables["household"]["household_support_channel"] = "asec" + person["person_support_channel"] = "asec" + exact_leaf = person.loc[0, "source_person_id"] + person.loc[3, "source_person_id"] = exact_leaf + if scope == "year": + person.loc[3, "source_household_id"] = person.loc[0, "source_household_id"] + person.loc[3, "source_year"] = person.loc[0, "source_year"] - 1 + donor = _fake_donor(monkeypatch, tmp_path, tables) + assert donor.source_person_ids[0] != donor.source_person_ids[3] + for index in (0, 3): + composite = json.loads(donor.source_person_ids[index]) + assert composite[0] == "us:asec" + assert composite[1] == ["int", str(person.loc[index, "source_year"])] + assert composite[2] == ["int", str(person.loc[index, "source_household_id"])] + assert composite[3] == ["str", exact_leaf] + seeds = derive_transport_seed(donor.source_person_ids, "synthetic-lineage") + assert seeds[0] != seeds[3] + + +@pytest.mark.parametrize( + "missing", ["source_year", "source_household_id", "source_person_id"] +) +def test_reader_refuses_partial_source_lineage(tmp_path, monkeypatch, missing): + tables = _tables() + tables["person"] = tables["person"].drop(columns=["person_source_id", missing]) + with pytest.raises(ValueError, match="source lineage requires"): + _fake_donor(monkeypatch, tmp_path, tables) + + +def test_reader_recodes_unmarried_reference_partner_symmetrically( + tmp_path, monkeypatch +): + tables = _tables() + person = tables["person"] + person.loc[[0, 1], "A_SPOUSE"] = 0 + person.loc[1, "A_EXPRRP"] = 13 + donor = _fake_donor(monkeypatch, tmp_path, tables) + pointers = donor.tables["person"]["partner_person_id"] + assert pointers.iloc[0] == person.loc[1, "person_id"] + assert pointers.iloc[1] == person.loc[0, "person_id"] + assert validate_concept_tables(donor.tables) == () + + +@pytest.mark.parametrize("conflict", ["multiple", "contradictory_spouse"]) +def test_reader_refuses_conflicting_unmarried_reference_partners( + tmp_path, monkeypatch, conflict +): + tables = _tables() + person = tables["person"] + person.loc[[0, 1], "A_SPOUSE"] = 0 + person.loc[1, "A_EXPRRP"] = 13 + if conflict == "multiple": + person.loc[2, "A_EXPRRP"] = 13 + else: + person.loc[0, "A_SPOUSE"] = 3 + person.loc[2, "A_SPOUSE"] = 1 + with pytest.raises(ValueError): + _fake_donor(monkeypatch, tmp_path, tables) + + +def test_cached_real_donor_size_only_when_present(): + pin = _fixture()["real_donor_pin"] + path = Path.home() / pin["cache_relative_path"] + if not path.is_file(): + pytest.skip("Pinned real donor is not locally cached") + assert path.stat().st_size == pin["size"] + + +@PROPERTY +@given( + ids=st.lists( + st.one_of(st.integers(min_value=0, max_value=2**100), st.text(min_size=1)), + max_size=30, + ), + stream=st.text(min_size=1), +) +def test_transport_seed_range_determinism_order_and_clone_invariance(ids, stream): + first = derive_transport_seed(ids, stream) + assert first.dtype == np.float64 + assert np.all((first >= 0) & (first < 1)) + np.testing.assert_array_equal(first, derive_transport_seed(ids, stream)) + np.testing.assert_array_equal(first[::-1], derive_transport_seed(ids[::-1], stream)) + np.testing.assert_array_equal( + np.repeat(first, 2), + derive_transport_seed(np.repeat(np.asarray(ids, dtype=object), 2), stream), + ) + + +def test_transport_seed_separates_streams_and_exact_typed_ids(): + ids = [2**53 + 1, 2**53 + 2, "0000000000000000000001", "1", 1] + first = derive_transport_seed(ids, "synthetic-a") + second = derive_transport_seed(ids, "synthetic-b") + assert len(np.unique(first)) == len(ids) + assert np.all(first != second) + + +@pytest.mark.parametrize("ids", [[None], [True], [1.0], [np.nan], [""], [[1]]]) +def test_transport_seed_refuses_inexact_or_missing_ids(ids): + with pytest.raises(ValueError, match="Source person IDs"): + derive_transport_seed(ids, "synthetic") + + +@pytest.mark.parametrize("stream", ["", None, 1]) +def test_transport_seed_requires_explicit_nonempty_stream(stream): + with pytest.raises(ValueError, match="stream"): + derive_transport_seed([1], stream) + + +@PROPERTY +@given( + values=st.lists( + st.floats(min_value=-1e8, max_value=1e8, allow_nan=False, allow_infinity=False), + min_size=1, + max_size=30, + ), + rate=st.floats( + min_value=0.01, max_value=100, allow_nan=False, allow_infinity=False + ), +) +def test_currency_bridge_matches_exactly_one_multiplication_and_preserves_inputs( + values, rate +): + tables = { + "person": pd.DataFrame( + { + "rental_income": values, + "age": np.arange(len(values)), + "employment_income": np.ones(len(values)), + } + ), + "household": pd.DataFrame( + {"household_id": [1], "rent": [10.0], "household_weight": [7.0]} + ), + } + before = {entity: table.copy(deep=True) for entity, table in tables.items()} + result = currency_bridge( + tables, {"person": ["rental_income"], "household": ["rent"]}, rate + ) + np.testing.assert_array_equal( + result["person"]["rental_income"], np.round(np.asarray(values) * rate, 2) + ) + np.testing.assert_array_equal( + result["household"]["rent"], np.round(np.array([10.0]) * rate, 2) + ) + pd.testing.assert_frame_equal( + result["person"].drop(columns="rental_income"), + before["person"].drop(columns="rental_income"), + ) + pd.testing.assert_frame_equal( + result["household"].drop(columns="rent"), + before["household"].drop(columns="rent"), + ) + for entity, table in before.items(): + pd.testing.assert_frame_equal(tables[entity], table) + + +@pytest.mark.parametrize( + "columns", + [ + {"person": ["age"]}, + {"person": ["unknown"]}, + {"person": ["rental_income", "rental_income"]}, + {"person": "rental_income"}, + {"missing": ["rent"]}, + ], +) +def test_currency_bridge_refuses_nonamount_missing_duplicate_or_unstructured_columns( + columns, +): + tables = { + "person": pd.DataFrame({"rental_income": [1.0], "age": [30]}), + "household": pd.DataFrame({"rent": [1.0]}), + } + with pytest.raises(ValueError): + currency_bridge(tables, columns, 1.25) + + +@pytest.mark.parametrize("rate", [0.0, -1.0, np.inf, np.nan]) +def test_currency_bridge_refuses_invalid_declared_rate(rate): + with pytest.raises(ValueError, match="finite and positive"): + currency_bridge( + {"person": pd.DataFrame({"rental_income": [1.0]})}, + {"person": ["rental_income"]}, + rate, + ) + + +@pytest.mark.parametrize("rate", [True, False, SignedScale(1.0, 2.0)]) +def test_currency_bridge_requires_one_numeric_scalar_rate(rate): + with pytest.raises(ValueError): + currency_bridge( + {"person": pd.DataFrame({"rental_income": [1.0]})}, + {"person": ["rental_income"]}, + rate, + ) + + +def test_currency_bridge_refuses_nonfinite_rounding_result(): + with pytest.raises(ValueError, match="finite"): + currency_bridge( + {"person": pd.DataFrame({"rental_income": [np.finfo(np.float64).max]})}, + {"person": ["rental_income"]}, + 1.0, + ) diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py new file mode 100644 index 0000000000..dccca8aa37 --- /dev/null +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py @@ -0,0 +1,484 @@ +"""Declared weighted transport distributions, without country engines.""" + +from math import fsum + +import numpy as np +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from microcosm.frame.transport import quantile_map + + +def _bands(): + # Synthetic test distribution; these are not country target values. + return [ + {"lower": 0.0, "upper": 10.0, "share": 0.25}, + {"lower": 10.0, "upper": 30.0, "share": 0.5}, + {"lower": 30.0, "upper": 90.0, "share": 0.25}, + ] + + +_samples = st.lists( + st.tuples(st.integers(-100, 100), st.integers(1, 40), st.integers(0, 2)), + min_size=1, + max_size=40, +) + + +def _reference_uniform(values, weights, bands, components): + """Deliberately quadratic independent definition of pooled weighted ranks.""" + result = [] + for value, component in zip(values, components, strict=True): + if value == 0: + result.append(0.0) + continue + peers = [ + (abs(peer), weight) + for peer, weight, peer_component in zip( + values, weights, components, strict=True + ) + if peer_component == component and peer * value > 0 + ] + rank = ( + fsum(weight for peer, weight in peers if peer < abs(value)) + + fsum(weight for peer, weight in peers if peer == abs(value)) / 2 + ) / fsum(weight for _, weight in peers) + start = 0.0 + for band in bands: + stop = start + band["share"] + if rank < stop: + within = (rank - start) / band["share"] + result.append( + np.sign(value) + * (band["lower"] + within * (band["upper"] - band["lower"])) + ) + break + start = stop + return np.array(result) + + +def test_synthetic_distribution_round_trip(): + values = np.array([1.0, 2.0, 3.0, 4.0]) + mapped = quantile_map(values, np.ones(4), _bands(), interpolation="uniform") + np.testing.assert_array_equal(mapped, [5.0, 15.0, 25.0, 60.0]) + np.testing.assert_array_equal( + quantile_map(mapped, np.ones(4), _bands(), interpolation="uniform"), mapped + ) + np.testing.assert_array_equal(values, [1.0, 2.0, 3.0, 4.0]) + + +def test_tied_values_share_a_pooled_midpoint_rank(): + values = [2.0, 1.0, 2.0] + weights = [1.0, 1.0, 2.0] + mapped = quantile_map(values, weights, _bands(), interpolation="uniform") + np.testing.assert_array_equal(mapped, [25.0, 5.0, 25.0]) + np.testing.assert_array_equal( + quantile_map(mapped, weights, _bands(), interpolation="uniform"), mapped + ) + + +def test_zero_and_sign_distributions_are_independent(): + bands = { + "positive": [{"lower": 10.0, "upper": 20.0, "share": 1.0}], + "negative": [{"lower": 30.0, "upper": 50.0, "share": 1.0}], + } + mapped = quantile_map( + [-1, -2, 0, 1, 2], [1, 1, 100, 1, 1], bands, interpolation="uniform" + ) + np.testing.assert_array_equal(mapped, [-35.0, -45.0, 0.0, 12.5, 17.5]) + + +def test_components_are_ranked_independently(): + mapped = quantile_map( + [1, 2, 100, 200], + [1, 1, 3, 3], + _bands(), + interpolation="uniform", + group=["a", "a", "b", "b"], + ) + np.testing.assert_array_equal(mapped, [10.0, 30.0, 10.0, 30.0]) + + +def test_entity_aggregation_maps_one_entity_weight_then_rescales_members(): + mapped = quantile_map( + [1, 3, 2, 0, 0], + [2, 2, 1, 1, 4], + [{"lower": 0, "upper": 90, "share": 1}], + interpolation="uniform", + group={"entity_ids": ["h1", "h1", "h2", "h2", "h3"]}, + ) + np.testing.assert_allclose(mapped, [15.0, 45.0, 15.0, 0.0, 0.0]) + assert mapped[1] / mapped[0] == 3 + + +def test_aggregation_accepts_components_and_preserves_cancelling_members(): + mapped = quantile_map( + [3, -1, 2, -2, 5, 0], + np.ones(6), + [{"lower": 10, "upper": 30, "share": 1}], + interpolation="uniform", + group={ + "entity_ids": ["a", "a", "b", "b", "c", "c"], + "components": [0, 0, 0, 0, 1, 1], + }, + ) + np.testing.assert_array_equal(mapped, [30.0, -10.0, 2.0, -2.0, 20.0, 0.0]) + + +@pytest.mark.parametrize("upper", [None, 100.0]) +def test_pareto_top_band_uses_explicit_alpha(upper): + bands = [{"lower": 10.0, "upper": upper, "share": 1.0}] + mapped = quantile_map( + [1.0, 2.0], [1, 1], bands, interpolation={"method": "pareto", "alpha": 2} + ) + tail = 0.0 if upper is None else (10 / upper) ** 2 + expected = 10 / np.sqrt(1 - np.array([0.25, 0.75]) * (1 - tail)) + np.testing.assert_allclose(mapped, expected) + np.testing.assert_array_equal( + quantile_map( + mapped, [1, 1], bands, interpolation={"method": "pareto", "alpha": 2} + ), + mapped, + ) + + +def test_pareto_only_changes_the_top_band(): + bands = [ + {"lower": 0, "upper": 10, "share": 0.5}, + {"lower": 10, "upper": None, "share": 0.5}, + ] + mapped = quantile_map( + [1, 2], [1, 1], bands, interpolation={"method": "pareto", "alpha": 1} + ) + np.testing.assert_array_equal(mapped, [5.0, 20.0]) + + +@given(_samples) +@settings(max_examples=100, deadline=None) +def test_rank_sign_zero_determinism_and_idempotence_properties(samples): + values, weights, components = map(np.array, zip(*samples, strict=True)) + mapped = quantile_map( + values, weights, _bands(), interpolation="uniform", group=components + ) + np.testing.assert_array_equal(np.sign(mapped), np.sign(values)) + np.testing.assert_array_equal(mapped[values == 0], values[values == 0]) + np.testing.assert_array_equal( + quantile_map( + values, weights, _bands(), interpolation="uniform", group=components + ), + mapped, + ) + np.testing.assert_array_equal( + quantile_map( + mapped, weights, _bands(), interpolation="uniform", group=components + ), + mapped, + ) + for component in set(components): + for sign in (-1, 1): + selected = (components == component) & (values * sign > 0) + order = np.argsort(values[selected], kind="stable") + assert np.all(np.diff(mapped[selected][order]) >= 0) + + +@given(_samples) +@settings(max_examples=100, deadline=None) +def test_weighted_band_shares_within_indivisible_tie_tolerance(samples): + values, weights, components = map(np.array, zip(*samples, strict=True)) + mapped = quantile_map( + values, weights, _bands(), interpolation="uniform", group=components + ) + for component in set(components): + for sign in (-1, 1): + selected = (components == component) & (values * sign > 0) + if not selected.any(): + continue + originals, magnitudes, mass = ( + abs(values[selected]), + abs(mapped[selected]), + weights[selected], + ) + total = mass.sum() + tolerance = ( + max(mass[originals == value].sum() for value in set(originals)) / total + ) + for position, band in enumerate(_bands()): + inside = (magnitudes >= band["lower"]) & ( + magnitudes <= band["upper"] + if position == 2 + else magnitudes < band["upper"] + ) + share = mass[inside].sum() / total + assert abs(share - band["share"]) <= tolerance + 1e-12 + + +@given(_samples) +@settings(max_examples=100, deadline=None) +def test_uniform_map_matches_independent_reference(samples): + values, weights, components = map(np.array, zip(*samples, strict=True)) + mapped = quantile_map( + values, weights, _bands(), interpolation="uniform", group=components + ) + np.testing.assert_allclose( + mapped, + _reference_uniform(values, weights, _bands(), components), + rtol=1e-12, + atol=1e-12, + ) + + +@given(_samples, st.integers(1, 100)) +@settings(max_examples=80, deadline=None) +def test_mapping_preserves_weights_and_is_invariant_to_their_common_scale( + samples, factor +): + values, weights, components = map(np.array, zip(*samples, strict=True)) + original_weights = weights.copy() + first = quantile_map( + values, weights, _bands(), interpolation="uniform", group=components + ) + second = quantile_map( + values, weights * factor, _bands(), interpolation="uniform", group=components + ) + np.testing.assert_allclose(second, first, rtol=1e-12, atol=1e-12) + np.testing.assert_array_equal(weights, original_weights) + + +@given(_samples) +@settings(max_examples=80, deadline=None) +def test_permuting_members_preserves_mapped_values(samples): + values, weights, components = map(np.array, zip(*samples, strict=True)) + first = quantile_map( + values, weights, _bands(), interpolation="uniform", group=components + ) + second = quantile_map( + values[::-1], + weights[::-1], + _bands(), + interpolation="uniform", + group=components[::-1], + ) + np.testing.assert_allclose(second[::-1], first, rtol=1e-12, atol=1e-12) + + +@given( + st.lists( + st.lists(st.integers(0, 50), min_size=1, max_size=5), min_size=1, max_size=12 + ) +) +@settings(max_examples=100, deadline=None) +def test_aggregate_pro_rata_totals_and_idempotence(entities): + values = np.array([value for members in entities for value in members], dtype=float) + entity_ids = np.array( + [entity for entity, members in enumerate(entities) for _ in members] + ) + group = {"entity_ids": entity_ids} + mapped = quantile_map( + values, np.ones(len(values)), _bands(), interpolation="uniform", group=group + ) + totals = np.array([fsum(members) for members in entities]) + expected = quantile_map( + totals, np.ones(len(entities)), _bands(), interpolation="uniform" + ) + actual = np.array( + [fsum(mapped[entity_ids == entity]) for entity in range(len(entities))] + ) + np.testing.assert_allclose(actual, expected, rtol=1e-12, atol=1e-12) + np.testing.assert_allclose( + quantile_map( + mapped, np.ones(len(values)), _bands(), interpolation="uniform", group=group + ), + mapped, + rtol=1e-12, + atol=1e-12, + ) + for entity, total in enumerate(totals): + if total: + selected = (entity_ids == entity) & (values != 0) + np.testing.assert_allclose( + mapped[selected] / values[selected], expected[entity] / total + ) + + +@given(st.integers(2, 10000), st.integers(1, 9)) +@settings(max_examples=100, deadline=None) +def test_tied_entity_totals_remain_tied_on_repeated_mapping(total, divider): + split = max(1, total // divider) + values = [float(split), float(total - split), float(total)] + group = {"entity_ids": [0, 0, 1]} + mapped = quantile_map( + values, [1, 1, 1], _bands(), interpolation="uniform", group=group + ) + np.testing.assert_allclose( + quantile_map(mapped, [1, 1, 1], _bands(), interpolation="uniform", group=group), + mapped, + rtol=1e-12, + atol=1e-12, + ) + + +@given( + st.lists(st.integers(1, 100), min_size=1, max_size=30), + st.floats(0.5, 5, allow_nan=False, allow_infinity=False), +) +@settings(max_examples=80, deadline=None) +def test_pareto_rank_and_idempotence_properties(values, alpha): + bands = [{"lower": 10, "upper": None, "share": 1}] + interpolation = {"method": "pareto", "alpha": alpha} + mapped = quantile_map( + values, np.ones(len(values)), bands, interpolation=interpolation + ) + order = np.argsort(values) + assert np.all(np.diff(mapped[order]) >= 0) + np.testing.assert_array_equal( + quantile_map(mapped, np.ones(len(values)), bands, interpolation=interpolation), + mapped, + ) + + +@pytest.mark.parametrize( + "values, weights", + [ + ([1, np.nan], [1, 1]), + ([np.inf], [1]), + ([1], [0]), + ([1], [-1]), + ([1], [np.inf]), + ([[1]], [[1]]), + ([1, 2], [1]), + ], +) +def test_invalid_values_and_weights_fail_closed(values, weights): + with pytest.raises(ValueError): + quantile_map(values, weights, _bands(), interpolation="uniform") + + +@pytest.mark.parametrize( + "bands", + [ + [], + [{"lower": -1, "upper": 1, "share": 1}], + [{"lower": 1, "upper": 1, "share": 1}], + [{"lower": 0, "upper": 1, "share": 0}], + [{"lower": 0, "upper": 1, "share": 0.5}], + [ + {"lower": 0, "upper": 1, "share": 0.5}, + {"lower": 2, "upper": 3, "share": 0.5}, + ], + [{"lower": 0, "upper": None, "share": 1}], + {"positive": _bands()}, + [{"lower": 0, "upper": 1, "share": 1, "other": 3}], + ], +) +def test_invalid_target_band_contract_fails_closed(bands): + with pytest.raises(ValueError): + quantile_map([1], [1], bands, interpolation="uniform") + + +@pytest.mark.parametrize( + "interpolation", + [ + "pareto", + "normal", + {"method": "pareto", "alpha": 0}, + {"method": "pareto", "alpha": np.inf}, + {"method": "uniform", "alpha": 2}, + {"method": "pareto"}, + ], +) +def test_invalid_interpolation_fails_closed(interpolation): + with pytest.raises(ValueError): + quantile_map([1], [1], _bands(), interpolation=interpolation) + + +def test_pareto_requires_a_positive_top_lower_bound(): + with pytest.raises(ValueError, match="positive lower"): + quantile_map( + [1], + [1], + [{"lower": 0, "upper": None, "share": 1}], + interpolation={"method": "pareto", "alpha": 2}, + ) + + +@pytest.mark.parametrize( + "group", + [ + [1], + [None, 1], + [[1], [2]], + {"components": [1, 2]}, + {"entity_ids": [1, 1], "components": [1, 2]}, + {"entity_ids": [1, 2], "other": [1, 2]}, + ], +) +def test_invalid_group_contract_fails_closed(group): + with pytest.raises(ValueError): + quantile_map([1, 2], [1, 1], _bands(), interpolation="uniform", group=group) + + +def test_aggregate_rejects_multiple_weights_for_one_entity(): + with pytest.raises(ValueError, match="one entity weight"): + quantile_map( + [1, 2], + [1, 2], + _bands(), + interpolation="uniform", + group={"entity_ids": [1, 1]}, + ) + + +def test_empty_population_is_supported(): + assert quantile_map([], [], _bands(), interpolation="uniform").shape == (0,) + assert quantile_map( + [], [], _bands(), interpolation="uniform", group={"entity_ids": []} + ).shape == (0,) + + +def test_aggregate_refuses_severe_cancellation_instead_of_distorting_small_member(): + values = [1e16, -(1e16 - 2), 1] + bands = [{"lower": 0, "upper": 36, "share": 1}] + with pytest.raises(ValueError, match="float64"): + quantile_map( + values, + [1, 1, 1], + bands, + interpolation="uniform", + group={"entity_ids": [0, 0, 0]}, + ) + + +def test_float64_rank_collapse_is_refused(): + bands = [{"lower": 1e16, "upper": np.nextafter(1e16, np.inf), "share": 1}] + with pytest.raises(ValueError, match="distinct ranks"): + quantile_map([1, 2, 3], [1, 1, 1], bands, interpolation="uniform") + + +def test_unrepresentable_weight_ratio_is_refused(): + with pytest.raises(ValueError, match="Weight ratios"): + quantile_map([1, 2], [1e308, 1e-308], _bands(), interpolation="uniform") + + +def test_pro_rata_mapping_avoids_overflow_in_the_common_scale_factor(): + bands = [{"lower": 0, "upper": 1e308, "share": 1}] + mapped = quantile_map( + [1e-308, 1e-308], + [1, 1], + bands, + interpolation="uniform", + group={"entity_ids": [0, 0]}, + ) + np.testing.assert_array_equal(mapped, [2.5e307, 2.5e307]) + + +def test_cancelling_members_refuse_an_unrepresentable_aggregate_total(): + bands = [{"lower": 10, "upper": 30, "share": 1}] + with pytest.raises(ValueError, match="exact totals"): + quantile_map( + [1e16, -9999999999999998.0], + [1, 1], + bands, + interpolation="uniform", + group={"entity_ids": [0, 0]}, + ) diff --git a/packages/microcosm-frame/tests/fixtures/transport/synthetic_donor.json b/packages/microcosm-frame/tests/fixtures/transport/synthetic_donor.json new file mode 100644 index 0000000000..c19959aeab --- /dev/null +++ b/packages/microcosm-frame/tests/fixtures/transport/synthetic_donor.json @@ -0,0 +1,236 @@ +{ + "description": "Synthetic US donor support records; these are never New Zealand microdata.", + "real_donor_pin": { + "sha256": "48b9d479fb4fd1c3537f9383ce4697d130b6f618658409d74f6233c43b994c7e", + "size": 462915783, + "cache_relative_path": ".cache/huggingface/hub/datasets--policyengine--populace-us/blobs/48b9d479fb4fd1c3537f9383ce4697d130b6f618658409d74f6233c43b994c7e" + }, + "tables": { + "person": [ + { + "person_id": 9007199254741001, + "person_household_id": 101, + "source_person_id": "0000000000000000000001", + "person_support_channel": "asec", + "A_LINENO": 1, + "A_SPOUSE": 2, + "PEPAR1": 0, + "PEPAR2": 0, + "A_EXPRRP": 1, + "age": 40, + "is_female": false, + "employment_income_before_lsr": 0, + "taxable_interest_income": 1, + "tax_exempt_interest_income": 2, + "pre_subsidy_rent": 25, + "bank_account_assets": 10, + "stock_assets": 20, + "bond_assets": 30, + "source_year": 2024, + "source_household_id": 101, + "person_source_id": "0000000000000000000001" + }, + { + "person_id": 9007199254741002, + "person_household_id": 101, + "source_person_id": "0000000000000000000002", + "person_support_channel": "asec", + "A_LINENO": 2, + "A_SPOUSE": 1, + "PEPAR1": 0, + "PEPAR2": 0, + "A_EXPRRP": 4, + "age": 38, + "is_female": true, + "employment_income_before_lsr": 10, + "taxable_interest_income": 2, + "tax_exempt_interest_income": 3, + "pre_subsidy_rent": 0, + "bank_account_assets": 0, + "stock_assets": 0, + "bond_assets": 0, + "source_year": 2024, + "source_household_id": 101, + "person_source_id": "0000000000000000000002" + }, + { + "person_id": 9007199254741003, + "person_household_id": 101, + "source_person_id": "0000000000000000000003", + "person_support_channel": "asec", + "A_LINENO": 3, + "A_SPOUSE": 0, + "PEPAR1": 1, + "PEPAR2": 2, + "A_EXPRRP": 5, + "age": 10, + "is_female": false, + "employment_income_before_lsr": 20, + "taxable_interest_income": 0, + "tax_exempt_interest_income": 0, + "pre_subsidy_rent": 0, + "bank_account_assets": 0, + "stock_assets": 0, + "bond_assets": 0, + "source_year": 2024, + "source_household_id": 101, + "person_source_id": "0000000000000000000003" + }, + { + "person_id": 9007199254741004, + "person_household_id": 202, + "source_person_id": "0000000000000000000004", + "person_support_channel": "puf_tax_detail", + "A_LINENO": 1, + "A_SPOUSE": 2, + "PEPAR1": 0, + "PEPAR2": 0, + "A_EXPRRP": 1, + "age": 45, + "is_female": false, + "employment_income_before_lsr": 30, + "taxable_interest_income": 3, + "tax_exempt_interest_income": 4, + "pre_subsidy_rent": 0, + "bank_account_assets": 30, + "stock_assets": 60, + "bond_assets": 90, + "source_year": 2024, + "source_household_id": 202, + "person_source_id": "0000000000000000000004" + }, + { + "person_id": 9007199254741005, + "person_household_id": 202, + "source_person_id": "0000000000000000000005", + "person_support_channel": "puf_tax_detail", + "A_LINENO": 2, + "A_SPOUSE": 1, + "PEPAR1": 0, + "PEPAR2": 0, + "A_EXPRRP": 4, + "age": 42, + "is_female": true, + "employment_income_before_lsr": 40, + "taxable_interest_income": 0, + "tax_exempt_interest_income": 0, + "pre_subsidy_rent": 0, + "bank_account_assets": 0, + "stock_assets": 0, + "bond_assets": 0, + "source_year": 2024, + "source_household_id": 202, + "person_source_id": "0000000000000000000005" + }, + { + "person_id": 9007199254741006, + "person_household_id": 202, + "source_person_id": "0000000000000000000006", + "person_support_channel": "puf_tax_detail", + "A_LINENO": 3, + "A_SPOUSE": 0, + "PEPAR1": 0, + "PEPAR2": 2, + "A_EXPRRP": 5, + "age": 12, + "is_female": true, + "employment_income_before_lsr": 50, + "taxable_interest_income": 0, + "tax_exempt_interest_income": 0, + "pre_subsidy_rent": 0, + "bank_account_assets": 0, + "stock_assets": 0, + "bond_assets": 0, + "source_year": 2024, + "source_household_id": 202, + "person_source_id": "0000000000000000000006" + } + ], + "household": [ + { + "household_id": 101, + "household_weight": 1, + "household_support_channel": "asec" + }, + { + "household_id": 202, + "household_weight": 3, + "household_support_channel": "puf_tax_detail" + } + ] + }, + "expected": { + "partner_person_id": [ + 9007199254741002, + 9007199254741001, + null, + 9007199254741005, + 9007199254741004, + null + ], + "parent_1_person_id": [ + null, + null, + 9007199254741001, + null, + null, + 9007199254741005 + ], + "parent_2_person_id": [ + null, + null, + 9007199254741002, + null, + null, + null + ], + "reference_person_id": [ + 9007199254741001, + 9007199254741004 + ], + "liquid_financial_assets": [ + 60, + 0, + 0, + 180, + 0, + 0 + ], + "employment_income": [ + 0, + 10, + 20, + 30, + 40, + 50 + ], + "mapped_employment_income": [ + 0, + 109.0909090909091, + 127.27272727272727, + 163.63636363636363, + 236.36363636363637, + 345.45454545454544 + ], + "mapped_liquid_financial_assets": [ + 125, + 0, + 0, + 250, + 0, + 0 + ] + }, + "target_bands": [ + { + "lower": 100, + "upper": 200, + "share": 0.5 + }, + { + "lower": 200, + "upper": 400, + "share": 0.5 + } + ] +} From 3522cad825b9b132839f0655a2deb7fa41102f76 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Wed, 7 Oct 2026 14:10:25 -0400 Subject: [PATCH 2/6] Read CPS pointers as the pinned donor writes them; fix pro rata ties (#1130) The donor reader refused the real file it is pinned to. Treat a line pointer that is missing or <= 0 as absent (the pinned donor writes -1 for an absent parent or cohabiting partner and 0 for an absent spouse), as the US build's _parent_person_ids does for PEPAR1/2, and use the -1 sentinel in the synthetic fixture. Take partners only from A_SPOUSE, or from PECOHAB on rows whose vintage carries it, and refuse the two naming different people. Drop the relationship-code-13 recovery: code 13 also covers roommates and housemates (docs/us-asec-census-person-columns.md, COVERAGE_PROGRESS.md). On the pinned donor only the 2024 rows carry PECOHAB. Take the reference person from A_EXPRRP 1 or 2, fall back to is_household_head where a household has none, and refuse the whole donor if any household's two disagree. On the pinned donor 56,943 households have an A_EXPRRP reference person and 297 take the fallback; 113 of the 56,943 disagree, so the reader refuses that file until the hub rules. All 113 are 2022/2023 rows, whose A_EXPRRP equals the build's line-1 derivation (asec_pool._with_relationship_recode); is_household_head is P_SEQ == 1 on every row. A new smoke test runs the full reader on the pinned blob, skips when it is absent, and asserts that exact refusal, so any earlier failure (the old -1 refusal, a pin mismatch) turns it red. The is_household_head override before decode is dropped: with the disagreement refused, it rewrote every flag to its own value. Fix aggregate quantile_map refusing ordinary members, which made the PR's own Hypothesis property fail (8 of 10 seeds; the engine-free CI job's one failure). Correcting only the largest member cannot reach a total that the other members leave exactly halfway between two of its sums (members [1, 2] mapped to 0.21), nor one the residual step overshoots by several units (13 members mapped to 8192, where the spacing halves). After the old correction fails, restart the largest member from the correctly rounded difference (fsum), and if needed move one other member by one unit in the last place: up to four members finer than the largest, largest first, then the smallest, each away from and then toward zero. Inputs the old correction accepted map bit-identically, same-sign members always reach the exact total, and the refusal path stays O(n). Pin the counterexamples and add a same-sign exact-total property. Add a fault-injected test for the pro rata correction bound, which the exact-total check otherwise always caught first. Document pro rata amplification of nearly cancelling mixed-sign members, the remaining mixed-sign refusals, and the hash-then-reopen gap. Update the ADR status and changelog to say what the reader does on the pinned file. Co-Authored-By: Claude Opus 5.5 --- changelog.d/nz-transport-operators.added.md | 11 +- docs/concept-schema-transport-adr.md | 30 +- .../src/microcosm/frame/transport.py | 265 +++++++++++++----- .../shared/test_transport_donor.py | 182 ++++++++++-- .../shared/test_transport_quantile.py | 145 +++++++++- .../fixtures/transport/synthetic_donor.json | 30 +- 6 files changed, 538 insertions(+), 125 deletions(-) diff --git a/changelog.d/nz-transport-operators.added.md b/changelog.d/nz-transport-operators.added.md index 75ad7b707f..3bb786ed7f 100644 --- a/changelog.d/nz-transport-operators.added.md +++ b/changelog.d/nz-transport-operators.added.md @@ -1,3 +1,8 @@ -Add engine-free transport operators for authenticated local US donor concepts, -stable source-ID seeds, one-step currency conversion, and weighted signed -quantile maps with component and entity-aggregate support. +Add engine-free transport operators: a reader for local populace-US donor +files pinned by size and SHA-256, stable source-ID seeds, one-step currency +conversion, and weighted signed quantile maps with component and +entity-aggregate support. Run on the donor its tests pin, the reader resolves +every line pointer but refuses the file: in 113 of its 57,240 households the +`A_EXPRRP` reference person is not the `is_household_head` person. Once that +is resolved, concept validation would also refuse two households whose +parent pointers form a cycle. diff --git a/docs/concept-schema-transport-adr.md b/docs/concept-schema-transport-adr.md index 13bf329550..fd22a2e667 100644 --- a/docs/concept-schema-transport-adr.md +++ b/docs/concept-schema-transport-adr.md @@ -1,10 +1,18 @@ # Engine-neutral concepts and transport (ADR) **Status:** accepted for the concept schema, the per-engine mappings and their -coverage reports. The transport operators now include a pinned local US donor -reader, source-ID seeds, a currency bridge and weighted quantile mapping; -the composed transport pipeline below remains a design, and no build has -migrated. **Date:** 2026-09-28. **Source:** Max's 27 September rulings on +coverage reports. `microcosm.frame.transport` adds a reader for local +populace-US donor files pinned by size and SHA-256, source-ID seeds, a +currency bridge and weighted quantile mapping. The reader does not yet load +the donor its tests pin. Run on that file (166,321 persons, 57,240 +households), it resolves every line pointer, then refuses the whole file: in +113 households the `A_EXPRRP` reference person is not the `is_household_head` +person, and the reader will not choose between them. All 113 are 2022 or 2023 +rows, whose `A_EXPRRP` matches the build's line-1 derivation rather than the +Census recode. Once that is resolved, concept validation would also refuse +two other households, copies of one source household in which two people name +each other as parents. The composed transport pipeline below remains a +design, and no build has migrated. **Date:** 2026-09-28. **Source:** Max's 27 September rulings on international populations and law-anchored concepts. **Amended 2026-10-07:** schema version 2 adds `liquid_financial_assets`, the first stock concept, for New Zealand's Accommodation Supplement cash-asset test. Only the Axiom New @@ -285,11 +293,15 @@ These hold for every input and are tested (Hypothesis properties unless noted): encoding it through the adapter's mapping, one build at a time. - The spec engine's `imputation.yaml` `concepts:` block is unrelated. It groups predictor columns for imputation, and its names are not concept ids. -- `microcosm.frame.transport` supplies the donor bank and the income - transformation operators (step 3). Its bank declares US donor provenance - and source-currency amounts; it does not declare destination residents. - The consuming builder supplies the destination frame declaration after - transport. Transport still needs a unit-construction step that builds engine +- `microcosm.frame.transport` supplies a reader for pinned local donor files + and the income transformation operators (step 3). The reader's bank + declares US donor provenance and source-currency amounts; it does not + declare destination residents. It produces no bank from the donor its + tests pin yet: that file stops at the reference-person disagreement in the + status line above, and two households with a parent-pointer cycle would + then fail concept validation. The consuming builder supplies the + destination frame declaration after transport. Transport still needs a + unit-construction step that builds engine units from the concept pointers (until it exists, 13 of New Zealand's 80 and 46 of Belgium's 105 module bindings, and the take-up and housing bindings on US and UK group entities, stay deferred), and a country pack with geography, diff --git a/packages/microcosm-frame/src/microcosm/frame/transport.py b/packages/microcosm-frame/src/microcosm/frame/transport.py index e37888a7cc..7e4dea2cd9 100644 --- a/packages/microcosm-frame/src/microcosm/frame/transport.py +++ b/packages/microcosm-frame/src/microcosm/frame/transport.py @@ -149,6 +149,38 @@ def _quantile_component( return mapped[inverse] +def _retain_exact_total( + result: np.ndarray, selected: np.ndarray, anchor: int, total: float +) -> bool: + """Try to correct ``result[anchor]`` in place so the members sum to total. + + The residual step uses the rounded current sum, which can miss by several + units of a small anchor's spacing (a total at a power of two halves its + spacing below it). If it and two one-unit steps fail, the anchor restarts + from the correctly rounded difference between the total and the other + members, then takes up to two more one-unit steps. + """ + for attempt in range(3): + current_total = fsum(result[selected]) + if current_total == total: + return True + if attempt == 0: + result[anchor] += total - current_total + else: + direction = np.inf if current_total < total else -np.inf + result[anchor] = np.nextafter(result[anchor], direction) + if fsum(result[selected]) == total: + return True + result[anchor] = fsum([total, *(-result[selected[selected != anchor]])]) + for _ in range(2): + current_total = fsum(result[selected]) + if current_total == total: + return True + direction = np.inf if current_total < total else -np.inf + result[anchor] = np.nextafter(result[anchor], direction) + return fsum(result[selected]) == total + + def _quantile_map_components( values: np.ndarray, weights: np.ndarray, @@ -211,11 +243,23 @@ def quantile_map( once, and members are rescaled pro rata by a positive factor. A floating rounding residual is assigned to the largest nonzero member so that tied entity totals remain exactly tied on repeated mapping. This also - preserves member signs for mixed-sign assets. A zero-total entity is - unchanged, including cancelling nonzero members. Component labels and - entity ids must be nonmissing hashable scalars. Mapping fails if float64 - cannot retain distinct ranks, signs, pro rata precision, or the mapped - aggregate totals (including severe mixed-sign cancellation). + preserves member signs for mixed-sign assets. Where the largest member + alone cannot reach the exact total (members ``[1, 2]`` whose total maps + to 0.21), a member with finer float64 spacing first moves by one unit in + the last place. A zero-total entity is unchanged, including cancelling + nonzero members. + The pro rata factor is the mapped total over the entity total and has no + upper bound: members that nearly cancel are amplified with it, so + members ``[1000, -999]`` whose total of 1 maps to 15 become ``[15000, + -14985]``. Mixed-sign members should not be aggregated when the entity + total is small relative to their magnitudes; ``liquid_financial_assets`` + (``lower=0``) cannot mix signs. Component labels and entity ids must be + nonmissing hashable scalars. Mapping fails if float64 cannot retain + distinct ranks, signs, pro rata precision, or the mapped aggregate totals + (including severe mixed-sign cancellation). Same-sign members always + reach the exact total. Mixed-sign members can still be refused when no + tried one-unit move reaches it, as for members ``[-5, 1, 1]`` whose total + of -3 maps to -20. """ values = np.asarray(values, dtype=np.float64) weights = np.asarray(weights, dtype=np.float64) @@ -298,18 +342,40 @@ def quantile_map( # cancellation may make the exact total unrepresentable. anchor = nonzero[np.argmax(np.abs(result[nonzero]))] original_anchor = result[anchor] - for attempt in range(3): - current_total = fsum(result[selected]) - if current_total == mapped[entity]: + retained = _retain_exact_total(result, selected, anchor, mapped[entity]) + # Smaller members can leave the exact total halfway between + # two sums the largest member reaches (members [1, 2] whose + # total maps to 0.21). A member with finer float64 spacing + # than that tie moving by one unit in the last place, a + # relative change near 1e-16, breaks it. Up to four members + # finer than the largest are tried, largest first, then the + # smallest member: the tie can sit a binade above the largest + # member, and the smallest has the finest spacing. Each moves + # away from zero and then toward it. + candidates = [] + if not retained: + others = nonzero[nonzero != anchor] + others = others[np.argsort(-np.abs(result[others]), kind="stable")] + finer = np.spacing(np.abs(result[others])) < np.spacing( + abs(original_anchor) + ) + candidates = list(dict.fromkeys([*others[finer][:4], others[-1]])) + for other, toward in [ + (other, side * np.sign(result[other])) + for other in candidates + for side in (np.inf, 0.0) + ]: + if retained: break - if attempt == 0: - result[anchor] += mapped[entity] - current_total - else: - direction = ( - np.inf if current_total < mapped[entity] else -np.inf - ) - result[anchor] = np.nextafter(result[anchor], direction) - if fsum(result[selected]) != mapped[entity]: + kept = result[other] + result[anchor] = original_anchor + result[other] = np.nextafter(kept, toward) + retained = result[other] != 0 and _retain_exact_total( + result, selected, anchor, mapped[entity] + ) + if not retained: + result[other] = kept + if not retained: raise ValueError( "Pro rata mapping cannot retain exact totals at float64 precision." ) @@ -350,19 +416,39 @@ class DonorBank: currency: str = "USD" -def _line_numbers(values: pd.Series, *, nullable: bool = False) -> np.ndarray: +def _line_numbers(values: pd.Series, *, pointer: bool = False) -> np.ndarray: numbers = pd.to_numeric(values, errors="raise").to_numpy( dtype=np.float64, na_value=np.nan ) - if nullable: + if pointer: + # A missing pointer is absent (rows from a vintage without PECOHAB + # leave it NaN). numbers = np.where(np.isnan(numbers), 0, numbers) if not np.isfinite(numbers).all() or np.any(numbers != np.floor(numbers)): raise ValueError(f"Donor {values.name} must contain integer line numbers.") - if np.any(numbers < (0 if nullable else 1)) or np.any(numbers >= 2**63): + if pointer: + # The pinned donor writes -1 for an absent parent or cohabiting + # partner and 0 for an absent spouse. Every pointer <= 0 is absent, + # as the US build's eligibility_inputs._parent_person_ids does for + # PEPAR1/2. + numbers = np.maximum(numbers, 0) + if np.any(numbers < (0 if pointer else 1)) or np.any(numbers >= 2**63): raise ValueError(f"Donor {values.name} contains invalid line numbers.") return numbers.astype(np.int64) +def _pointed_rows( + values: pd.Series, roster: pd.MultiIndex, household_ids: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """Each line pointer's presence and its target row in the household roster.""" + pointed_line = _line_numbers(values, pointer=True) + present = pointed_line != 0 + rows = roster.get_indexer(pd.MultiIndex.from_arrays([household_ids, pointed_line])) + if np.any(present & (rows < 0)): + raise ValueError(f"Donor {values.name} points to an unknown household line.") + return present, rows + + def _donor_extras(person: pd.DataFrame, household: pd.DataFrame) -> dict: """Translate donor-only columns before crossing the concept decode seam.""" required = ( @@ -371,6 +457,7 @@ def _donor_extras(person: pd.DataFrame, household: pd.DataFrame) -> dict: "PEPAR1", "PEPAR2", "A_EXPRRP", + "is_household_head", "bank_account_assets", "stock_assets", "bond_assets", @@ -395,62 +482,70 @@ def _donor_extras(person: pd.DataFrame, household: pd.DataFrame) -> dict: if (person[PERSON_ID_COLUMN] > np.iinfo(np.int64).max).any(): raise ValueError("Donor person IDs must fit the concept pointer's int64 dtype.") person_ids = person[PERSON_ID_COLUMN].to_numpy(dtype=np.int64) - pointers = {} - for source, destination in ( - ("A_SPOUSE", "partner_person_id"), - ("PEPAR1", "parent_1_person_id"), - ("PEPAR2", "parent_2_person_id"), - ): - pointed_line = _line_numbers(person[source], nullable=True) - present = pointed_line != 0 - rows = roster.get_indexer( - pd.MultiIndex.from_arrays([household_ids, pointed_line]) + + def pointer(present: np.ndarray, rows: np.ndarray) -> pd.arrays.IntegerArray: + values = pd.array([pd.NA] * len(person), dtype="Int64") + values[present] = person_ids[rows[present]] + return values + + # Partners come only from the spouse and cohabiting-partner line pointers. + # Relationship code 13 is never read: the Census recode uses it for + # roommates and housemates as well as unmarried partners, and the + # ACS-to-CPS crosswalk maps all three to it, so it does not identify a + # partner. Rows from a vintage without PECOHAB leave it NaN, so their + # cohabiting couples get no partner pointer. + partner_present, partner_rows = _pointed_rows( + person["A_SPOUSE"], roster, household_ids + ) + if "PECOHAB" in person: + cohabiting, cohabitant_rows = _pointed_rows( + person["PECOHAB"], roster, household_ids ) - if np.any(present & (rows < 0)): - raise ValueError(f"Donor {source} points to an unknown household line.") - pointer = pd.array([pd.NA] * len(person), dtype="Int64") - pointer[present] = person_ids[rows[present]] - pointers[destination] = pointer + if np.any(partner_present & cohabiting & (partner_rows != cohabitant_rows)): + raise ValueError("Donor A_SPOUSE and PECOHAB name different partners.") + partner_rows = np.where(partner_present, partner_rows, cohabitant_rows) + partner_present = partner_present | cohabiting + pointers = { + "partner_person_id": pointer(partner_present, partner_rows), + "parent_1_person_id": pointer( + *_pointed_rows(person["PEPAR1"], roster, household_ids) + ), + "parent_2_person_id": pointer( + *_pointed_rows(person["PEPAR2"], roster, household_ids) + ), + } # A sole co-resident parent is always parent 1 in the concept contract. move = ( pointers["parent_1_person_id"].isna() & ~pointers["parent_2_person_id"].isna() ) pointers["parent_1_person_id"][move] = pointers["parent_2_person_id"][move] pointers["parent_2_person_id"][move] = pd.NA - relationship = _line_numbers(person["A_EXPRRP"]) - reference = np.isin(relationship, (1, 2)) - reference_households = pd.Index(household_ids[reference]) - if not reference_households.is_unique: + # The reference person is the A_EXPRRP reference person (codes 1 and 2). + # A household without one takes its is_household_head person. Where a + # household has an A_EXPRRP reference person, is_household_head must flag + # exactly that person: the reader refuses any disagreement rather than + # choosing between the two. + survey_reference = np.isin(_line_numbers(person["A_EXPRRP"]), (1, 2)) + if not pd.Index(household_ids[survey_reference]).is_unique: raise ValueError("Donor household has multiple reference persons.") - reference_rows = reference_households.get_indexer(household[HOUSEHOLD_ID_COLUMN]) + head_flags = person["is_household_head"] + if not pd.api.types.is_bool_dtype(head_flags.dtype) or head_flags.isna().any(): + raise ValueError("Donor is_household_head must be boolean with no nulls.") + head = head_flags.to_numpy(dtype=bool) + if not pd.Index(household_ids[head]).is_unique: + raise ValueError("Donor household has multiple is_household_head persons.") + surveyed = pd.Index(household_ids[survey_reference]).get_indexer(household_ids) >= 0 + if np.any(surveyed & (head != survey_reference)): + raise ValueError( + "Donor A_EXPRRP reference person and is_household_head disagree." + ) + reference = survey_reference | (head & ~surveyed) + reference_rows = pd.Index(household_ids[reference]).get_indexer( + household[HOUSEHOLD_ID_COLUMN] + ) if np.any(reference_rows < 0): raise ValueError("Donor household has no reference person.") reference_ids = person_ids[reference][reference_rows] - # The CPS relationship recode identifies the reference person's unmarried - # partner as 13 (the ACS-to-CPS crosswalk uses the same code). A_SPOUSE - # alone does not capture these cohabiting couples. - unmarried_rows = np.flatnonzero(relationship == 13) - if not pd.Index(household_ids[unmarried_rows]).is_unique: - raise ValueError("Donor household has multiple unmarried partners.") - reference_positions = np.flatnonzero(reference) - for partner_row in unmarried_rows: - reference_row = reference_positions[ - reference_households.get_indexer([household_ids[partner_row]])[0] - ] - partners = pointers["partner_person_id"] - for own_row, other_row in ( - (partner_row, reference_row), - (reference_row, partner_row), - ): - if ( - not pd.isna(partners[own_row]) - and partners[own_row] != person_ids[other_row] - ): - raise ValueError( - "Donor unmarried-partner roster conflicts with spouse pointers." - ) - partners[partner_row] = person_ids[reference_row] - partners[reference_row] = person_ids[partner_row] # The donor's three leaves are not invertible through an encoding mapping: # their sum is the reviewed donor-side stock, without inventing a split. leaves = person[["bank_account_assets", "stock_assets", "bond_assets"]].to_numpy( @@ -525,12 +620,37 @@ def _donor_extras(person: pd.DataFrame, household: pd.DataFrame) -> dict: def read_populace_us_donor(path: str | Path, *, sha256: str, size: int) -> DonorBank: """Authenticate a local donor H5, then decode its person/household inputs. - Size and SHA-256 are mandatory and verified before any HDF parsing. There - is no download or alternate donor fallback. CPS pointers are resolved by - (household id, line number), and invalid/asymmetric relations fail concept - validation. Donor-only enrichment is staged before ``decode``; all work - after that seam reads concept names. The result declares transport content - with a US donor and unchanged source-currency amounts and design weights. + Size and SHA-256 are mandatory and verified before any HDF parsing. The + file is hashed through one handle and then reopened by path, so a file + swapped or rewritten in between would be read unverified; the reader + assumes a local cache that nothing rewrites while it runs. (PyTables' core + driver can parse the hashed bytes from memory, but only under a filename + that does not exist on disk.) There is no download or alternate donor + fallback. + + CPS line pointers are resolved by (household id, line number). A pointer + that is missing or ``<= 0`` is absent: the donor the tests pin writes + ``-1`` for an absent parent or cohabiting partner and ``0`` for an absent + spouse. The partner pointer comes from ``A_SPOUSE``, or from ``PECOHAB`` + on rows whose source vintage carries it; the two must agree when both are + set. In a donor that pools vintages, cohabiting partners are therefore + recovered only for some source years (only the 2024 rows of the donor the + tests pin carry ``PECOHAB``). Relationship code 13 is never read, because + it also covers roommates and housemates. + + The household reference person is the person whose ``A_EXPRRP`` is 1 or + 2; ``A_EXPRRP`` must be present on every row. A household with no such + person takes the one person its ``is_household_head`` flag marks. Where a + household has an ``A_EXPRRP`` reference person, ``is_household_head`` + must mark exactly that person, and the reader refuses the whole donor if + any household's two disagree. Dangling pointers, spouse and cohabitant + conflicts and reference-person conflicts are refused while staging; + asymmetric, self, cyclic and other invalid relations fail concept + validation. + + Donor-only enrichment is staged before ``decode``; all work after that + seam reads concept names. The result declares transport content with a + US donor and unchanged source-currency amounts and design weights. """ path = Path(path) if isinstance(size, bool) or not isinstance(size, int) or size < 0: @@ -552,11 +672,6 @@ def read_populace_us_donor(path: str | Path, *, sha256: str, size: int) -> Donor for entity in ("person", "household") } staged = _donor_extras(raw["person"], raw["household"]) - # Use the validated roster as the input mapping's householder flag, so - # stale engine flags cannot override the survey reference-person pointer. - raw["person"]["is_household_head"] = raw["person"][PERSON_ID_COLUMN].isin( - staged["reference_person_id"] - ) tables = POLICYENGINE_US_CONCEPT_MAPPING.decode(raw) del raw for name, values in staged["pointers"].items(): diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py index 7b68fab16c..f7109f946b 100644 --- a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py @@ -236,10 +236,18 @@ def test_reader_pin_mismatch_property_precedes_hdf_parse( [ ("dangling_partner", "unknown household line"), ("dangling_parent", "unknown household line"), + ("dangling_cohabitant", "unknown household line"), + ("noninteger_pointer", "integer line numbers"), + ("infinite_pointer", "integer line numbers"), ("asymmetric_partner", "Partners must point at each other"), + ("asymmetric_cohabitant", "Partners must point at each other"), + ("spouse_cohabitant_conflict", "name different partners"), ("duplicate_line", "unique within a household"), ("multiple_reference", "multiple reference persons"), ("missing_reference", "no reference person"), + ("multiple_heads", "multiple is_household_head persons"), + ("head_disagrees", "reference person and is_household_head disagree"), + ("head_absent", "reference person and is_household_head disagree"), ("channel_mismatch", "support channels disagree"), ], ) @@ -247,20 +255,91 @@ def test_reader_refuses_invalid_rosters(tmp_path, monkeypatch, change, message): tables = _tables() person = tables["person"] mutations = { - "dangling_partner": (0, "A_SPOUSE", 99), - "dangling_parent": (2, "PEPAR1", 99), - "asymmetric_partner": (1, "A_SPOUSE", 0), - "duplicate_line": (1, "A_LINENO", 1), - "multiple_reference": (1, "A_EXPRRP", 1), - "missing_reference": (0, "A_EXPRRP", 5), - "channel_mismatch": (0, "person_support_channel", "another-channel"), + "dangling_partner": [(0, "A_SPOUSE", 99)], + "dangling_parent": [(2, "PEPAR1", 99)], + "dangling_cohabitant": [(2, "PECOHAB", 99)], + # Only whole nonpositive or missing pointers are absent; fractional + # or infinite ones still fail. + "noninteger_pointer": [(2, "PEPAR1", -1.5)], + "infinite_pointer": [(2, "PEPAR1", -np.inf)], + "asymmetric_partner": [(1, "A_SPOUSE", 0)], + "asymmetric_cohabitant": [ + (3, "A_SPOUSE", 0), + (4, "A_SPOUSE", 0), + (3, "PECOHAB", 2), + ], + "spouse_cohabitant_conflict": [(0, "PECOHAB", 3)], + "duplicate_line": [(1, "A_LINENO", 1)], + "multiple_reference": [(1, "A_EXPRRP", 1)], + # Neither an A_EXPRRP reference person nor a household head. + "missing_reference": [(0, "A_EXPRRP", 5), (0, "is_household_head", False)], + "multiple_heads": [(4, "is_household_head", True)], + "head_disagrees": [ + (0, "is_household_head", False), + (1, "is_household_head", True), + ], + "head_absent": [(0, "is_household_head", False)], + "channel_mismatch": [(0, "person_support_channel", "another-channel")], } - row, column, value = mutations[change] - person.loc[row, column] = value + for row, column, value in mutations[change]: + if isinstance(value, float) and person[column].dtype.kind == "i": + person[column] = person[column].astype(np.float64) + person.loc[row, column] = value with pytest.raises(ValueError, match=message): _fake_donor(monkeypatch, tmp_path, tables) +@pytest.mark.parametrize( + "flags", + [ + pd.Series([1, 0, 0, 1, 0, 0]), + pd.Series([True, False, None] * 2, dtype="boolean"), + ], + ids=["integer", "nullable"], +) +def test_reader_requires_a_boolean_household_head_flag(tmp_path, monkeypatch, flags): + tables = _tables() + tables["person"]["is_household_head"] = flags + with pytest.raises(ValueError, match="is_household_head must be boolean"): + _fake_donor(monkeypatch, tmp_path, tables) + + +@pytest.mark.parametrize("absent", [-1, 0, np.nan]) +def test_reader_reads_nonpositive_or_missing_pointers_as_absent( + tmp_path, monkeypatch, absent +): + tables = _tables() + person = tables["person"] + for column in ("A_SPOUSE", "PEPAR1", "PEPAR2", "PECOHAB"): + values = person[column].astype(np.float64) + person[column] = values.where(values > 0, absent) + donor = _fake_donor(monkeypatch, tmp_path, tables) + expected = _fixture()["expected"] + for column in ("partner_person_id", "parent_1_person_id", "parent_2_person_id"): + pd.testing.assert_series_equal( + donor.tables["person"][column], + pd.Series(pd.array(expected[column], dtype="Int64"), name=column), + ) + + +def test_reader_takes_household_head_when_no_a_exprrp_reference_person( + tmp_path, monkeypatch +): + tables = _tables() + person = tables["person"] + # Household 202 has no A_EXPRRP reference person, and its head flag marks + # line 2 rather than line 1, so only the fallback can name that person. + person.loc[3, "A_EXPRRP"] = 3 + person.loc[3, "is_household_head"] = False + person.loc[4, "is_household_head"] = True + donor = _fake_donor(monkeypatch, tmp_path, tables) + assert donor.tables["household"]["reference_person_id"].tolist() == [ + person.loc[0, "person_id"], + person.loc[4, "person_id"], + ] + assert validate_concept_tables(donor.tables) == () + + @READER_PROPERTY @given( offset=st.integers(min_value=2**53 + 1, max_value=2**63 - 7), @@ -295,7 +374,7 @@ def test_pointer_recode_is_exact_bijective_and_partner_symmetric( parents = [ lookup[(row.person_household_id, line)] for line in (row.PEPAR1, row.PEPAR2) - if line + if line > 0 ] for position, column in enumerate(("parent_1_person_id", "parent_2_person_id")): observed = getattr(recovered, column) @@ -381,35 +460,62 @@ def test_reader_refuses_partial_source_lineage(tmp_path, monkeypatch, missing): _fake_donor(monkeypatch, tmp_path, tables) -def test_reader_recodes_unmarried_reference_partner_symmetrically( - tmp_path, monkeypatch +@pytest.mark.parametrize("housemates", [1, 2]) +def test_reader_never_reads_relationship_code_13_as_a_partner( + tmp_path, monkeypatch, housemates ): tables = _tables() person = tables["person"] - person.loc[[0, 1], "A_SPOUSE"] = 0 - person.loc[1, "A_EXPRRP"] = 13 + # Code 13 also covers roommates and housemates. Household 202's reference + # person shares the home with one or two code-13 housemates and no spouse + # or cohabiting-partner pointer, so nobody in it has a partner. + person.loc[[3, 4], "A_SPOUSE"] = 0 + person.loc[4, "A_EXPRRP"] = 13 + if housemates == 2: + person.loc[5, ["A_EXPRRP", "PEPAR2"]] = [13, -1] donor = _fake_donor(monkeypatch, tmp_path, tables) - pointers = donor.tables["person"]["partner_person_id"] - assert pointers.iloc[0] == person.loc[1, "person_id"] - assert pointers.iloc[1] == person.loc[0, "person_id"] + partners = donor.tables["person"]["partner_person_id"] + assert partners.iloc[3:].isna().all() + assert partners.iloc[:2].tolist() == [ + person.loc[1, "person_id"], + person.loc[0, "person_id"], + ] assert validate_concept_tables(donor.tables) == () -@pytest.mark.parametrize("conflict", ["multiple", "contradictory_spouse"]) -def test_reader_refuses_conflicting_unmarried_reference_partners( - tmp_path, monkeypatch, conflict +def test_reader_recodes_cohabiting_partner_pointers_symmetrically( + tmp_path, monkeypatch ): tables = _tables() person = tables["person"] - person.loc[[0, 1], "A_SPOUSE"] = 0 - person.loc[1, "A_EXPRRP"] = 13 - if conflict == "multiple": - person.loc[2, "A_EXPRRP"] = 13 + # Household 202's couple cohabits: no spouse pointer, a PECOHAB pair. + person.loc[[3, 4], "A_SPOUSE"] = 0 + person.loc[[3, 4], "PECOHAB"] = [2, 1] + donor = _fake_donor(monkeypatch, tmp_path, tables) + partners = donor.tables["person"]["partner_person_id"] + assert partners.iloc[3] == person.loc[4, "person_id"] + assert partners.iloc[4] == person.loc[3, "person_id"] + assert validate_concept_tables(donor.tables) == () + + +@pytest.mark.parametrize("cohabitation", ["agrees_with_spouse", "column_absent"]) +def test_reader_spouse_pointers_hold_with_or_without_pecohab( + tmp_path, monkeypatch, cohabitation +): + tables = _tables() + if cohabitation == "agrees_with_spouse": + # PECOHAB may repeat a spouse pointer; only different partners fail. + tables["person"].loc[[0, 1], "PECOHAB"] = [2, 1] else: - person.loc[0, "A_SPOUSE"] = 3 - person.loc[2, "A_SPOUSE"] = 1 - with pytest.raises(ValueError): - _fake_donor(monkeypatch, tmp_path, tables) + tables["person"] = tables["person"].drop(columns="PECOHAB") + donor = _fake_donor(monkeypatch, tmp_path, tables) + pd.testing.assert_series_equal( + donor.tables["person"]["partner_person_id"], + pd.Series( + pd.array(_fixture()["expected"]["partner_person_id"], dtype="Int64"), + name="partner_person_id", + ), + ) def test_cached_real_donor_size_only_when_present(): @@ -420,6 +526,26 @@ def test_cached_real_donor_size_only_when_present(): assert path.stat().st_size == pin["size"] +def test_pinned_real_donor_reaches_the_reference_person_ruling_when_cached(): + # About 10 s: hashes the 463 MB pin, reads both tables and stages every + # line pointer, which must all resolve, before the reference-person + # check. #1130 is waiting on a ruling: in 113 pinned households (all + # 2022/2023 rows, whose A_EXPRRP the build derived from line 1) the + # A_EXPRRP reference person is not the is_household_head (P_SEQ == 1) + # person, so the reader refuses the file. Once that is ruled, assert the + # bank instead; one source household's two copies also hold a + # parent-pointer cycle that concept validation refuses. + pin = _fixture()["real_donor_pin"] + path = Path.home() / pin["cache_relative_path"] + if not path.is_file(): + pytest.skip("Pinned real donor is not locally cached") + with pytest.raises( + ValueError, + match=r"^Donor A_EXPRRP reference person and is_household_head disagree\.$", + ): + read_populace_us_donor(path, sha256=pin["sha256"], size=pin["size"]) + + @PROPERTY @given( ids=st.lists( diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py index dccca8aa37..de92a71478 100644 --- a/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py @@ -4,9 +4,10 @@ import numpy as np import pytest -from hypothesis import given, settings +from hypothesis import example, given, settings from hypothesis import strategies as st +from microcosm.frame import transport from microcosm.frame.transport import quantile_map @@ -19,6 +20,18 @@ def _bands(): ] +def _one_entity(values, total): + # One entity whose total maps to the given value (the band's midpoint). + bands = [{"lower": 0.0, "upper": 2 * abs(total), "share": 1}] + return quantile_map( + values, + np.ones(len(values)), + {"positive": bands, "negative": bands}, + interpolation="uniform", + group={"entity_ids": np.zeros(len(values), dtype=int)}, + ) + + _samples = st.lists( st.tuples(st.integers(-100, 100), st.integers(1, 40), st.integers(0, 2)), min_size=1, @@ -267,6 +280,10 @@ def test_permuting_members_preserves_mapped_values(samples): st.lists(st.integers(0, 50), min_size=1, max_size=5), min_size=1, max_size=12 ) ) +# Members [1, 2] whose total maps to a total the largest member alone cannot +# reach; the second is the counterexample #1130's engine-free CI job found. +@example(entities=[[1], [1], [1, 2]]) +@example(entities=[[1], [1], [1], [1], [1], [1], [1], [4], [1, 2]]) @settings(max_examples=100, deadline=None) def test_aggregate_pro_rata_totals_and_idempotence(entities): values = np.array([value for members in entities for value in members], dtype=float) @@ -472,6 +489,132 @@ def test_pro_rata_mapping_avoids_overflow_in_the_common_scale_factor(): np.testing.assert_array_equal(mapped, [2.5e307, 2.5e307]) +def test_pro_rata_correction_beyond_float64_rounding_is_refused(monkeypatch): + # No correctly rounded entity total has been found to move the anchor past + # the bound, so a mis-rounded total stands in for a pro rata step that + # loses more than rounding. Correcting the anchor still restores the + # exact total, so only the correction bound refuses the result. + members = [1.0, 3.0] + + def misrounded(items): + items = list(items) + total = fsum(items) + return total * (1 + 1e-9) if items == members else total + + monkeypatch.setattr(transport, "fsum", misrounded) + with pytest.raises(ValueError, match="correction exceeds float64 rounding"): + quantile_map( + members, + [1, 1], + [{"lower": 0, "upper": 120, "share": 1}], + interpolation="uniform", + group={"entity_ids": [0, 0]}, + ) + + +def test_pro_rata_moves_one_other_member_when_the_largest_cannot_reach_the_total(): + # [1, 2] scaled to a total of 0.21: every value of the larger member + # leaves the sum halfway between 0.21 and a float64 neighbour, and + # rounding to even never lands on 0.21, so the smaller member moves by + # one unit in the last place. + values = np.array([1.0, 2.0]) + mapped = quantile_map( + values, + [1, 1], + [{"lower": 0, "upper": 0.42, "share": 1}], + interpolation="uniform", + group={"entity_ids": [0, 0]}, + ) + pro_rata = values / 3 * 0.21 + assert fsum(mapped) == 0.21 + assert mapped[0] == np.nextafter(pro_rata[0], np.inf) + assert mapped[1] == pro_rata[1] + + +_amounts = st.one_of( + st.integers(1, 50).map(float), st.integers(1, 10**6).map(lambda cents: cents / 100) +) + + +@given( + # Whole amounts and amounts in cents. Entities of one to three members + # produce the halfway ties most often; larger ones cover household sizes. + st.lists( + st.one_of( + st.lists(_amounts, min_size=1, max_size=3), + st.lists(_amounts, min_size=4, max_size=16), + ), + min_size=1, + max_size=20, + ), + st.sampled_from([1.0, -1.0]), +) +@example(entities=[[1.0], [1.0], [1.0, 2.0]], sign=1.0) +@example(entities=[[1.0], [1.0], [1.0, 2.0]], sign=-1.0) +@example( + entities=[[1.0], [2.0], [3.0], [1.62, 1.53], [5.0], [6.0], [7.0], [8.0], [9.0]], + sign=1.0, +) +@settings(max_examples=100, deadline=None) +def test_same_sign_members_always_retain_exact_entity_totals(entities, sign): + values = sign * np.array([value for members in entities for value in members]) + entity_ids = np.array( + [entity for entity, members in enumerate(entities) for _ in members] + ) + mapped = quantile_map( + values, + np.ones(len(values)), + _bands(), + interpolation="uniform", + group={"entity_ids": entity_ids}, + ) + totals = np.array([fsum(values[entity_ids == e]) for e in range(len(entities))]) + expected = quantile_map( + totals, np.ones(len(entities)), _bands(), interpolation="uniform" + ) + for entity, total in enumerate(totals): + selected = entity_ids == entity + assert fsum(mapped[selected]) == expected[entity] + pro_rata = values[selected] / total * expected[entity] + moved = np.abs(mapped[selected] - pro_rata) / np.spacing(np.abs(pro_rata)) + # Only the largest member absorbs more than one unit in the last place. + assert (np.delete(moved, np.argmax(np.abs(pro_rata))) <= 1).all() + + +def test_pro_rata_reaches_a_power_of_two_total_far_from_the_residual_step(): + # Thirteen members scaled to 8192: the total's spacing halves below it, + # and the residual step from the rounded sum misses by more than the two + # one-unit steps can recover. The largest member alone can still reach + # the total, so every other member keeps its exact pro rata value. + values = [575.99, 442.46, 26.4, 365.04, 974.36, 860.0, 783.83] + values += [715.86, 762.91, 803.26, 810.82, 704.65, 512.47] + mapped = _one_entity(values, 8192.0) + pro_rata = np.array(values) / fsum(values) * 8192.0 + largest = np.argmax(values) + assert fsum(mapped) == 8192.0 + np.testing.assert_array_equal( + np.delete(mapped, largest), np.delete(pro_rata, largest) + ) + + +@pytest.mark.parametrize("sign", [1.0, -1.0]) +def test_pro_rata_moves_the_smallest_member_when_the_tie_sits_a_binade_up(sign): + # [7, 3.57] scaled to 1.51: the larger member falls just below 1, but the + # anchor the total needs sits halfway between two values above 1, where + # the other member's spacing is no finer than the larger member's own. + values = sign * np.array([7.0, 3.57]) + mapped = _one_entity(values, sign * 1.51) + pro_rata = values / fsum(values) * (sign * 1.51) + assert fsum(mapped) == sign * 1.51 + assert (np.abs(mapped - pro_rata) <= np.spacing(np.abs(pro_rata))).all() + + +def test_mixed_sign_members_mapped_once_map_to_themselves_again(): + first = _one_entity([30.59, -90.38], -763.16) + assert fsum(first) == -763.16 + np.testing.assert_array_equal(_one_entity(first, -763.16), first) + + def test_cancelling_members_refuse_an_unrepresentable_aggregate_total(): bands = [{"lower": 10, "upper": 30, "share": 1}] with pytest.raises(ValueError, match="exact totals"): diff --git a/packages/microcosm-frame/tests/fixtures/transport/synthetic_donor.json b/packages/microcosm-frame/tests/fixtures/transport/synthetic_donor.json index c19959aeab..8142107314 100644 --- a/packages/microcosm-frame/tests/fixtures/transport/synthetic_donor.json +++ b/packages/microcosm-frame/tests/fixtures/transport/synthetic_donor.json @@ -14,9 +14,11 @@ "person_support_channel": "asec", "A_LINENO": 1, "A_SPOUSE": 2, - "PEPAR1": 0, - "PEPAR2": 0, + "PEPAR1": -1, + "PEPAR2": -1, + "PECOHAB": -1, "A_EXPRRP": 1, + "is_household_head": true, "age": 40, "is_female": false, "employment_income_before_lsr": 0, @@ -37,9 +39,11 @@ "person_support_channel": "asec", "A_LINENO": 2, "A_SPOUSE": 1, - "PEPAR1": 0, - "PEPAR2": 0, + "PEPAR1": -1, + "PEPAR2": -1, + "PECOHAB": -1, "A_EXPRRP": 4, + "is_household_head": false, "age": 38, "is_female": true, "employment_income_before_lsr": 10, @@ -62,7 +66,9 @@ "A_SPOUSE": 0, "PEPAR1": 1, "PEPAR2": 2, + "PECOHAB": -1, "A_EXPRRP": 5, + "is_household_head": false, "age": 10, "is_female": false, "employment_income_before_lsr": 20, @@ -83,9 +89,11 @@ "person_support_channel": "puf_tax_detail", "A_LINENO": 1, "A_SPOUSE": 2, - "PEPAR1": 0, - "PEPAR2": 0, + "PEPAR1": -1, + "PEPAR2": -1, + "PECOHAB": null, "A_EXPRRP": 1, + "is_household_head": true, "age": 45, "is_female": false, "employment_income_before_lsr": 30, @@ -106,9 +114,11 @@ "person_support_channel": "puf_tax_detail", "A_LINENO": 2, "A_SPOUSE": 1, - "PEPAR1": 0, - "PEPAR2": 0, + "PEPAR1": -1, + "PEPAR2": -1, + "PECOHAB": null, "A_EXPRRP": 4, + "is_household_head": false, "age": 42, "is_female": true, "employment_income_before_lsr": 40, @@ -129,9 +139,11 @@ "person_support_channel": "puf_tax_detail", "A_LINENO": 3, "A_SPOUSE": 0, - "PEPAR1": 0, + "PEPAR1": -1, "PEPAR2": 2, + "PECOHAB": null, "A_EXPRRP": 5, + "is_household_head": false, "age": 12, "is_female": true, "employment_income_before_lsr": 50, From 9df282d3b5053fb81ca3df0072b7fb01d8367590 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Wed, 7 Oct 2026 15:08:28 -0400 Subject: [PATCH 3/6] Apply NZ hub head and parent-cycle rulings to donor reader (#1130) Require exactly one is_household_head person per donor household and use that person as its reference. Remove A_EXPRRP selection, fallback and disagreement refusal; document the P_SEQ and Census-recode evidence. Repair reciprocal two-person parent cycles by dropping the edge naming a younger parent, or both edges for equal ages or an absent age column. Check the original graph for longer cycles before mutation, preserve all other concept refusals, compact surviving parents and expose the dropped directed-edge count as DonorBank.dropped_parent_cycle_edges. Update regression tests, the cached full-reader smoke, ADR and changelog. The pinned donor now returns 57,240 households, 166,321 persons and two dropped parent edges with no concept violations. Targeted donor tests: 84 passed, no skips. Targeted quantile tests: 60 passed. Ruff and CI test layout/graph acceptance verification pass. Co-Authored-By: GPT-6.1 Sol --- changelog.d/nz-transport-operators.added.md | 12 +- docs/concept-schema-transport-adr.md | 42 ++-- .../src/microcosm/frame/transport.py | 116 ++++++--- .../shared/test_transport_donor.py | 225 +++++++++++++++--- 4 files changed, 318 insertions(+), 77 deletions(-) diff --git a/changelog.d/nz-transport-operators.added.md b/changelog.d/nz-transport-operators.added.md index 3bb786ed7f..71765ccb2d 100644 --- a/changelog.d/nz-transport-operators.added.md +++ b/changelog.d/nz-transport-operators.added.md @@ -1,8 +1,10 @@ Add engine-free transport operators: a reader for local populace-US donor files pinned by size and SHA-256, stable source-ID seeds, one-step currency conversion, and weighted signed quantile maps with component and -entity-aggregate support. Run on the donor its tests pin, the reader resolves -every line pointer but refuses the file: in 113 of its 57,240 households the -`A_EXPRRP` reference person is not the `is_household_head` person. Once that -is resolved, concept validation would also refuse two households whose -parent pointers form a cycle. +entity-aggregate support. The reader requires exactly one +`is_household_head` reference person per household. For a two-person parent +cycle, it drops the edge naming a younger parent, or both edges when ages are +equal or the donor has no age column, and returns the dropped-edge count. +Longer cycles and other inconsistent relationships are refused. +The full reader smoke test loads the pinned donor's 166,321 persons in +57,240 households, dropping two parent edges. diff --git a/docs/concept-schema-transport-adr.md b/docs/concept-schema-transport-adr.md index fd22a2e667..cd5ebb9436 100644 --- a/docs/concept-schema-transport-adr.md +++ b/docs/concept-schema-transport-adr.md @@ -3,16 +3,15 @@ **Status:** accepted for the concept schema, the per-engine mappings and their coverage reports. `microcosm.frame.transport` adds a reader for local populace-US donor files pinned by size and SHA-256, source-ID seeds, a -currency bridge and weighted quantile mapping. The reader does not yet load -the donor its tests pin. Run on that file (166,321 persons, 57,240 -households), it resolves every line pointer, then refuses the whole file: in -113 households the `A_EXPRRP` reference person is not the `is_household_head` -person, and the reader will not choose between them. All 113 are 2022 or 2023 -rows, whose `A_EXPRRP` matches the build's line-1 derivation rather than the -Census recode. Once that is resolved, concept validation would also refuse -two other households, copies of one source household in which two people name -each other as parents. The composed transport pipeline below remains a -design, and no build has migrated. **Date:** 2026-09-28. **Source:** Max's 27 September rulings on +currency bridge and weighted quantile mapping. The reader selects the donor's +`is_household_head` as the reference person and repairs two-person parent +cycles under the age rule below, following the NZ hub's 2026-10-07 defaults +recorded for D1 method card d1043. The full reader smoke test returns a valid +bank from its pinned donor: 166,321 persons, 57,240 households, 57,240 heads, +74,572 partner pointers, 62,096 `parent_1_person_id` pointers and 43,731 +`parent_2_person_id` pointers, after dropping two parent edges. Concept +validation reports no violations. The composed transport pipeline below +remains a design, and no build has migrated. **Date:** 2026-09-28. **Source:** Max's 27 September rulings on international populations and law-anchored concepts. **Amended 2026-10-07:** schema version 2 adds `liquid_financial_assets`, the first stock concept, for New Zealand's Accommodation Supplement cash-asset test. Only the Axiom New @@ -296,12 +295,23 @@ These hold for every input and are tested (Hypothesis properties unless noted): - `microcosm.frame.transport` supplies a reader for pinned local donor files and the income transformation operators (step 3). The reader's bank declares US donor provenance and source-currency amounts; it does not - declare destination residents. It produces no bank from the donor its - tests pin yet: that file stops at the reference-person disagreement in the - status line above, and two households with a parent-pointer cycle would - then fail concept validation. The consuming builder supplies the - destination frame declaration after transport. Transport still needs a - unit-construction step that builds engine + declare destination residents. Each household must have exactly one person + whose boolean `is_household_head` flag is true; this supplies its + `reference_person_id`. `A_EXPRRP` is not used to select it. On the pinned + donor, `is_household_head` equals `P_SEQ == 1` on every row and is the US + engine's own head. It agrees with the Census `A_EXPRRP` reference person in + all 19,753 households from 2024, including 222 where that person is not + line 1. The 2022/2023 `A_EXPRRP` values instead come from the build's line-1 + derivation. +- If two people name each other as parents, the reader drops the edge whose + named parent is younger than the child and keeps the other edge. If their + ages are equal or the donor has no age column, it drops both edges. A + present age column must satisfy the existing concept contract: whole, + nonmissing ages in completed years, from 0 through 130. + `DonorBank.dropped_parent_cycle_edges` exposes the number of dropped edges. + Longer cycles and other inconsistent relationships still fail validation. +- The consuming builder supplies the destination frame declaration after + transport. Transport still needs a unit-construction step that builds engine units from the concept pointers (until it exists, 13 of New Zealand's 80 and 46 of Belgium's 105 module bindings, and the take-up and housing bindings on US and UK group entities, stay deferred), and a country pack with geography, diff --git a/packages/microcosm-frame/src/microcosm/frame/transport.py b/packages/microcosm-frame/src/microcosm/frame/transport.py index 7e4dea2cd9..66a732bbf6 100644 --- a/packages/microcosm-frame/src/microcosm/frame/transport.py +++ b/packages/microcosm-frame/src/microcosm/frame/transport.py @@ -405,6 +405,8 @@ class DonorBank: ``currency`` describes source amounts; the destination producer supplies its own frame declaration after transport. No destination country or year is inferred by this reader. + ``dropped_parent_cycle_edges`` counts directed parent pointers removed by + the reader's age-based two-person-cycle rule, including support copies. """ tables: Mapping[str, pd.DataFrame] @@ -414,6 +416,7 @@ class DonorBank: content_basis: ContentBasis = ContentBasis.TRANSPORT donor_country: str = "us" currency: str = "USD" + dropped_parent_cycle_edges: int = 0 def _line_numbers(values: pd.Series, *, pointer: bool = False) -> np.ndarray: @@ -456,7 +459,6 @@ def _donor_extras(person: pd.DataFrame, household: pd.DataFrame) -> dict: "A_SPOUSE", "PEPAR1", "PEPAR2", - "A_EXPRRP", "is_household_head", "bank_account_assets", "stock_assets", @@ -520,32 +522,20 @@ def pointer(present: np.ndarray, rows: np.ndarray) -> pd.arrays.IntegerArray: ) pointers["parent_1_person_id"][move] = pointers["parent_2_person_id"][move] pointers["parent_2_person_id"][move] = pd.NA - # The reference person is the A_EXPRRP reference person (codes 1 and 2). - # A household without one takes its is_household_head person. Where a - # household has an A_EXPRRP reference person, is_household_head must flag - # exactly that person: the reader refuses any disagreement rather than - # choosing between the two. - survey_reference = np.isin(_line_numbers(person["A_EXPRRP"]), (1, 2)) - if not pd.Index(household_ids[survey_reference]).is_unique: - raise ValueError("Donor household has multiple reference persons.") + # The donor's US-engine head (P_SEQ == 1) is the reference person in + # every vintage. A_EXPRRP is not read: older rows derive it from line 1. head_flags = person["is_household_head"] if not pd.api.types.is_bool_dtype(head_flags.dtype) or head_flags.isna().any(): raise ValueError("Donor is_household_head must be boolean with no nulls.") head = head_flags.to_numpy(dtype=bool) if not pd.Index(household_ids[head]).is_unique: raise ValueError("Donor household has multiple is_household_head persons.") - surveyed = pd.Index(household_ids[survey_reference]).get_indexer(household_ids) >= 0 - if np.any(surveyed & (head != survey_reference)): - raise ValueError( - "Donor A_EXPRRP reference person and is_household_head disagree." - ) - reference = survey_reference | (head & ~surveyed) - reference_rows = pd.Index(household_ids[reference]).get_indexer( + reference_rows = pd.Index(household_ids[head]).get_indexer( household[HOUSEHOLD_ID_COLUMN] ) if np.any(reference_rows < 0): - raise ValueError("Donor household has no reference person.") - reference_ids = person_ids[reference][reference_rows] + raise ValueError("Donor household has no is_household_head person.") + reference_ids = person_ids[head][reference_rows] # The donor's three leaves are not invertible through an encoding mapping: # their sum is the reviewed donor-side stock, without inventing a split. leaves = person[["bank_account_assets", "stock_assets", "bond_assets"]].to_numpy( @@ -617,6 +607,66 @@ def pointer(present: np.ndarray, rows: np.ndarray) -> pd.arrays.IntegerArray: } +def _drop_two_person_parent_cycles(person: pd.DataFrame) -> int: + """Repair reciprocal parent edges only after other concepts validate. + + Check the original graph before mutation: a reciprocal edge may also + belong to a longer cycle, which must still refuse. Removing its direct + reverse edge for a reachability search exposes any such longer cycle. + Searches stay within a household and run only for reciprocal edges. + """ + columns = ("parent_1_person_id", "parent_2_person_id") + ids = pd.Index(person[PERSON_ID_COLUMN]) + parents = np.full((len(person), len(columns)), -1, dtype=np.int64) + for slot, column in enumerate(columns): + present = person[column].notna().to_numpy() + parents[present, slot] = ids.get_indexer( + person[column][present].to_numpy(dtype=np.int64) + ) + rows = np.arange(len(person)) + reciprocal = np.zeros(parents.shape, dtype=bool) + for slot in range(len(columns)): + present = parents[:, slot] >= 0 + reciprocal[present, slot] = np.any( + parents[parents[present, slot]] == rows[present, None], axis=1 + ) + edges = np.argwhere(reciprocal) + if not len(edges): + return 0 + for child, slot in edges: + parent = parents[child, slot] + pending = [parent] + visited = {parent} + while pending: + current = pending.pop() + for target in parents[current]: + if target < 0 or (current == parent and target == child): + continue + if target == child: + # Leave the original violations for the caller to report. + return 0 + if target not in visited: + visited.add(target) + pending.append(target) + ages = ( + person["age"].to_numpy(dtype=np.float64, na_value=np.nan) + if "age" in person + else np.full(len(person), np.nan) + ) + children, slots = edges.T + targets = parents[children, slots] + known = np.isfinite(ages[children]) & np.isfinite(ages[targets]) + drop = ~known | (ages[targets] <= ages[children]) + for slot, column in enumerate(columns): + person.iloc[ + children[drop & (slots == slot)], person.columns.get_loc(column) + ] = pd.NA + move = person[columns[0]].isna() & person[columns[1]].notna() + person.loc[move, columns[0]] = person.loc[move, columns[1]] + person.loc[move, columns[1]] = pd.NA + return int(drop.sum()) + + def read_populace_us_donor(path: str | Path, *, sha256: str, size: int) -> DonorBank: """Authenticate a local donor H5, then decode its person/household inputs. @@ -638,15 +688,21 @@ def read_populace_us_donor(path: str | Path, *, sha256: str, size: int) -> Donor tests pin carry ``PECOHAB``). Relationship code 13 is never read, because it also covers roommates and housemates. - The household reference person is the person whose ``A_EXPRRP`` is 1 or - 2; ``A_EXPRRP`` must be present on every row. A household with no such - person takes the one person its ``is_household_head`` flag marks. Where a - household has an ``A_EXPRRP`` reference person, ``is_household_head`` - must mark exactly that person, and the reader refuses the whole donor if - any household's two disagree. Dangling pointers, spouse and cohabitant - conflicts and reference-person conflicts are refused while staging; - asymmetric, self, cyclic and other invalid relations fail concept - validation. + The reference person is the donor's ``is_household_head`` person in + every household; exactly one boolean head flag is required. The US build + sets it from ``P_SEQ == 1`` (``relationship_inputs.py``). On all 19,753 + households from 2024 in the pinned donor, that head agrees with the Census + reference person, including 222 whose reference person is not line 1. + The 2022/2023 ``A_EXPRRP`` recode instead derives from line 1 + (``asec_pool.py``), so ``A_EXPRRP`` is not used. + + When two people name each other as parents, drop the edge naming the + younger parent and keep the other. Equal ages or an absent age column + drop both edges. The bank's ``dropped_parent_cycle_edges`` counts removed + directed pointers. A remaining sole parent is compacted into parent 1. + Longer cycles, including ones sharing reciprocal edges, and every other + concept inconsistency still refuse. A present age column must satisfy + the concept schema's whole, nonmissing-age contract. Donor-only enrichment is staged before ``decode``; all work after that seam reads concept names. The result declares transport content with a @@ -679,6 +735,11 @@ def read_populace_us_donor(path: str | Path, *, sha256: str, size: int) -> Donor tables["person"]["liquid_financial_assets"] = staged["liquid_financial_assets"] tables["household"]["reference_person_id"] = staged["reference_person_id"] violations = validate_concept_tables(tables) + dropped = 0 + if violations and all(v.code == "parent_cycle" for v in violations): + dropped = _drop_two_person_parent_cycles(tables["person"]) + if dropped: + violations = validate_concept_tables(tables) if violations: details = "; ".join(f"{v.entity}.{v.column}: {v.message}" for v in violations) raise ValueError(f"Invalid donor concepts: {details}") @@ -687,6 +748,7 @@ def read_populace_us_donor(path: str | Path, *, sha256: str, size: int) -> Donor staged["weights"], staged["support_strata"], staged["source_person_ids"], + dropped_parent_cycle_edges=dropped, ) diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py index f7109f946b..3f89358c66 100644 --- a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py @@ -130,6 +130,7 @@ def spy_decode(mapping, tables): assert donor.content_basis is ContentBasis.TRANSPORT assert donor.donor_country == "us" assert donor.currency == "USD" + assert donor.dropped_parent_cycle_edges == 0 assert validate_concept_tables(donor.tables) == () assert "A_LINENO" not in person and "bank_account_assets" not in person assert "household_weight" not in donor.tables["household"] @@ -243,11 +244,9 @@ def test_reader_pin_mismatch_property_precedes_hdf_parse( ("asymmetric_cohabitant", "Partners must point at each other"), ("spouse_cohabitant_conflict", "name different partners"), ("duplicate_line", "unique within a household"), - ("multiple_reference", "multiple reference persons"), - ("missing_reference", "no reference person"), ("multiple_heads", "multiple is_household_head persons"), - ("head_disagrees", "reference person and is_household_head disagree"), - ("head_absent", "reference person and is_household_head disagree"), + ("missing_head_with_reference", "no is_household_head"), + ("missing_head_without_reference", "no is_household_head"), ("channel_mismatch", "support channels disagree"), ], ) @@ -270,15 +269,12 @@ def test_reader_refuses_invalid_rosters(tmp_path, monkeypatch, change, message): ], "spouse_cohabitant_conflict": [(0, "PECOHAB", 3)], "duplicate_line": [(1, "A_LINENO", 1)], - "multiple_reference": [(1, "A_EXPRRP", 1)], - # Neither an A_EXPRRP reference person nor a household head. - "missing_reference": [(0, "A_EXPRRP", 5), (0, "is_household_head", False)], "multiple_heads": [(4, "is_household_head", True)], - "head_disagrees": [ + "missing_head_with_reference": [(0, "is_household_head", False)], + "missing_head_without_reference": [ + (0, "A_EXPRRP", 5), (0, "is_household_head", False), - (1, "is_household_head", True), ], - "head_absent": [(0, "is_household_head", False)], "channel_mismatch": [(0, "person_support_channel", "another-channel")], } for row, column, value in mutations[change]: @@ -322,16 +318,22 @@ def test_reader_reads_nonpositive_or_missing_pointers_as_absent( ) -def test_reader_takes_household_head_when_no_a_exprrp_reference_person( - tmp_path, monkeypatch +@pytest.mark.parametrize("recode", ["different", "absent", "duplicate", "no_column"]) +def test_reader_uses_household_head_regardless_of_a_exprrp( + tmp_path, monkeypatch, recode ): tables = _tables() person = tables["person"] - # Household 202 has no A_EXPRRP reference person, and its head flag marks - # line 2 rather than line 1, so only the fallback can name that person. - person.loc[3, "A_EXPRRP"] = 3 + # Household 202's flagged head is line 2. The survey recode may instead + # name line 1, name nobody or name multiple persons; none chooses the head. person.loc[3, "is_household_head"] = False person.loc[4, "is_household_head"] = True + if recode == "absent": + person.loc[3, "A_EXPRRP"] = 3 + elif recode == "duplicate": + person.loc[4, "A_EXPRRP"] = 2 + elif recode == "no_column": + tables["person"] = person.drop(columns="A_EXPRRP") donor = _fake_donor(monkeypatch, tmp_path, tables) assert donor.tables["household"]["reference_person_id"].tolist() == [ person.loc[0, "person_id"], @@ -340,6 +342,169 @@ def test_reader_takes_household_head_when_no_a_exprrp_reference_person( assert validate_concept_tables(donor.tables) == () +@pytest.mark.parametrize("younger_row", [4, 5]) +@pytest.mark.parametrize("parent_column", ["PEPAR1", "PEPAR2"]) +def test_reader_drops_only_the_younger_parent_edge_of_a_two_person_cycle( + tmp_path, monkeypatch, younger_row, parent_column +): + tables = _tables() + person = tables["person"] + person.loc[[3, 4], "A_SPOUSE"] = 0 + person.loc[[4, 5], ["PEPAR1", "PEPAR2"]] = -1 + person.loc[[4, 5], parent_column] = [3, 2] + older_row = 9 - younger_row + person.loc[younger_row, "age"] = 31 + person.loc[older_row, "age"] = 51 + donor = _fake_donor(monkeypatch, tmp_path, tables) + parents = donor.tables["person"]["parent_1_person_id"] + assert parents.iloc[younger_row] == person.loc[older_row, "person_id"] + assert pd.isna(parents.iloc[older_row]) + assert donor.tables["person"]["parent_2_person_id"].iloc[[4, 5]].isna().all() + assert donor.dropped_parent_cycle_edges == 1 + assert validate_concept_tables(donor.tables) == () + + +@pytest.mark.parametrize("ages", ["equal", "no_column"]) +def test_reader_drops_both_two_person_cycle_edges_when_ages_cannot_order_them( + tmp_path, monkeypatch, ages +): + tables = _tables() + person = tables["person"] + person.loc[[3, 4], "A_SPOUSE"] = 0 + person.loc[[4, 5], ["PEPAR1", "PEPAR2"]] = -1 + # One edge starts in parent 2, so either slot participates in the rule. + person.loc[4, "PEPAR1"] = 3 + person.loc[5, "PEPAR2"] = 2 + if ages == "equal": + person.loc[[4, 5], "age"] = 31 + else: + tables["person"] = person.drop(columns="age") + donor = _fake_donor(monkeypatch, tmp_path, tables) + for column in ("parent_1_person_id", "parent_2_person_id"): + assert donor.tables["person"][column].iloc[[4, 5]].isna().all() + assert donor.dropped_parent_cycle_edges == 2 + assert validate_concept_tables(donor.tables) == () + + +@pytest.mark.parametrize("missing_rows", [(4,), (5,), (4, 5)]) +def test_reader_cycle_rule_preserves_the_nonmissing_age_contract( + tmp_path, monkeypatch, missing_rows +): + tables = _tables() + person = tables["person"] + person.loc[[3, 4], "A_SPOUSE"] = 0 + person.loc[[4, 5], ["PEPAR1", "PEPAR2"]] = -1 + person.loc[4, "PEPAR1"] = 3 + person.loc[5, "PEPAR2"] = 2 + # An absent age column is permitted, but null cells in a present column + # violate the concept's whole, nonmissing-age contract before repair. + person["age"] = person["age"].astype(np.float64) + person.loc[list(missing_rows), "age"] = np.nan + with pytest.raises( + ValueError, match="fact:person.age decodes only from whole numbers" + ): + _fake_donor(monkeypatch, tmp_path, tables) + + +def test_reader_compacts_the_other_parent_after_dropping_a_cycle_edge( + tmp_path, monkeypatch +): + tables = _tables() + person = tables["person"] + person.loc[[3, 4], "A_SPOUSE"] = 0 + person.loc[4, ["age", "PEPAR1", "PEPAR2"]] = [31, 3, -1] + person.loc[5, ["age", "PEPAR1", "PEPAR2"]] = [51, 2, 1] + donor = _fake_donor(monkeypatch, tmp_path, tables) + parents = donor.tables["person"] + assert parents["parent_1_person_id"].iloc[4] == person.loc[5, "person_id"] + assert parents["parent_1_person_id"].iloc[5] == person.loc[3, "person_id"] + assert parents["parent_2_person_id"].iloc[[4, 5]].isna().all() + assert donor.dropped_parent_cycle_edges == 1 + assert validate_concept_tables(donor.tables) == () + + +@pytest.mark.parametrize("reverse_persons", [False, True]) +def test_reader_repairs_two_person_cycles_that_share_one_person( + tmp_path, monkeypatch, reverse_persons +): + tables = _tables() + person = tables["person"] + person.loc[[3, 4], "A_SPOUSE"] = 0 + person.loc[[3, 4, 5], "age"] = [60, 40, 20] + # Line 1 <-> 2 <-> 3 is two reciprocal pairs without a longer cycle. + person.loc[[3, 4, 5], "PEPAR1"] = [2, 1, 2] + person.loc[[3, 4, 5], "PEPAR2"] = [-1, 3, -1] + if reverse_persons: + tables["person"] = person.iloc[::-1].reset_index(drop=True) + donor = _fake_donor(monkeypatch, tmp_path, tables) + parents = donor.tables["person"].set_index("person_id") + assert pd.isna(parents.loc[person.loc[3, "person_id"], "parent_1_person_id"]) + assert ( + parents.loc[person.loc[4, "person_id"], "parent_1_person_id"] + == person.loc[3, "person_id"] + ) + assert ( + parents.loc[person.loc[5, "person_id"], "parent_1_person_id"] + == person.loc[4, "person_id"] + ) + assert ( + parents.loc[person.loc[[3, 4, 5], "person_id"], "parent_2_person_id"] + .isna() + .all() + ) + assert donor.dropped_parent_cycle_edges == 2 + assert validate_concept_tables(donor.tables) == () + + +@pytest.mark.parametrize("also_two_person_cycle", [False, True]) +def test_reader_still_refuses_a_three_person_parent_cycle( + tmp_path, monkeypatch, also_two_person_cycle +): + tables = _tables() + person = tables["person"] + person.loc[[3, 4], "A_SPOUSE"] = 0 + person.loc[[3, 4, 5], "PEPAR1"] = [2, 3, 1] + person.loc[[3, 4, 5], "PEPAR2"] = -1 + if also_two_person_cycle: + # Line 1 -> 2 -> 3 -> 1 must still refuse even though dropping the + # younger-parent edge of line 1 <-> 2 would erase the longer cycle. + person.loc[4, "PEPAR2"] = 1 + with pytest.raises(ValueError, match="Parent pointers form a cycle"): + _fake_donor(monkeypatch, tmp_path, tables) + + +@pytest.mark.parametrize( + ("inconsistency", "message"), + [ + ("self_parent", "points at the person"), + ("duplicate_parents", "The two parents must differ"), + ("partner_parent", "partner cannot also be their parent"), + ], +) +def test_reader_cycle_rule_preserves_other_parent_refusals( + tmp_path, monkeypatch, inconsistency, message +): + tables = _tables() + person = tables["person"] + if inconsistency == "self_parent": + # A self edge remains an inconsistency even alongside a repairable + # reciprocal cycle. Repair must not silently delete the self edge. + person.loc[[3, 4], "A_SPOUSE"] = 0 + person.loc[4, ["age", "PEPAR1", "PEPAR2"]] = [51, 3, 2] + person.loc[5, ["age", "PEPAR1", "PEPAR2"]] = [31, 2, -1] + elif inconsistency == "duplicate_parents": + # The age repair would remove both copies of line 2 -> 3 and hide + # the duplicate-parent defect unless it is checked before repair. + person.loc[[3, 4], "A_SPOUSE"] = 0 + person.loc[4, ["age", "PEPAR1", "PEPAR2"]] = [51, 3, 3] + person.loc[5, ["age", "PEPAR1", "PEPAR2"]] = [31, 2, -1] + else: + person.loc[3, "PEPAR1"] = 2 + person.loc[4, "PEPAR1"] = 1 + with pytest.raises(ValueError, match=message): + _fake_donor(monkeypatch, tmp_path, tables) + + @READER_PROPERTY @given( offset=st.integers(min_value=2**53 + 1, max_value=2**63 - 7), @@ -526,24 +691,26 @@ def test_cached_real_donor_size_only_when_present(): assert path.stat().st_size == pin["size"] -def test_pinned_real_donor_reaches_the_reference_person_ruling_when_cached(): - # About 10 s: hashes the 463 MB pin, reads both tables and stages every - # line pointer, which must all resolve, before the reference-person - # check. #1130 is waiting on a ruling: in 113 pinned households (all - # 2022/2023 rows, whose A_EXPRRP the build derived from line 1) the - # A_EXPRRP reference person is not the is_household_head (P_SEQ == 1) - # person, so the reader refuses the file. Once that is ruled, assert the - # bank instead; one source household's two copies also hold a - # parent-pointer cycle that concept validation refuses. +def test_pinned_real_donor_reads_the_complete_bank_when_cached(): + # Hash and read the 463 MB pin. All households use is_household_head; + # the ASEC and PUF copies of source household 9424 each drop one + # reciprocal parent edge that names a parent younger than the child. pin = _fixture()["real_donor_pin"] path = Path.home() / pin["cache_relative_path"] if not path.is_file(): pytest.skip("Pinned real donor is not locally cached") - with pytest.raises( - ValueError, - match=r"^Donor A_EXPRRP reference person and is_household_head disagree\.$", - ): - read_populace_us_donor(path, sha256=pin["sha256"], size=pin["size"]) + donor = read_populace_us_donor(path, sha256=pin["sha256"], size=pin["size"]) + household = donor.tables["household"] + person = donor.tables["person"] + assert len(household) == 57_240 + assert len(person) == 166_321 + assert household["reference_person_id"].notna().sum() == 57_240 + assert household["reference_person_id"].is_unique + assert person["partner_person_id"].notna().sum() == 74_572 + assert person["parent_1_person_id"].notna().sum() == 62_096 + assert person["parent_2_person_id"].notna().sum() == 43_731 + assert donor.dropped_parent_cycle_edges == 2 + assert validate_concept_tables(donor.tables) == () @PROPERTY From ce4fd0e2ad96bd24239751d41002e47330ca9b8d Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Thu, 8 Oct 2026 11:29:50 -0400 Subject: [PATCH 4/6] Evaluate the truncated Pareto inverse without cancellation; refuse inexact lines (#1130) The delta review at 9df282d3 found that the top band's truncated-Pareto inverse lost precision for small alpha: it formed 1 - (lower/upper)**alpha, which cancels, so alpha=1e-18 on [10, 100] mapped the median to 10 instead of 31.6227766 and alpha=1e-16 was 4% off. _pareto_inverse now writes the truncated mass as -expm1(-alpha * w) with w = log(upper) - log(lower) and the offset as -log1p(-fraction * mass) / alpha, uses the exact alpha -> 0 limit (log-uniform, lower * exp(fraction * w)) when alpha * w underflows, clamps the offset to w, and forms lower * exp(offset) through log space only where the plain product would overflow inside a finite band. The untruncated path is the same expression with mass 1. Donor line numbers and pointers are validated as float64 (pointer columns can hold NaN), which holds every integer exactly only below 2**53; 2**53 + 1 rounded onto line 2**53. _line_numbers now refuses any magnitude >= 2**53 in the column's own dtype before conversion (infinities still fall to the integer check). CPS rosters are far below that. Tests (written by the previous GPT-6.1 Sol lane attempt, which stopped at the disk floor before any fix; the fix and the line-number unit test are the hub's): small-alpha medians against a 100-digit Decimal reference, a Hypothesis property against that reference over alpha in [1e-18, 50] and bands up to 400 decades (rtol 1e-12), the alpha -> 0 log-uniform limit at 1e-18, 1e-24 and 1e-100, and _line_numbers refusing 2**53, 2**53 + 1 and -(2**53) - 1 while keeping 2**53 - 1 exact. On 9df282d3's transport.py the six Pareto cases and both line-number cases fail; with the fix all pass. The lane's HDF-level dangling-pointer test passes on both versions (the reader already refuses that fixture for another reason) and is kept as coverage. test_transport_quantile.py 66 passed; test_transport_donor.py 90 passed, including the full read of the pinned real donor (no skips). ruff check and format clean; tools/ci_test_plan.py verify ok. Co-Authored-By: GPT-6.1 Sol Co-Authored-By: Claude Opus 5.5 --- .../src/microcosm/frame/transport.py | 59 +++++++++++-- .../shared/test_transport_donor.py | 37 +++++++++ .../shared/test_transport_quantile.py | 82 ++++++++++++++++++- 3 files changed, 168 insertions(+), 10 deletions(-) diff --git a/packages/microcosm-frame/src/microcosm/frame/transport.py b/packages/microcosm-frame/src/microcosm/frame/transport.py index 66a732bbf6..964801f465 100644 --- a/packages/microcosm-frame/src/microcosm/frame/transport.py +++ b/packages/microcosm-frame/src/microcosm/frame/transport.py @@ -104,6 +104,44 @@ def _quantile_label_codes(labels: Sequence[object], size: int, name: str) -> np. return codes +def _pareto_inverse( + fraction: np.ndarray, lower: float, upper: float | None, alpha: float +) -> np.ndarray: + """Inverse CDF of a Pareto(alpha) band starting at ``lower``. + + A finite ``upper`` truncates the distribution to ``[lower, upper]``, whose + CDF is ``(1 - (lower/x)**alpha) / (1 - (lower/upper)**alpha)``. Writing + ``w = log(upper) - log(lower)``, the truncated mass is + ``1 - (lower/upper)**alpha = -expm1(-alpha * w)`` and the inverse is + ``lower * exp(-log1p(-fraction * mass) / alpha)``. ``expm1`` and ``log1p`` + keep both steps accurate when ``alpha * w`` is tiny, where subtracting + ``(lower/upper)**alpha`` from 1 cancels: at ``alpha=1e-18`` on ``[10, 100]`` + the naive form returns 10 for the median instead of 31.62. As ``alpha`` + tends to 0 the truncated Pareto tends to the log-uniform distribution, + ``lower * exp(fraction * w)``, which is used when ``alpha * w`` underflows. + The result is ``lower * exp(offset)``; where ``exp(offset)`` alone would + overflow inside a finite band, it is formed as ``exp(log(lower) + offset)``. + """ + + with np.errstate(over="ignore", invalid="ignore", divide="ignore"): + if upper is None: + offset = -np.log1p(-fraction) / alpha + else: + width = np.log(upper) - np.log(lower) + mass = -np.expm1(-alpha * width) + if mass == 0: + offset = fraction * width + else: + # offset <= width exactly; clamp the last-bit rounding above it. + offset = np.minimum(-np.log1p(-fraction * mass) / alpha, width) + result = lower * np.exp(offset) + if upper is not None: + result = np.where( + np.isfinite(result), result, np.exp(np.log(lower) + offset) + ) + return result + + def _quantile_component( magnitudes: np.ndarray, weights: np.ndarray, @@ -131,12 +169,7 @@ def _quantile_component( fraction = (ranks[selected] - start) / (boundaries[position] - start) if method == "pareto" and position == len(bands) - 1: assert alpha is not None - # A finite upper bound declares a truncated Pareto distribution. - tail = 0.0 if band.upper is None else (band.lower / band.upper) ** alpha - with np.errstate(over="ignore", invalid="ignore", divide="ignore"): - mapped[selected] = band.lower * np.exp( - -np.log1p(-fraction * (1 - tail)) / alpha - ) + mapped[selected] = _pareto_inverse(fraction, band.lower, band.upper, alpha) else: assert band.upper is not None mapped[selected] = band.lower + fraction * (band.upper - band.lower) @@ -419,10 +452,18 @@ class DonorBank: dropped_parent_cycle_edges: int = 0 +#: Line numbers and pointers are validated as float64 (pointer columns can hold +#: NaN), which represents every integer exactly only below 2**53. A larger value +#: could round onto a different roster line, so it is refused before conversion. +_MAX_EXACT_LINE = 2**53 + + def _line_numbers(values: pd.Series, *, pointer: bool = False) -> np.ndarray: - numbers = pd.to_numeric(values, errors="raise").to_numpy( - dtype=np.float64, na_value=np.nan - ) + native = pd.to_numeric(values, errors="raise") + numbers = native.to_numpy(dtype=np.float64, na_value=np.nan) + # Compare in the column's own dtype; infinities fall to the integer check. + if ((native.abs() >= _MAX_EXACT_LINE) & ~np.isinf(numbers)).any(): + raise ValueError(f"Donor {values.name} contains invalid line numbers.") if pointer: # A missing pointer is absent (rows from a vintage without PECOHAB # leave it NaN). diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py index 3f89358c66..101c7e2a6e 100644 --- a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py @@ -19,6 +19,7 @@ from microcosm.frame.concept_mapping import ConceptMapping from microcosm.frame.concepts import ContentBasis, validate_concept_tables from microcosm.frame.transport import ( + _line_numbers, currency_bridge, derive_transport_seed, quantile_map, @@ -285,6 +286,42 @@ def test_reader_refuses_invalid_rosters(tmp_path, monkeypatch, change, message): _fake_donor(monkeypatch, tmp_path, tables) +@pytest.mark.parametrize("pointer", [False, True]) +def test_line_numbers_float64_cannot_hold_exactly_are_refused(pointer): + # 2**53 + 1 rounds to 2**53 in float64, so it would point at line 2**53. + name = "PEPAR1" if pointer else "A_LINENO" + for value in (2**53, 2**53 + 1, -(2**53) - 1): + with pytest.raises(ValueError, match="invalid line numbers"): + _line_numbers(pd.Series([1, value], name=name), pointer=pointer) + largest = _line_numbers(pd.Series([1, 2**53 - 1], name=name), pointer=pointer) + assert largest.tolist() == [1, 2**53 - 1] + + +@pytest.mark.parametrize("column", ["A_SPOUSE", "PECOHAB", "PEPAR1", "PEPAR2"]) +def test_reader_refuses_dangling_pointers_that_float64_would_round_to_a_roster_line( + tmp_path, column +): + tables = _tables() + person = tables["person"] + existing_line = 2**53 + person.loc[0, "A_LINENO"] = existing_line + for pointer_column in ("A_SPOUSE", "PECOHAB", "PEPAR1", "PEPAR2"): + # Keep the HDF input's integer identities exact before reader validation. + pointers = person[pointer_column].fillna(0).astype(np.int64) + person[pointer_column] = pointers.mask(pointers == 1, existing_line) + row = 1 if column in ("A_SPOUSE", "PECOHAB") else 2 + if column == "PECOHAB": + person.loc[row, "A_SPOUSE"] = 0 + elif column == "PEPAR2": + person.loc[row, "PEPAR1"] = -1 + person.loc[row, column] = existing_line + 1 + assert person.loc[row, column] != person.loc[0, "A_LINENO"] + path = tmp_path / "oversized-roster.h5" + pin = _write_h5(path, tables) + with pytest.raises(ValueError, match="unknown household line|invalid line numbers"): + read_populace_us_donor(path, **pin) + + @pytest.mark.parametrize( "flags", [ diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py index de92a71478..98453daa70 100644 --- a/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py @@ -1,6 +1,8 @@ """Declared weighted transport distributions, without country engines.""" -from math import fsum +from decimal import Decimal, localcontext +from fractions import Fraction +from math import fsum, log10 import numpy as np import pytest @@ -71,6 +73,37 @@ def _reference_uniform(values, weights, bands, components): return np.array(result) +def _reference_truncated_pareto(lower, upper, alpha, rank): + """Evaluate the declared CDF directly, independently of float64 arithmetic.""" + with localcontext() as context: + context.prec = 100 + lower, upper, alpha = map(Decimal.from_float, (lower, upper, alpha)) + rank = Decimal(rank.numerator) / Decimal(rank.denominator) + tail = (lower / upper) ** alpha + return float(lower * (1 - rank * (1 - tail)) ** (-1 / alpha)) + + +@st.composite +def _truncated_pareto_cases(draw): + lower_exponent = draw(st.integers(-200, 200)) + width_exponent = draw(st.integers(1, min(400, 300 - lower_exponent))) + alpha = draw( + st.one_of( + st.sampled_from([1e-18, 1e-16, 1e-12, 0.1, 1.0, 50.0]), + st.floats(-18, log10(50), allow_nan=False, allow_infinity=False).map( + lambda exponent: min(50.0, 10**exponent) + ), + ) + ) + weights = draw(st.tuples(st.integers(1, 100), st.integers(1, 100))) + return ( + 10.0**lower_exponent, + 10.0 ** (lower_exponent + width_exponent), + alpha, + weights, + ) + + def test_synthetic_distribution_round_trip(): values = np.array([1.0, 2.0, 3.0, 4.0]) mapped = quantile_map(values, np.ones(4), _bands(), interpolation="uniform") @@ -156,6 +189,53 @@ def test_pareto_top_band_uses_explicit_alpha(upper): ) +@pytest.mark.parametrize("alpha", [1e-18, 1e-16]) +def test_truncated_pareto_small_alpha_midpoint_matches_decimal(alpha): + mapped = quantile_map( + [1], + [1], + [{"lower": 10, "upper": 100, "share": 1}], + interpolation={"method": "pareto", "alpha": alpha}, + ) + expected = _reference_truncated_pareto(10.0, 100.0, alpha, Fraction(1, 2)) + np.testing.assert_allclose(mapped, [expected], rtol=1e-12, atol=0) + + +@given(_truncated_pareto_cases()) +@example(case=(10.0, 100.0, 1e-18, (1, 1))) +@example(case=(1e-200, 1e200, 1e-18, (1, 3))) +@example(case=(1e100, 1e300, 50.0, (100, 1))) +@settings(max_examples=100, deadline=None) +def test_truncated_pareto_matches_high_precision_reference(case): + lower, upper, alpha, weights = case + mapped = quantile_map( + [1, 2], + weights, + [{"lower": lower, "upper": upper, "share": 1}], + interpolation={"method": "pareto", "alpha": alpha}, + ) + total = sum(weights) + ranks = ( + Fraction(weights[0], 2 * total), + Fraction(2 * weights[0] + weights[1], 2 * total), + ) + expected = [ + _reference_truncated_pareto(lower, upper, alpha, rank) for rank in ranks + ] + np.testing.assert_allclose(mapped, expected, rtol=1e-12, atol=0) + + +@pytest.mark.parametrize("alpha", [1e-18, 1e-24, 1e-100]) +def test_truncated_pareto_tends_to_log_uniform_as_alpha_tends_to_zero(alpha): + mapped = quantile_map( + [1], + [1], + [{"lower": 10, "upper": 100, "share": 1}], + interpolation={"method": "pareto", "alpha": alpha}, + ) + np.testing.assert_allclose(mapped, [10 * (100 / 10) ** 0.5], rtol=1e-12, atol=0) + + def test_pareto_only_changes_the_top_band(): bands = [ {"lower": 0, "upper": 10, "share": 0.5}, From 833999b28ca84d72f6dd83ebfa6506fa88d431ce Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Thu, 8 Oct 2026 17:59:21 -0400 Subject: [PATCH 5/6] Close the subnormal-alpha gap and the review's accuracy nits (#1130) The delta review at ce4fd0e2 (in-session, two lenses plus a verifier, against mpmath at 120 digits) found one blocking defect in the same small-alpha class: the log-uniform limit ran only when -expm1(-alpha * w) was exactly 0, so a subnormal alpha * w kept only a few bits of the mass and returned silently wrong values (alpha=1e-320 on [10, 100]: 3.8e-4 relative; a single record at alpha=5e-324: 27.18 instead of 31.62). The limit now applies once alpha * w < 2**-60, where the dropped term is below 2**-61. The same review's accuracy nits, all fixed here: - ranks near 1 with mass near 1: log1p(-fraction * mass) cancelled (1.4e-3 relative at fraction = 1 - 2**-53 on [1e-10, 1e20], alpha=0.5); above fraction * mass = 1/2 the inverse now computes log((1 - fraction) + fraction * exp(-alpha * w)), where 1 - fraction is exact; - the band width is log1p((upper - lower) / lower) for narrow bands (upper - lower is exact there) and log(upper / lower) otherwise, and truncated results are clipped to [lower, upper] (they could exceed upper by ~1e-13 relative); - the untruncated path also falls back to log space where lower * exp(offset) alone overflows (lower=1e-300, alpha=1e-3 was refused although the quantiles are finite); - _line_numbers guards only values >= 2**53: only positive values can round onto a roster line, abs() wrapped int64's minimum, and negative pointers are absent as before 9df282d3 (negative lines stay invalid); - the HDF-level dangling-pointer test now repoints only row 0's household, uses line 2**53 - 1 with a pointer of 2**53, and matches the pointer column's own message, so each of A_SPOUSE, PECOHAB, PEPAR1 and PEPAR2 is reached (all four fail on 9df282d3); - the high-precision property test's tolerance is condition-aware (1e-12 + 4 * kappa * 2**-53, kappa = d log x / d log rank), because quantile_map rounds ranks to float64; the review's counterexample (1e-200, 1e200, alpha=10**-2.1, weights (99, 1)) is now an example. Accuracy on an 8,000-point grid against an mpmath reference (alpha down to 3e-324, bands up to 600 decades, ranks to 1 - 3e-16): no non-finite or out-of-band value, and every error within the condition-aware tolerance (worst 0.17 of it). On ce4fd0e2's transport.py, 10 of the new or extended cases fail; the four HDF dangling-pointer cases fail on 9df282d3. Refusals that depend on float64 last bits (pro-rata exact totals, distinct mapped ranks) can move in either direction for a small share of synthetic mixed-sign inputs; they still fail closed. test_transport_quantile.py 77 passed; test_transport_donor.py 90 passed, including the full read of the pinned real donor. ruff check and format clean. Co-Authored-By: Claude Opus 5.5 --- .../src/microcosm/frame/transport.py | 75 ++++++++++++------ .../shared/test_transport_donor.py | 32 ++++++-- .../shared/test_transport_quantile.py | 76 ++++++++++++++++++- 3 files changed, 153 insertions(+), 30 deletions(-) diff --git a/packages/microcosm-frame/src/microcosm/frame/transport.py b/packages/microcosm-frame/src/microcosm/frame/transport.py index 964801f465..8455e20251 100644 --- a/packages/microcosm-frame/src/microcosm/frame/transport.py +++ b/packages/microcosm-frame/src/microcosm/frame/transport.py @@ -104,6 +104,24 @@ def _quantile_label_codes(labels: Sequence[object], size: int, name: str) -> np. return codes +#: Below this alpha * log(upper / lower) the truncated Pareto equals its +#: log-uniform limit to within float64 rounding: the dropped term is about +#: alpha * w / 2 relative, under 2**-61 here. +_LOG_UNIFORM_LIMIT = 2.0**-60 + + +def _log_width(lower: float, upper: float) -> float: + """``log(upper / lower)`` without cancellation or overflow.""" + + ratio = upper / lower + if np.isfinite(ratio) and ratio <= 2: + # upper - lower is exact here (Sterbenz), so log1p keeps narrow bands. + return float(np.log1p((upper - lower) / lower)) + if np.isfinite(ratio): + return float(np.log(ratio)) + return float(np.log(upper) - np.log(lower)) + + def _pareto_inverse( fraction: np.ndarray, lower: float, upper: float | None, alpha: float ) -> np.ndarray: @@ -111,34 +129,44 @@ def _pareto_inverse( A finite ``upper`` truncates the distribution to ``[lower, upper]``, whose CDF is ``(1 - (lower/x)**alpha) / (1 - (lower/upper)**alpha)``. Writing - ``w = log(upper) - log(lower)``, the truncated mass is - ``1 - (lower/upper)**alpha = -expm1(-alpha * w)`` and the inverse is - ``lower * exp(-log1p(-fraction * mass) / alpha)``. ``expm1`` and ``log1p`` - keep both steps accurate when ``alpha * w`` is tiny, where subtracting - ``(lower/upper)**alpha`` from 1 cancels: at ``alpha=1e-18`` on ``[10, 100]`` - the naive form returns 10 for the median instead of 31.62. As ``alpha`` - tends to 0 the truncated Pareto tends to the log-uniform distribution, - ``lower * exp(fraction * w)``, which is used when ``alpha * w`` underflows. - The result is ``lower * exp(offset)``; where ``exp(offset)`` alone would - overflow inside a finite band, it is formed as ``exp(log(lower) + offset)``. + ``w = log(upper / lower)``, the truncated mass is + ``m = 1 - (lower/upper)**alpha = -expm1(-alpha * w)`` and the inverse is + ``lower * exp(-log(1 - fraction * m) / alpha)``. + + - ``expm1`` keeps ``m`` accurate when ``alpha * w`` is tiny, where + subtracting ``(lower/upper)**alpha`` from 1 cancels (at ``alpha=1e-18`` + on ``[10, 100]`` the naive form returns 10 for the median, not 31.62). + - ``log(1 - fraction * m)`` is ``log1p(-fraction * m)`` while + ``fraction * m <= 1/2``; above that it is computed as + ``log((1 - fraction) + fraction * exp(-alpha * w))``, where + ``1 - fraction`` is exact, so ranks near 1 keep their accuracy. + - As ``alpha`` tends to 0 the truncated Pareto tends to the log-uniform + distribution ``lower * exp(fraction * w)``. It is used once + ``alpha * w < 2**-60`` (the dropped term is below 2**-61), because + ``expm1`` of a subnormal argument keeps only a few bits. + - The result is ``lower * exp(offset)``, formed through log space where that + product alone would overflow, and a truncated result is clipped to + ``[lower, upper]`` (the computed ``w`` can exceed the true one by an ulp). """ with np.errstate(over="ignore", invalid="ignore", divide="ignore"): if upper is None: offset = -np.log1p(-fraction) / alpha else: - width = np.log(upper) - np.log(lower) - mass = -np.expm1(-alpha * width) - if mass == 0: + width = _log_width(lower, upper) + scaled = alpha * width + if scaled < _LOG_UNIFORM_LIMIT: offset = fraction * width else: - # offset <= width exactly; clamp the last-bit rounding above it. - offset = np.minimum(-np.log1p(-fraction * mass) / alpha, width) + mass = -np.expm1(-scaled) + product = fraction * mass + near_one = -np.log((1.0 - fraction) + fraction * np.exp(-scaled)) + offset = np.where(product > 0.5, near_one, -np.log1p(-product)) / alpha + offset = np.minimum(offset, width) result = lower * np.exp(offset) - if upper is not None: - result = np.where( - np.isfinite(result), result, np.exp(np.log(lower) + offset) - ) + result = np.where(np.isfinite(result), result, np.exp(np.log(lower) + offset)) + if upper is not None: + result = np.clip(result, lower, upper) return result @@ -454,15 +482,18 @@ class DonorBank: #: Line numbers and pointers are validated as float64 (pointer columns can hold #: NaN), which represents every integer exactly only below 2**53. A larger value -#: could round onto a different roster line, so it is refused before conversion. +#: could round onto a different roster line, so a value at or above it is +#: refused before conversion. _MAX_EXACT_LINE = 2**53 def _line_numbers(values: pd.Series, *, pointer: bool = False) -> np.ndarray: native = pd.to_numeric(values, errors="raise") numbers = native.to_numpy(dtype=np.float64, na_value=np.nan) - # Compare in the column's own dtype; infinities fall to the integer check. - if ((native.abs() >= _MAX_EXACT_LINE) & ~np.isinf(numbers)).any(): + # Compare in the column's own dtype (no abs(), which wraps int64's minimum); + # infinities fall to the integer check. Only positive values can round onto a + # roster line: every pointer <= 0 is absent and every line must be >= 1. + if ((native >= _MAX_EXACT_LINE) & ~np.isinf(numbers)).any(): raise ValueError(f"Donor {values.name} contains invalid line numbers.") if pointer: # A missing pointer is absent (rows from a vintage without PECOHAB diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py index 101c7e2a6e..0582905818 100644 --- a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py @@ -290,11 +290,24 @@ def test_reader_refuses_invalid_rosters(tmp_path, monkeypatch, change, message): def test_line_numbers_float64_cannot_hold_exactly_are_refused(pointer): # 2**53 + 1 rounds to 2**53 in float64, so it would point at line 2**53. name = "PEPAR1" if pointer else "A_LINENO" - for value in (2**53, 2**53 + 1, -(2**53) - 1): + for value in (2**53, 2**53 + 1, np.iinfo(np.int64).max): with pytest.raises(ValueError, match="invalid line numbers"): _line_numbers(pd.Series([1, value], name=name), pointer=pointer) + for value in (np.uint64(2**63),): + with pytest.raises(ValueError, match="invalid line numbers"): + _line_numbers( + pd.Series([1, value], dtype=np.uint64, name=name), pointer=pointer + ) largest = _line_numbers(pd.Series([1, 2**53 - 1], name=name), pointer=pointer) assert largest.tolist() == [1, 2**53 - 1] + # Negative values cannot round onto a roster line: a pointer <= 0 is absent + # (int64's minimum included), and a line below 1 is invalid as before. + negative = pd.Series([1, -(2**53) - 1, np.iinfo(np.int64).min], name=name) + if pointer: + assert _line_numbers(negative, pointer=True).tolist() == [1, 0, 0] + else: + with pytest.raises(ValueError, match="invalid line numbers"): + _line_numbers(negative) @pytest.mark.parametrize("column", ["A_SPOUSE", "PECOHAB", "PEPAR1", "PEPAR2"]) @@ -303,12 +316,19 @@ def test_reader_refuses_dangling_pointers_that_float64_would_round_to_a_roster_l ): tables = _tables() person = tables["person"] - existing_line = 2**53 + # The largest exactly representable line, and a pointer one above it: in + # float64 the pointer is still distinct, but nothing at or above 2**53 can be + # trusted, so the reader must refuse the pointer column itself. + existing_line = 2**53 - 1 person.loc[0, "A_LINENO"] = existing_line + household = person["person_household_id"] == person.loc[0, "person_household_id"] for pointer_column in ("A_SPOUSE", "PECOHAB", "PEPAR1", "PEPAR2"): - # Keep the HDF input's integer identities exact before reader validation. + # Keep the HDF input's integer identities exact before reader validation, + # and repoint only row 0's household at its new line. pointers = person[pointer_column].fillna(0).astype(np.int64) - person[pointer_column] = pointers.mask(pointers == 1, existing_line) + person[pointer_column] = pointers.mask( + household & (pointers == 1), existing_line + ) row = 1 if column in ("A_SPOUSE", "PECOHAB") else 2 if column == "PECOHAB": person.loc[row, "A_SPOUSE"] = 0 @@ -318,7 +338,9 @@ def test_reader_refuses_dangling_pointers_that_float64_would_round_to_a_roster_l assert person.loc[row, column] != person.loc[0, "A_LINENO"] path = tmp_path / "oversized-roster.h5" pin = _write_h5(path, tables) - with pytest.raises(ValueError, match="unknown household line|invalid line numbers"): + with pytest.raises( + ValueError, match=f"Donor {column} contains invalid line numbers" + ): read_populace_us_donor(path, **pin) diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py index 98453daa70..53acbf2ed4 100644 --- a/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py @@ -10,7 +10,7 @@ from hypothesis import strategies as st from microcosm.frame import transport -from microcosm.frame.transport import quantile_map +from microcosm.frame.transport import _pareto_inverse, quantile_map def _bands(): @@ -83,6 +83,20 @@ def _reference_truncated_pareto(lower, upper, alpha, rank): return float(lower * (1 - rank * (1 - tail)) ** (-1 / alpha)) +def _rank_condition(lower, upper, alpha, rank): + """d log(x) / d log(rank) of the truncated inverse at ``rank``. + + ``quantile_map`` rounds each rank to float64, so its output can differ from + the exact-rank reference by about this factor times 2**-53. + """ + with localcontext() as context: + context.prec = 100 + lower, upper, alpha = map(Decimal.from_float, (lower, upper, alpha)) + rank = Decimal(rank.numerator) / Decimal(rank.denominator) + mass = 1 - (lower / upper) ** alpha + return float(rank * mass / (alpha * (1 - rank * mass))) + + @st.composite def _truncated_pareto_cases(draw): lower_exponent = draw(st.integers(-200, 200)) @@ -205,6 +219,7 @@ def test_truncated_pareto_small_alpha_midpoint_matches_decimal(alpha): @example(case=(10.0, 100.0, 1e-18, (1, 1))) @example(case=(1e-200, 1e200, 1e-18, (1, 3))) @example(case=(1e100, 1e300, 50.0, (100, 1))) +@example(case=(1e-200, 1e200, 10**-2.1, (99, 1))) @settings(max_examples=100, deadline=None) def test_truncated_pareto_matches_high_precision_reference(case): lower, upper, alpha, weights = case @@ -222,10 +237,14 @@ def test_truncated_pareto_matches_high_precision_reference(case): expected = [ _reference_truncated_pareto(lower, upper, alpha, rank) for rank in ranks ] - np.testing.assert_allclose(mapped, expected, rtol=1e-12, atol=0) + for value, reference, rank in zip(mapped, expected, ranks, strict=True): + tolerance = 1e-12 + 4 * _rank_condition(lower, upper, alpha, rank) * 2.0**-53 + np.testing.assert_allclose(value, reference, rtol=tolerance, atol=0) -@pytest.mark.parametrize("alpha", [1e-18, 1e-24, 1e-100]) +@pytest.mark.parametrize( + "alpha", [1e-18, 1e-24, 1e-100, 1e-310, 1e-315, 1e-320, 1e-322, 5e-324] +) def test_truncated_pareto_tends_to_log_uniform_as_alpha_tends_to_zero(alpha): mapped = quantile_map( [1], @@ -236,6 +255,57 @@ def test_truncated_pareto_tends_to_log_uniform_as_alpha_tends_to_zero(alpha): np.testing.assert_allclose(mapped, [10 * (100 / 10) ** 0.5], rtol=1e-12, atol=0) +@pytest.mark.parametrize("alpha", [1e-310, 1e-320, 5e-324]) +def test_subnormal_alpha_quartiles_follow_the_log_uniform_limit(alpha): + # Below alpha * w = 2**-60 the inverse is the closed-form limit + # lower * (upper / lower)**rank; expm1 of a subnormal argument is not. + mapped = quantile_map( + [1, 2, 3, 4], + [1, 1, 1, 1], + [{"lower": 10, "upper": 100, "share": 1}], + interpolation={"method": "pareto", "alpha": alpha}, + ) + expected = [10 * 10 ** (rank / 8) for rank in (1, 3, 5, 7)] + np.testing.assert_allclose(mapped, expected, rtol=1e-12, atol=0) + + +def test_ranks_near_one_keep_their_accuracy(): + # 1 - fraction * mass cancels as fraction -> 1 with mass near 1; the inverse + # computes log((1 - fraction) + fraction * exp(-alpha * w)) there instead. + fraction = 1 - 2.0**-53 + lower, upper, alpha = 1e-10, 1e20, 0.5 + mapped = _pareto_inverse(np.array([fraction]), lower, upper, alpha) + with localcontext() as context: + context.prec = 100 + lo, hi, a, f = map(Decimal.from_float, (lower, upper, alpha, fraction)) + expected = float(lo * (1 - f * (1 - (lo / hi) ** a)) ** (-1 / a)) + np.testing.assert_allclose(mapped, [expected], rtol=1e-13, atol=0) + + +def test_truncated_values_never_exceed_the_band(): + mapped = quantile_map( + [1, 2], + [1e16, 1], + [{"lower": 1e-200, "upper": 1e100, "share": 1}], + interpolation={"method": "pareto", "alpha": 1e-30}, + ) + assert (mapped <= 1e100).all() and (mapped >= 1e-200).all() + + +def test_untruncated_band_reaches_finite_values_through_log_space(): + # lower * exp(offset) overflows here although the true quantile is finite. + mapped = quantile_map( + [1, 2], + [1, 1], + [{"lower": 1e-300, "upper": None, "share": 1}], + interpolation={"method": "pareto", "alpha": 1e-3}, + ) + expected = [ + np.exp(np.log(1e-300) - 1000 * np.log1p(-rank)) for rank in (0.25, 0.75) + ] + np.testing.assert_allclose(mapped, expected, rtol=1e-12, atol=0) + + def test_pareto_only_changes_the_top_band(): bands = [ {"lower": 0, "upper": 10, "share": 0.5}, From 2089808f97bb095ba53a84eb9b5f65c608626f77 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Thu, 8 Oct 2026 18:16:30 -0400 Subject: [PATCH 6/6] Pin the narrow-band width and the clip with their own tests (#1130) The delta review at 833999b2 (APPROVE_WITH_NITS) found that the narrow-band width fix and the final clip were each untested on their own. A narrow-band test ([1e300, 1e300 * (1 + 1e-14)], alpha 2, four records) now requires quantile_map to accept the band and match the high-precision reference; with the naive log(upper) - log(lower) width it refuses the distinct ranks. A clip test requires _pareto_inverse(1 - 2**-53, 1e-300, 1e-299, 0.1) <= 1e-299; without the clip it returns 1.0000000000000003e-299. The HDF pointer test is renamed to what it checks (pointer columns at or above 2**53, refused by their own name) and its comment says 9df282d3 refused the same input with a different message. test_transport_quantile.py 79 passed; test_transport_donor.py 90 passed. Co-Authored-By: Claude Opus 5.5 --- .../shared/test_transport_donor.py | 12 +++++----- .../shared/test_transport_quantile.py | 23 +++++++++++++++++++ 2 files changed, 29 insertions(+), 6 deletions(-) diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py index 0582905818..475a4b692d 100644 --- a/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py @@ -311,14 +311,14 @@ def test_line_numbers_float64_cannot_hold_exactly_are_refused(pointer): @pytest.mark.parametrize("column", ["A_SPOUSE", "PECOHAB", "PEPAR1", "PEPAR2"]) -def test_reader_refuses_dangling_pointers_that_float64_would_round_to_a_roster_line( - tmp_path, column -): +def test_reader_refuses_pointer_columns_at_or_above_2_53(tmp_path, column): tables = _tables() person = tables["person"] - # The largest exactly representable line, and a pointer one above it: in - # float64 the pointer is still distinct, but nothing at or above 2**53 can be - # trusted, so the reader must refuse the pointer column itself. + # The largest exactly representable line, and a pointer one above it. A real + # rounding collision needs a line at or above 2**53 too, which the A_LINENO + # guard refuses first, so this checks the pointer guard as a second line of + # defence: each pointer column is refused by its own name. (9df282d3 also + # refused this input, but as an unknown household line.) existing_line = 2**53 - 1 person.loc[0, "A_LINENO"] = existing_line household = person["person_household_id"] == person.loc[0, "person_household_id"] diff --git a/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py index 53acbf2ed4..a6eeddeba3 100644 --- a/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py @@ -292,6 +292,29 @@ def test_truncated_values_never_exceed_the_band(): assert (mapped <= 1e100).all() and (mapped >= 1e-200).all() +def test_narrow_band_maps_without_refusal_and_stays_accurate(): + # log(upper) - log(lower) loses the width of a band this narrow; log1p of + # the exact difference keeps it, so distinct ranks stay distinct. + lower, upper, alpha = 1e300, 1e300 * (1 + 1e-14), 2.0 + mapped = quantile_map( + [1, 2, 3, 4], + [1, 1, 1, 1], + [{"lower": lower, "upper": upper, "share": 1}], + interpolation={"method": "pareto", "alpha": alpha}, + ) + assert ((mapped >= lower) & (mapped <= upper)).all() + expected = [ + _reference_truncated_pareto(lower, upper, alpha, Fraction(rank, 8)) + for rank in (1, 3, 5, 7) + ] + np.testing.assert_allclose(mapped, expected, rtol=1e-12, atol=0) + + +def test_rounding_above_the_band_is_clipped(): + mapped = _pareto_inverse(np.array([1 - 2.0**-53]), 1e-300, 1e-299, 0.1) + assert mapped[0] <= 1e-299 + + def test_untruncated_band_reaches_finite_values_through_log_space(): # lower * exp(offset) overflows here although the true quantile is finite. mapped = quantile_map(