diff --git a/changelog.d/nz-transport-operators.added.md b/changelog.d/nz-transport-operators.added.md new file mode 100644 index 0000000000..71765ccb2d --- /dev/null +++ b/changelog.d/nz-transport-operators.added.md @@ -0,0 +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. 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 3a34d1b265..cd5ebb9436 100644 --- a/docs/concept-schema-transport-adr.md +++ b/docs/concept-schema-transport-adr.md @@ -1,9 +1,17 @@ # 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. `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 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 @@ -284,8 +292,26 @@ 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 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. 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 new file mode 100644 index 0000000000..8455e20251 --- /dev/null +++ b/packages/microcosm-frame/src/microcosm/frame/transport.py @@ -0,0 +1,906 @@ +"""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 + + +#: 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: + """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 / 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 = _log_width(lower, upper) + scaled = alpha * width + if scaled < _LOG_UNIFORM_LIMIT: + offset = fraction * width + else: + 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) + 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 + + +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 + 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) + 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 _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, + 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. 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) + 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] + 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 + 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." + ) + 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. + ``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] + weights: np.ndarray + support_strata: np.ndarray + source_person_ids: np.ndarray + content_basis: ContentBasis = ContentBasis.TRANSPORT + donor_country: str = "us" + currency: str = "USD" + 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 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 (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 + # 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 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 = ( + "A_LINENO", + "A_SPOUSE", + "PEPAR1", + "PEPAR2", + "is_household_head", + "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) + + 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(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 + # 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.") + 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 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( + 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 _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. + + 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 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 + 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"]) + 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) + 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}") + return DonorBank( + MappingProxyType(tables), + staged["weights"], + staged["support_strata"], + staged["source_person_ids"], + dropped_parent_cycle_edges=dropped, + ) + + +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..475a4b692d --- /dev/null +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_donor.py @@ -0,0 +1,910 @@ +"""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 ( + _line_numbers, + 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 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"] + + +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"), + ("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_heads", "multiple is_household_head persons"), + ("missing_head_with_reference", "no is_household_head"), + ("missing_head_without_reference", "no is_household_head"), + ("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)], + "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_heads": [(4, "is_household_head", True)], + "missing_head_with_reference": [(0, "is_household_head", False)], + "missing_head_without_reference": [ + (0, "A_EXPRRP", 5), + (0, "is_household_head", False), + ], + "channel_mismatch": [(0, "person_support_channel", "another-channel")], + } + 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("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, 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"]) +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. 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"] + for pointer_column in ("A_SPOUSE", "PECOHAB", "PEPAR1", "PEPAR2"): + # 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( + household & (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=f"Donor {column} contains invalid line numbers" + ): + read_populace_us_donor(path, **pin) + + +@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), + ) + + +@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'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"], + person.loc[4, "person_id"], + ] + 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), + 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 > 0 + ] + 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) + + +@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"] + # 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) + 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) == () + + +def test_reader_recodes_cohabiting_partner_pointers_symmetrically( + tmp_path, monkeypatch +): + tables = _tables() + person = tables["person"] + # 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: + 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(): + 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"] + + +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") + 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 +@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..a6eeddeba3 --- /dev/null +++ b/packages/microcosm-frame/tests/engine_free/shared/test_transport_quantile.py @@ -0,0 +1,800 @@ +"""Declared weighted transport distributions, without country engines.""" + +from decimal import Decimal, localcontext +from fractions import Fraction +from math import fsum, log10 + +import numpy as np +import pytest +from hypothesis import example, given, settings +from hypothesis import strategies as st + +from microcosm.frame import transport +from microcosm.frame.transport import _pareto_inverse, 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}, + ] + + +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, + 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 _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)) + + +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)) + 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") + 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, + ) + + +@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))) +@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 + 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 + ] + 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, 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], + [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) + + +@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_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( + [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}, + {"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 + ) +) +# 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) + 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_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"): + 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..8142107314 --- /dev/null +++ b/packages/microcosm-frame/tests/fixtures/transport/synthetic_donor.json @@ -0,0 +1,248 @@ +{ + "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": -1, + "PEPAR2": -1, + "PECOHAB": -1, + "A_EXPRRP": 1, + "is_household_head": true, + "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": -1, + "PEPAR2": -1, + "PECOHAB": -1, + "A_EXPRRP": 4, + "is_household_head": false, + "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, + "PECOHAB": -1, + "A_EXPRRP": 5, + "is_household_head": false, + "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": -1, + "PEPAR2": -1, + "PECOHAB": null, + "A_EXPRRP": 1, + "is_household_head": true, + "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": -1, + "PEPAR2": -1, + "PECOHAB": null, + "A_EXPRRP": 4, + "is_household_head": false, + "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": -1, + "PEPAR2": 2, + "PECOHAB": null, + "A_EXPRRP": 5, + "is_household_head": false, + "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 + } + ] +}