From 61e8857e41ee724ad4aefc642ababd219b10519f Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 10:45:41 -0400 Subject: [PATCH 01/17] SM1: differential mortality normalized within age-sex-year cells Adds the opt-in package populace_dynamics.mortality_ses. - Relative risks are specified over three dimensions: lifetime-earnings group (primary), education, and race and Hispanic origin. They are log-additive and indexed by sex, age and year. Years past the table hold the last year's coefficients. A coefficient path is pluggable. - The unclassified rule is `refuse` by default; `reference` contributes zero log risk. - A pure normalizer rescales probabilities so that each single-age x sex x year cell keeps the baseline's expected deaths. Fixed rows keep their fixed probability, and the others get min(1, q_base x RR x k), with iterative capping. It logs the residual and flags infeasible cells. - DifferentialMortality wraps a year-aware baseline. Its fixed rows are DI-origin persons at baseline q, so the DI-aware adapter sees the baseline rate for them. It exposes bands and __call__(frame, context), and deliberately no probabilities attribute. - apply_differential_mortality draws apply_mortality's uniforms. Nothing is wired into a registered run, runner, config or gates.yaml. Tests are on invented data: Hypothesis properties for cell conservation, feasibility, bit-identity at unit risk, bounds, DI composition, intended monotonicity, permutation and weight invariance, and determinism; plus a differential test against the DI adapter's _net_of_di_origin arithmetic. Built by a Subfleet build lane (GPT-6.1 Sol) from the lane brief; that lane's sandbox could not write .git, so the lane session commits it. Co-Authored-By: Claude Opus 5.5 --- .../mortality_ses/__init__.py | 24 + src/populace_dynamics/mortality_ses/model.py | 217 +++++ .../mortality_ses/normalizer.py | 342 ++++++++ .../mortality_ses/specification.py | 254 ++++++ src/populace_dynamics/mortality_ses/step.py | 63 ++ tests/mortality_ses/test_model_step.py | 489 +++++++++++ tests/mortality_ses/test_normalizer.py | 264 ++++++ tests/mortality_ses/test_properties.py | 762 ++++++++++++++++++ tests/mortality_ses/test_specification.py | 223 +++++ 9 files changed, 2638 insertions(+) create mode 100644 src/populace_dynamics/mortality_ses/__init__.py create mode 100644 src/populace_dynamics/mortality_ses/model.py create mode 100644 src/populace_dynamics/mortality_ses/normalizer.py create mode 100644 src/populace_dynamics/mortality_ses/specification.py create mode 100644 src/populace_dynamics/mortality_ses/step.py create mode 100644 tests/mortality_ses/test_model_step.py create mode 100644 tests/mortality_ses/test_normalizer.py create mode 100644 tests/mortality_ses/test_properties.py create mode 100644 tests/mortality_ses/test_specification.py diff --git a/src/populace_dynamics/mortality_ses/__init__.py b/src/populace_dynamics/mortality_ses/__init__.py new file mode 100644 index 00000000..ef0718c2 --- /dev/null +++ b/src/populace_dynamics/mortality_ses/__init__.py @@ -0,0 +1,24 @@ +"""Opt-in differential mortality; no fitted gradients or projection wiring.""" + +from .model import DifferentialMortality +from .normalizer import CELL_LOG_COLUMNS, normalize_probabilities +from .specification import ( + CoefficientPath, + HoldLastCoefficientPath, + RelativeRiskSpec, + RiskDimension, + UnclassifiedValueError, +) +from .step import apply_differential_mortality + +__all__ = [ + "CELL_LOG_COLUMNS", + "CoefficientPath", + "DifferentialMortality", + "HoldLastCoefficientPath", + "RelativeRiskSpec", + "RiskDimension", + "UnclassifiedValueError", + "apply_differential_mortality", + "normalize_probabilities", +] diff --git a/src/populace_dynamics/mortality_ses/model.py b/src/populace_dynamics/mortality_ses/model.py new file mode 100644 index 00000000..ee41b389 --- /dev/null +++ b/src/populace_dynamics/mortality_ses/model.py @@ -0,0 +1,217 @@ +"""Pure, year-aware differential mortality around a single-age baseline.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import ClassVar, Protocol + +import numpy as np +import pandas as pd + +from populace_dynamics.engine.di_entitlement_rates import ( + MAX_AGE, + SINGLE_YEAR_AGE_BANDS, + DIEntitlementRates, + DIEntitlementSpec, +) +from populace_dynamics.engine.loop import PeriodContext + +from .normalizer import normalize_probabilities +from .specification import RelativeRiskSpec, UnclassifiedValueError + + +class YearAwareMortality(Protocol): + """The context-aware baseline protocol; no yearless probability API.""" + + bands: tuple[tuple[int, int], ...] + + def __call__( + self, frame: pd.DataFrame, context: PeriodContext + ) -> np.ndarray: ... + + +def _integer_values(series: pd.Series, column: str) -> np.ndarray: + """Refuse missing, fractional and non-finite integer inputs.""" + try: + values = series.to_numpy(dtype=np.float64, na_value=np.nan) + except (TypeError, ValueError) as error: + raise ValueError(f"column {column!r} must contain integers") from error + if ( + not np.isfinite(values).all() + or (values != np.floor(values)).any() + or (values < np.iinfo(np.int64).min).any() + or (values >= float(np.iinfo(np.int64).max)).any() + ): + raise ValueError(f"column {column!r} must contain finite integers") + return values.astype(np.int64) + + +def _weights(frame: pd.DataFrame, column: str | None) -> np.ndarray: + if column is None: + return np.ones(len(frame), dtype=np.float64) + if column not in frame: + raise ValueError( + f"mortality frame is missing weight column {column!r}" + ) + try: + values = frame[column].to_numpy(dtype=np.float64, na_value=np.nan) + except (TypeError, ValueError) as error: + raise ValueError( + "mortality weights must be finite and >= 0" + ) from error + if not np.isfinite(values).all() or (values < 0).any(): + raise ValueError("mortality weights must be finite and >= 0") + return values + + +@dataclass(frozen=True) +class DifferentialMortality: + """Redistribute baseline deaths within single age, sex and year cells. + + ``weight_column`` is required explicitly; ``None`` counts persons equally. + If birth year is present, uncapped age must equal ``year - 1 - birth_year`` + regardless of the DI adapter's optional age validation. Both this model + and the adapter then cap cell ages at 120. DI-origin rows retain baseline + q so the adapter, using the same DI configuration, owns their death rates. + + This model deliberately exposes no ``probabilities`` attribute: the DI + adapter dispatches that attribute without a context or a projection year. + ``evaluate`` additionally returns the pure normalizer's cell accounting. + """ + + baseline: YearAwareMortality + specification: RelativeRiskSpec + weight_column: str | None = field(kw_only=True) + di_configuration: DIEntitlementRates | DIEntitlementSpec | None = field( + default=None, kw_only=True + ) + bands: ClassVar[tuple[tuple[int, int], ...]] = SINGLE_YEAR_AGE_BANDS + + def __post_init__(self) -> None: + if hasattr(self, "probabilities"): + raise ValueError( + "DifferentialMortality must not expose a probabilities " + "attribute; its interface requires context" + ) + if hasattr(self.baseline, "probabilities"): + raise ValueError( + "baseline must use (frame, context), without a probabilities " + "attribute" + ) + if not callable(self.baseline): + raise TypeError("baseline must be callable as (frame, context)") + if tuple(getattr(self.baseline, "bands", ())) != self.bands: + raise ValueError("baseline bands must be SINGLE_YEAR_AGE_BANDS") + if self.weight_column is not None and not isinstance( + self.weight_column, str + ): + raise TypeError("weight_column must be a column name or None") + if self.di_configuration is not None and not isinstance( + self.di_configuration, (DIEntitlementRates, DIEntitlementSpec) + ): + raise TypeError("di_configuration must be DI rates or their spec") + + def _fixed_rows(self, frame: pd.DataFrame) -> np.ndarray: + columns = {"di_entitled", "di_conversion_year"} + present = { + column + for column in frame.columns + if isinstance(column, str) and column.startswith("di_") + } + if self.di_configuration is None: + if present: + raise ValueError( + "DI columns require di_configuration matching the adapter" + ) + return np.zeros(len(frame), dtype=bool) + missing = columns.difference(frame.columns) + if missing: + raise ValueError( + f"mortality frame is missing columns {sorted(missing)}" + ) + entitled = frame["di_entitled"] + if entitled.isna().any(): + raise ValueError("column 'di_entitled' has missing values") + if not ( + pd.api.types.is_bool_dtype(entitled.dtype) + or entitled.isin([True, False]).all() + ): + raise ValueError("column 'di_entitled' must be boolean") + conversion = frame["di_conversion_year"] + converted = conversion.notna().to_numpy(dtype=bool) + _integer_values(conversion.loc[converted], "di_conversion_year") + spec = ( + self.di_configuration.spec + if isinstance(self.di_configuration, DIEntitlementRates) + else self.di_configuration + ) + return entitled.to_numpy(dtype=bool) | ( + converted & (spec.post_conversion_mortality == "di_origin") + ) + + def evaluate( + self, frame: pd.DataFrame, context: PeriodContext + ) -> tuple[np.ndarray, pd.DataFrame]: + """Return probabilities and cell accounting without consuming RNG.""" + missing = {"age", "sex"}.difference(frame.columns) + if missing: + raise ValueError( + f"mortality frame is missing columns {sorted(missing)}" + ) + year = int(context.year) + if year != context.year: + raise ValueError("mortality context year must be an integer") + age = _integer_values(frame["age"], "age") + if (age < 0).any(): + raise ValueError("mortality ages must be non-negative") + if "birth_year" in frame: + birth_year = _integer_values(frame["birth_year"], "birth_year") + if (age != year - 1 - birth_year).any(): + raise ValueError( + "mortality age must equal context.year - 1 - birth_year " + "before capping at 120" + ) + ages = np.minimum(age, MAX_AGE) + sex = frame["sex"].astype(str).to_numpy() + if not np.isin(sex, ("female", "male")).all(): + raise ValueError("mortality sex must be 'female' or 'male'") + weights = _weights(frame, self.weight_column) + fixed = self._fixed_rows(frame) + baseline = np.asarray(self.baseline(frame, context), dtype=np.float64) + if baseline.shape != (len(frame),): + raise ValueError("baseline mortality returned the wrong shape") + if ( + not np.isfinite(baseline).all() + or ((baseline < 0) | (baseline > 1)).any() + ): + raise ValueError("baseline mortality must lie in [0, 1]") + cells = list(zip([year] * len(frame), sex, ages, strict=True)) + try: + risk, unclassified = self.specification.relative_risks( + frame, year, ages + ) + except UnclassifiedValueError as error: + _, error.cell_log = normalize_probabilities( + baseline, + np.ones(len(frame)), + fixed, + baseline, + cells, + weights, + unclassified_counts=error.unclassified_counts, + ) + raise + return normalize_probabilities( + baseline, + risk, + fixed, + baseline, + cells, + weights, + unclassified_counts=unclassified, + ) + + def __call__( + self, frame: pd.DataFrame, context: PeriodContext + ) -> np.ndarray: + return self.evaluate(frame, context)[0] diff --git a/src/populace_dynamics/mortality_ses/normalizer.py b/src/populace_dynamics/mortality_ses/normalizer.py new file mode 100644 index 00000000..e009e8ad --- /dev/null +++ b/src/populace_dynamics/mortality_ses/normalizer.py @@ -0,0 +1,342 @@ +"""Deterministic cell accounting for differential mortality. + +The coefficients used here carry no empirical claim. This module only +redistributes each cell's baseline expected deaths, subject to prescribed +fixed probabilities and the probability bounds. +""" + +from __future__ import annotations + +from collections.abc import Mapping, Sequence +from decimal import Decimal, localcontext +from math import fsum + +import numpy as np +import pandas as pd + +CELL_LOG_COLUMNS = ( + "year", + "sex", + "age", + "persons", + "fixed_persons", + "weight", + "fixed_weight", + "baseline_expected_deaths", + "applied_expected_deaths", + "factor", + "clipped_persons", + "clipping_residual", + "infeasible", + "unclassified_refuse", + "unclassified_reference", +) + + +def _stable_sum(values: np.ndarray) -> np.longdouble: + """Sum nonnegative values without a row-order-dependent reduction. + + Scaling before ``fsum`` avoids an overflowing intermediate sum. + ``fsum`` gives the same result after a row permutation. Cells whose + products need a wider exponent range use the Decimal path below. + """ + if not values.size: + return np.longdouble(0) + maximum = np.max(values) + if maximum == 0: + return np.longdouble(0) + return maximum * np.longdouble(fsum(float(x / maximum) for x in values)) + + +def _decimal_cell(base, rr, w, is_fixed, prescribed): + """Handle extreme finite float64 inputs without platform assumptions. + + On some supported platforms ``longdouble`` has float64's exponent + range. Decimal's exponent range retains positive q*RR support even + when that product would underflow, and avoids overflowing products. + This slower path is restricted to cells containing extreme magnitudes. + """ + with localcontext() as decimal_context: + # Enough to retain exact products of three finite binary64 inputs, + # including subnormals, before the normalization division. + decimal_context.prec = 2400 + zero, one = Decimal(0), Decimal(1) + base, rr, w, prescribed = ( + [Decimal.from_float(float(value)) for value in values] + for values in (base, rr, w, prescribed) + ) + + def total(values): + return sum(sorted(values), zero) + + free = [i for i, fixed in enumerate(is_fixed) if not fixed] + fixed = [i for i, fixed in enumerate(is_fixed) if fixed] + result = [ + prescribed[i] if is_fixed[i] else zero for i in range(len(w)) + ] + baseline_deaths = total( + wi * qi for wi, qi in zip(w, base, strict=True) + ) + fixed_deaths = total(w[i] * prescribed[i] for i in fixed) + fixed_at_baseline = all(prescribed[i] == base[i] for i in fixed) + target = ( + total(w[i] * base[i] for i in free) + if fixed_at_baseline + else baseline_deaths - fixed_deaths + ) + support = [i for i in free if base[i] > zero] + capacity = total(w[i] for i in support) + capped = set() + infeasible = False + factor = one + if fixed_at_baseline and ( + not free or all(rr[i] == rr[free[0]] for i in free) + ): + result = base.copy() + factor = one / rr[free[0]] if free else one + elif target < zero: + infeasible, factor = True, zero + elif target == zero: + factor = zero + elif target >= capacity: + infeasible = target > capacity + for i in support: + result[i] = one + capped = set(support) + factor = max( + (one / (base[i] * rr[i]) for i in support), default=one + ) + else: + raw = [qi * ri for qi, ri in zip(base, rr, strict=True)] + active = {i for i in support if w[i] > zero} + remaining = target + while active: + factor = remaining / total(w[i] * raw[i] for i in active) + next_capped = {i for i in active if raw[i] * factor > one} + if not next_capped: + break + capped |= next_capped + active -= next_capped + remaining = target - total(w[i] for i in capped) + for i in free: + result[i] = min(one, raw[i] * factor) + capped = {i for i in free if raw[i] * factor > one} + probabilities = np.asarray([float(value) for value in result]) + applied = total( + wi * Decimal.from_float(float(qi)) + for wi, qi in zip(w, probabilities, strict=True) + ) + residual = baseline_deaths - applied if infeasible or capped else zero + return probabilities, { + "weight": float(total(w)), + "fixed_weight": float(total(w[i] for i in fixed)), + "baseline_expected_deaths": float(baseline_deaths), + "applied_expected_deaths": float(applied), + "factor": float(factor), + "clipped_persons": len(capped), + "clipping_residual": float(residual), + "infeasible": infeasible, + } + + +def _normalize_cell(base, rr, w, is_fixed, prescribed): + # Within these magnitudes, coefficients, weighted products, and capped + # scaling products remain representable even on float64-only platforms. + if any( + np.any((values > 0) & (values < 1e-80)) or np.any(values > 1e80) + for values in (base, rr, w, prescribed) + ): + return _decimal_cell(base, rr, w, is_fixed, prescribed) + free = ~is_fixed + result = np.zeros(len(base), dtype=np.longdouble) + result[is_fixed] = prescribed[is_fixed] + baseline_deaths = _stable_sum(w * base) + fixed_deaths = _stable_sum(w[is_fixed] * prescribed[is_fixed]) + fixed_at_baseline = np.array_equal(prescribed[is_fixed], base[is_fixed]) + # Avoid cancellation of small free totals when fixed rows dominate. + target = ( + _stable_sum(w[free] * base[free]) + if fixed_at_baseline + else baseline_deaths - fixed_deaths + ) + support = free & (base > 0) + capacity = _stable_sum(w[support]) + capped = np.zeros(len(base), dtype=bool) + infeasible = False + factor = np.longdouble(1) + if fixed_at_baseline and ( + not np.any(free) or np.all(rr[free] == rr[free][0]) + ): + result = base.copy() + if np.any(free): + factor = np.longdouble(1) / rr[free][0] + elif target < 0: + infeasible = True + factor = np.longdouble(0) + elif target == 0: + factor = np.longdouble(0) + elif target >= capacity: + # Baseline fixed rows are feasible by construction. Independent + # floating reductions can differ by a rounding bit at capacity. + infeasible = bool(target > capacity) and not fixed_at_baseline + result[support] = 1 + capped[support] = True + factor = ( + np.max(1 / (base[support] * rr[support])) + if np.any(support) + else np.longdouble(1) + ) + else: + raw = base * rr + active = support & (w > 0) + remaining = target + while np.any(active): + denominator = _stable_sum(w[active] * raw[active]) + factor = remaining / denominator + next_capped = active & (raw * factor > 1) + if not np.any(next_capped): + break + capped |= next_capped + active &= ~next_capped + remaining = target - _stable_sum(w[capped]) + result[free] = np.minimum(1, raw[free] * factor) + capped[free] = raw[free] * factor > 1 + + probabilities = result.astype(float) + applied_deaths = _stable_sum(w * probabilities) + residual = baseline_deaths - applied_deaths + if not infeasible and not np.any(capped): + residual = np.longdouble(0) + return probabilities, { + "weight": _stable_sum(w), + "fixed_weight": _stable_sum(w[is_fixed]), + "baseline_expected_deaths": baseline_deaths, + "applied_expected_deaths": applied_deaths, + "factor": factor, + "clipped_persons": int(np.count_nonzero(capped)), + "clipping_residual": residual, + "infeasible": infeasible, + } + + +def _vector(value, name: str, size: int | None = None) -> np.ndarray: + array = np.asarray(value, dtype=float) + if array.ndim != 1 or (size is not None and len(array) != size): + raise ValueError(f"{name} must be a one-dimensional matching array") + if not np.all(np.isfinite(array)): + raise ValueError(f"{name} must contain finite values") + return array + + +def normalize_probabilities( + q_base: np.ndarray, + relative_risk: np.ndarray, + fixed_mask: np.ndarray, + fixed_probabilities: np.ndarray, + cell_ids: Sequence[tuple[int, str, int]], + weights: np.ndarray, + *, + unclassified_counts: Mapping[str, np.ndarray] | None = None, +) -> tuple[np.ndarray, pd.DataFrame]: + """Normalize risks to baseline weighted expected deaths in each cell. + + Free rows receive ``min(1, q_base * relative_risk * factor)``. The + factor is recomputed after each capped group until the remaining + target can be allocated without further capping. Zero-weight rows + receive the same factor as other rows but do not affect its value. + + Fixed probabilities above the entire baseline target are infeasible; + free probabilities are then zero. A target above the free rows' + positive-probability support is also infeasible; those rows saturate + at one. The signed residual is baseline minus applied deaths. On + feasible cells without capping it is explicitly zero: arithmetic + roundoff is covered by the accounting tolerance, rather than being + described as a clipping loss. + + All relative risks must be finite and strictly positive, including + those on fixed rows. Probabilities lie in [0, 1], weights are finite + and nonnegative, and each cell id is a (year, sex, age) triple. + Optional per-row unclassified counts are summed into each cell log. + No input array is mutated and no randomness is consumed. + """ + baseline = _vector(q_base, "baseline probabilities") + size = len(baseline) + risk = _vector(relative_risk, "relative risks", size) + fixed_q = _vector(fixed_probabilities, "fixed probabilities", size) + weight = _vector(weights, "weights", size) + fixed = np.asarray(fixed_mask) + if fixed.ndim != 1 or len(fixed) != size or fixed.dtype.kind != "b": + raise ValueError("fixed_mask must be a matching boolean array") + if np.any((baseline < 0) | (baseline > 1)): + raise ValueError("baseline probabilities must lie in [0, 1]") + if np.any((fixed_q < 0) | (fixed_q > 1)): + raise ValueError("fixed probabilities must lie in [0, 1]") + if np.any(risk <= 0): + raise ValueError("relative risks must be strictly positive") + if np.any(weight < 0): + raise ValueError("weights must be nonnegative") + + ids = list(cell_ids) + if len(ids) != size: + raise ValueError("cell_ids must match the input arrays") + groups: dict[tuple, list[int]] = {} + for index, cell in enumerate(ids): + if not isinstance(cell, (tuple, list, np.ndarray)) or len(cell) != 3: + raise ValueError("cell ids must be (year, sex, age) triples") + key = tuple(cell) + if any(pd.isna(value) for value in key): + raise ValueError("cell ids must not contain missing values") + try: + groups.setdefault(key, []).append(index) + except TypeError as error: + raise ValueError( + "cell ids must contain hashable values" + ) from error + + counts = { + "unclassified_refuse": np.zeros(size, dtype=np.int64), + "unclassified_reference": np.zeros(size, dtype=np.int64), + } + for name, values in (unclassified_counts or {}).items(): + if not isinstance(name, str) or ( + name in CELL_LOG_COLUMNS and name not in counts + ): + raise ValueError( + "unclassified count names must be distinct log columns" + ) + array = _vector(values, f"{name} counts", size) + if np.any(array < 0) or np.any(array != np.floor(array)): + raise ValueError( + "unclassified counts must be nonnegative integers" + ) + counts[name] = array + + out = baseline.copy() + logs = [] + for cell in sorted(groups): + positions = np.asarray(groups[cell], dtype=np.int64) + base = baseline[positions].astype(np.longdouble) + rr = risk[positions].astype(np.longdouble) + w = weight[positions].astype(np.longdouble) + is_fixed = fixed[positions] + prescribed = fixed_q[positions].astype(np.longdouble) + out[positions], accounting = _normalize_cell( + base, rr, w, is_fixed, prescribed + ) + row = dict( + zip(("year", "sex", "age"), cell, strict=True), + persons=len(positions), + fixed_persons=int(np.count_nonzero(is_fixed)), + **accounting, + ) + row.update( + { + name: int(fsum(values[positions])) + for name, values in counts.items() + } + ) + logs.append(row) + columns = list(CELL_LOG_COLUMNS) + [ + name for name in counts if name not in CELL_LOG_COLUMNS + ] + return out, pd.DataFrame(logs, columns=columns) diff --git a/src/populace_dynamics/mortality_ses/specification.py b/src/populace_dynamics/mortality_ses/specification.py new file mode 100644 index 00000000..b3d916d8 --- /dev/null +++ b/src/populace_dynamics/mortality_ses/specification.py @@ -0,0 +1,254 @@ +"""Explicit, year-indexed log-additive mortality relative risks. + +No empirical coefficients or category definitions are supplied here. SM2 owns +those inputs; the mechanism requires all three named dimensions. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Protocol + +import numpy as np +import pandas as pd + +DIMENSIONS = ("lifetime_earnings", "education", "race_hispanic_origin") +SEXES = ("female", "male") +MAX_AGE = 120 +KEYS = ["dimension", "category", "sex", "age"] +LOG_RISK = "log_relative_risk" + + +class UnclassifiedValueError(ValueError): + """A refusal carrying row counts; model evaluation adds its cell log.""" + + def __init__(self, counts: dict[str, np.ndarray]): + self.unclassified_counts = counts + self.cell_log: pd.DataFrame | None = None + super().__init__( + "unclassified mortality attributes under refuse: " + f"{int(counts['unclassified_refuse'].sum())} contributions" + ) + + +@dataclass(frozen=True) +class RiskDimension: + """Ordered categories and a zero-log-risk reference for one dimension.""" + + name: str + column: str + categories: tuple[str, ...] + reference: str + unclassified_values: tuple[str, ...] = () + + def __post_init__(self) -> None: + object.__setattr__(self, "categories", tuple(self.categories)) + object.__setattr__( + self, "unclassified_values", tuple(self.unclassified_values) + ) + if self.name not in DIMENSIONS or not self.column: + raise ValueError("dimension needs a known name and frame column") + if ( + not self.categories + or any(not isinstance(c, str) or not c for c in self.categories) + or len(set(self.categories)) != len(self.categories) + or self.reference not in self.categories + ): + raise ValueError("categories must be unique with a reference") + if set(self.categories) & set(self.unclassified_values): + raise ValueError("unclassified values overlap categories") + + +class CoefficientPath(Protocol): + """SM2 extension point: complete coefficient tables for named years. + + ``coefficients`` returns dimension/category/sex/age/log_relative_risk + rows (or inclusive age_lower/age_upper bands). Every table is validated, + including a custom path's output at evaluation time. + """ + + @property + def years(self) -> tuple[int, ...]: ... + + def coefficients(self, year: int) -> pd.DataFrame: ... + + +class HoldLastCoefficientPath: + """A step-held table; after its last year the final gradient persists.""" + + def __init__(self, table: pd.DataFrame): + if "year" not in table or table.empty: + raise ValueError("coefficient table needs nonempty year coverage") + years = table["year"].to_numpy(dtype=float) + if ( + not np.isfinite(years).all() + or (years != np.floor(years)).any() + or (years < np.iinfo(np.int64).min).any() + or (years >= float(np.iinfo(np.int64).max)).any() + ): + raise ValueError("coefficient years must be finite integers") + self._table = table.copy(deep=True) + self._table["year"] = years.astype(int) + self._years = tuple(sorted(set(years.astype(int)))) + + @property + def years(self) -> tuple[int, ...]: + return self._years + + def coefficients(self, year: int) -> pd.DataFrame: + if year < self.years[0]: + raise ValueError("year precedes first mortality coefficient year") + selected = max(y for y in self.years if y <= year) + return ( + self._table.loc[self._table["year"] == selected] + .drop(columns="year") + .copy(deep=True) + ) + + +def _expanded(table: pd.DataFrame) -> pd.DataFrame: + missing = set(KEYS[:-1] + [LOG_RISK]) - set(table.columns) + if missing: + raise ValueError( + f"coefficient table missing columns {sorted(missing)}" + ) + if "age" in table: + out = table.copy(deep=True) + elif {"age_lower", "age_upper"} <= set(table.columns): + rows = [] + for row in table.to_dict("records"): + lower, upper = row.pop("age_lower"), row.pop("age_upper") + if ( + not np.isfinite([lower, upper]).all() + or int(lower) != lower + or int(upper) != upper + or not 0 <= lower <= upper <= MAX_AGE + ): + raise ValueError("invalid inclusive coefficient age band") + rows.extend( + dict(row, age=age) for age in range(int(lower), int(upper) + 1) + ) + out = pd.DataFrame(rows, columns=[*KEYS, LOG_RISK]) + else: + raise ValueError("coefficient table needs age or inclusive age bands") + age = out["age"].to_numpy(dtype=float) + if ( + not np.isfinite(age).all() + or (age != np.floor(age)).any() + or ((age < 0) | (age > MAX_AGE)).any() + ): + raise ValueError("coefficient ages must be integers 0-120") + out["age"] = age.astype(int) + return out + + +@dataclass(frozen=True) +class RelativeRiskSpec: + """Validated joint specification, with explicit unclassified handling.""" + + dimensions: tuple[RiskDimension, ...] + path: CoefficientPath + unclassified_rule: str = "refuse" + + def __post_init__(self) -> None: + object.__setattr__(self, "dimensions", tuple(self.dimensions)) + if len(self.dimensions) != len(DIMENSIONS) or { + d.name for d in self.dimensions + } != set(DIMENSIONS): + raise ValueError("specification must cover all three dimensions") + if len({d.column for d in self.dimensions}) != len(DIMENSIONS): + raise ValueError("each dimension needs its own frame column") + if self.unclassified_rule not in ("refuse", "reference"): + raise ValueError("unclassified rule must be refuse or reference") + years = tuple(self.path.years) + if not years or tuple(sorted(set(years))) != years: + raise ValueError("path years must be nonempty, unique and ordered") + if any(not isinstance(y, (int, np.integer)) for y in years): + raise ValueError("path years must be integers") + for year in years: + self._validated_table(year) + + def _validated_table(self, year: int) -> pd.Series: + if year < self.path.years[0]: + raise ValueError("year precedes first mortality coefficient year") + table = _expanded(self.path.coefficients(year)) + if table[KEYS].isna().any().any() or table.duplicated(KEYS).any(): + raise ValueError("coefficient keys missing or duplicated") + expected = { + (dimension.name, category, sex, age) + for dimension in self.dimensions + for category in dimension.categories + for sex in SEXES + for age in range(MAX_AGE + 1) + } + if set(table[KEYS].itertuples(index=False, name=None)) != expected: + raise ValueError("coefficient coverage or categories do not match") + values = table[LOG_RISK].to_numpy(dtype=float) + with np.errstate(over="ignore", under="ignore", invalid="ignore"): + risks = np.exp(values) + if not np.isfinite(risks).all() or (risks <= 0).any(): + raise ValueError("every relative risk must be finite and positive") + for dimension in self.dimensions: + reference = (table["dimension"] == dimension.name) & ( + table["category"] == dimension.reference + ) + if (values[reference] != 0).any(): + raise ValueError("reference category must have zero log risk") + return table.set_index(KEYS)[LOG_RISK] + + def relative_risks( + self, frame: pd.DataFrame, year: int, ages: np.ndarray + ) -> tuple[np.ndarray, dict[str, np.ndarray]]: + """Return log-additive risks and explicit per-dimension row counts.""" + required = {d.column for d in self.dimensions} | {"sex"} + if required - set(frame.columns): + raise ValueError( + "mortality frame missing dimension or sex columns" + ) + sex = frame["sex"].to_numpy() + if not frame["sex"].isin(SEXES).all(): + raise ValueError("unknown mortality sex") + ages = np.asarray(ages) + if ( + ages.shape != (len(frame),) + or not np.isin(ages, np.arange(MAX_AGE + 1)).all() + ): + raise ValueError("risk ages must be single ages 0-120") + coefficients = self._validated_table(year) + total = np.zeros(len(frame)) + counts = { + "unclassified_refuse": np.zeros(len(frame), dtype=int), + "unclassified_reference": np.zeros(len(frame), dtype=int), + } + for dimension in self.dimensions: + values = frame[dimension.column] + unclassified = ( + values.isna() | values.isin(dimension.unclassified_values) + ).to_numpy() + if ( + ~(values.isin(dimension.categories).to_numpy() | unclassified) + ).any(): + raise ValueError(f"unknown category in {dimension.name}") + key = f"unclassified_{self.unclassified_rule}" + counts[key] += unclassified + for rule in ("refuse", "reference"): + counts[f"{dimension.name}_unclassified_{rule}"] = ( + unclassified.astype(int) + if rule == self.unclassified_rule + else np.zeros(len(frame), dtype=int) + ) + categories = values.where(~unclassified, dimension.reference) + index = pd.MultiIndex.from_arrays( + [[dimension.name] * len(frame), categories, sex, ages], + names=KEYS, + ) + total += coefficients.reindex(index).to_numpy(dtype=float) + if counts["unclassified_refuse"].any(): + raise UnclassifiedValueError(counts) + with np.errstate(over="ignore", under="ignore", invalid="ignore"): + risks = np.exp(total) + if not np.isfinite(risks).all() or (risks <= 0).any(): + raise ValueError( + "combined relative risks must be finite and positive" + ) + return risks, counts diff --git a/src/populace_dynamics/mortality_ses/step.py b/src/populace_dynamics/mortality_ses/step.py new file mode 100644 index 00000000..f1fc30b1 --- /dev/null +++ b/src/populace_dynamics/mortality_ses/step.py @@ -0,0 +1,63 @@ +"""Context-aware mortality frame step using the engine's mortality streams.""" + +from __future__ import annotations + +from collections.abc import MutableMapping + +import numpy as np +import pandas as pd + +from populace_dynamics.engine.loop import PeriodContext +from populace_dynamics.engine.rng import ProjectionModule + +from .model import DifferentialMortality + + +def apply_differential_mortality( + frame: pd.DataFrame, + context: PeriodContext, + rng: np.random.Generator, + *, + model: DifferentialMortality, + weight_column: str | None, + death_log: MutableMapping[int, pd.DataFrame] | None = None, + cell_log: MutableMapping[int, pd.DataFrame] | None = None, +) -> pd.DataFrame: + """Draw the engine's uniforms and return survivors in person-ID order. + + As in ``engine.steps.apply_mortality``, uniforms are drawn in stable + person-ID order before evaluating the mortality model. With a registry, + each person uses ``ProjectionModule.MORTALITY`` and the injected generator + is untouched. Logs are mappings indexed by projection year. + """ + if weight_column != model.weight_column: + raise ValueError("step weight_column must match model.weight_column") + if "person_id" not in frame: + raise ValueError("mortality frame is missing column 'person_id'") + ordered = frame.sort_values("person_id", kind="stable") + if context.rng_registry is None: + uniform = rng.random(len(ordered)) + else: + uniform = np.asarray( + [ + context.person_generator( + ProjectionModule.MORTALITY, person_id + ).random() + for person_id in ordered["person_id"] + ] + ) + probability, cells = model.evaluate(ordered, context) + death = uniform < probability + year = int(context.year) + if cell_log is not None: + cell_log[year] = cells + if death_log is not None: + death_log[year] = pd.DataFrame( + { + "person_id": ordered["person_id"].to_numpy()[death], + "sex": ordered["sex"].to_numpy()[death], + "start_age": ordered["age"].to_numpy()[death], + "q_applied": probability[death], + } + ) + return ordered.loc[~death].reset_index(drop=True) diff --git a/tests/mortality_ses/test_model_step.py b/tests/mortality_ses/test_model_step.py new file mode 100644 index 00000000..cc735d51 --- /dev/null +++ b/tests/mortality_ses/test_model_step.py @@ -0,0 +1,489 @@ +"""INVENTED populations, categories, coefficients and death probabilities.""" + +from copy import deepcopy +from dataclasses import replace + +import numpy as np +import pandas as pd +import pytest + +from populace_dynamics.baselines.base import BaselineYearAwareMortality +from populace_dynamics.engine import di_entitlement +from populace_dynamics.engine.di_entitlement import apply_di_aware_mortality +from populace_dynamics.engine.di_entitlement_rates import ( + MAX_AGE, + SINGLE_YEAR_AGE_BANDS, + DIEntitlementRates, + DIEntitlementSpec, +) +from populace_dynamics.engine.loop import PeriodContext +from populace_dynamics.engine.rng import ProjectionRNGRegistry +from populace_dynamics.engine.steps import ( + AgeSexMortalityModel, + apply_mortality, +) +from populace_dynamics.mortality_ses.model import DifferentialMortality +from populace_dynamics.mortality_ses.specification import ( + DIMENSIONS, + HoldLastCoefficientPath, + RelativeRiskSpec, + RiskDimension, + UnclassifiedValueError, +) +from populace_dynamics.mortality_ses.step import apply_differential_mortality + + +def invented_specification(risk=1.0, rule="refuse"): + dimensions = tuple( + RiskDimension( + name, + name, + ("invented_reference", "invented_higher_risk"), + "invented_reference", + ("invented_unclassified",), + ) + for name in DIMENSIONS + ) + table = pd.DataFrame( + [ + { + "year": 2020, + "dimension": dimension.name, + "category": category, + "sex": sex, + "age_lower": 0, + "age_upper": MAX_AGE, + "log_relative_risk": ( + np.log(risk) if category == "invented_higher_risk" else 0.0 + ), + } + for dimension in dimensions + for category in dimension.categories + for sex in ("female", "male") + ] + ) + return RelativeRiskSpec(dimensions, HoldLastCoefficientPath(table), rule) + + +def invented_baseline(q=0.2): + return BaselineYearAwareMortality( + "invented_baseline", + {2020: {sex: np.full(MAX_AGE + 1, q) for sex in ("female", "male")}}, + ) + + +def invented_frame(*, di=False): + frame = pd.DataFrame( + { + "person_id": [8, 2, 5, 1, 9, 3], + "age": [40] * 6, + "sex": ["female"] * 6, + "birth_year": [1979] * 6, + "year": [2019] * 6, + "weight": [1.0, 2.0, 1.0, 3.0, 1.0, 2.0], + **{name: ["invented_reference"] * 6 for name in DIMENSIONS}, + } + ) + frame.loc[[1, 2, 4], "lifetime_earnings"] = "invented_higher_risk" + if di: + frame["di_entitled"] = [True, False, False, False, False, False] + frame["di_conversion_year"] = pd.array( + [None, 2018, None, None, None, None], dtype="Int64" + ) + for column in ("di_award_year", "di_award_age", "di_recovery_year"): + frame[column] = pd.array([None] * len(frame), dtype="Int64") + frame["di_converted"] = [False, True, False, False, False, False] + frame["di_event"] = "none" + return frame + + +def invented_context(frame, *, registry=False): + return PeriodContext( + 1, + 2020, + 0, + {}, + rng_registry=ProjectionRNGRegistry(19, 25) if registry else None, + person_ordinals={ + person: index + for index, person in enumerate(sorted(frame["person_id"])) + }, + ) + + +def invented_rates(death=0.1, q=0.2): + shape = (2, MAX_AGE + 1) + return DIEntitlementRates( + spec=DIEntitlementSpec(), + incidence=np.zeros(shape), + recovery_attained=np.zeros(shape), + death_attained=np.full(shape, death), + population_reference_death=np.full(shape, q), + recovery_select=None, + death_select=None, + recovery_level_factor=1.0, + death_level_factor=1.0, + ) + + +@pytest.mark.parametrize("registry", [False, True]) +def test_identity_probabilities_and_step_deaths(registry): + frame = invented_frame() + context = invented_context(frame, registry=registry) + baseline = invented_baseline() + model = DifferentialMortality( + baseline, invented_specification(), weight_column="weight" + ) + np.testing.assert_array_equal( + model(frame, context), baseline(frame, context) + ) + assert not hasattr(model, "probabilities") + equivalent = AgeSexMortalityModel( + SINGLE_YEAR_AGE_BANDS, + { + (AgeSexMortalityModel.band_label(a, b), sex): 0.2 + for a, b in SINGLE_YEAR_AGE_BANDS + for sex in ("female", "male") + }, + ) + first_rng = np.random.default_rng(72) + second_rng = np.random.default_rng(72) + deaths, cells = {}, {} + actual = apply_differential_mortality( + frame, + context, + first_rng, + model=model, + weight_column="weight", + death_log=deaths, + cell_log=cells, + ) + expected = apply_mortality(frame, context, second_rng, model=equivalent) + pd.testing.assert_frame_equal(actual, expected) + assert first_rng.bit_generator.state == second_rng.bit_generator.state + assert set(actual["person_id"]) | set(deaths[2020]["person_id"]) == set( + frame["person_id"] + ) + assert cells[2020]["clipping_residual"].eq(0).all() + + +@pytest.mark.parametrize("registry", [False, True]) +def test_di_adapter_identity(registry): + frame = invented_frame(di=True) + context = invented_context(frame, registry=registry) + baseline = invented_baseline() + rates = invented_rates() + model = DifferentialMortality( + baseline, + invented_specification(), + weight_column="weight", + di_configuration=rates, + ) + actual_log, expected_log = {}, {} + actual_cells, expected_cells = {}, {} + first_rng, second_rng = np.random.default_rng(5), np.random.default_rng(5) + actual = apply_di_aware_mortality( + frame, + context, + first_rng, + population_model=model, + rates=rates, + weight_column="weight", + death_log=actual_log, + cell_log=actual_cells, + ) + expected = apply_di_aware_mortality( + frame, + context, + second_rng, + population_model=baseline, + rates=rates, + weight_column="weight", + death_log=expected_log, + cell_log=expected_cells, + ) + pd.testing.assert_frame_equal(actual, expected) + pd.testing.assert_frame_equal(actual_log[2020], expected_log[2020]) + pd.testing.assert_frame_equal(actual_cells[2020], expected_cells[2020]) + assert first_rng.bit_generator.state == second_rng.bit_generator.state + + +def test_model_pure_and_no_randomness(monkeypatch): + frame = invented_frame() + original = frame.copy(deep=True) + context = invented_context(frame, registry=True) + rng = np.random.default_rng(5) + state = deepcopy(rng.bit_generator.state) + + def refuse_generator(*args, **kwargs): + raise AssertionError("pure model touched a generator") + + monkeypatch.setattr(PeriodContext, "person_generator", refuse_generator) + monkeypatch.setattr(np.random, "default_rng", refuse_generator) + model = DifferentialMortality( + invented_baseline(), invented_specification(3), weight_column="weight" + ) + first, first_log = model.evaluate(frame, context) + second, second_log = model.evaluate(frame, context) + np.testing.assert_array_equal(first, second) + pd.testing.assert_frame_equal(first_log, second_log) + pd.testing.assert_frame_equal(frame, original) + assert rng.bit_generator.state == state + + +def test_age_check_always_precedes_cap_even_if_di_check_disabled(): + frame = invented_frame(di=True) + context = invented_context(frame) + model = DifferentialMortality( + invented_baseline(), + invented_specification(), + weight_column=None, + di_configuration=DIEntitlementSpec(validate_age_column=False), + ) + frame.loc[0, "age"] += 1 + with pytest.raises(ValueError, match="context.year - 1 - birth_year"): + model(frame, context) + frame["age"] = 130 + frame["birth_year"] = 1889 + _, cells = model.evaluate(frame, context) + assert cells["age"].tolist() == [120] + + +@pytest.mark.parametrize("post_conversion", ["di_origin", "population"]) +def test_fixed_origin_matches_configured_adapter_rule(post_conversion): + frame = invented_frame(di=True) + context = invented_context(frame) + model = DifferentialMortality( + invented_baseline(), + invented_specification(3), + weight_column="weight", + di_configuration=DIEntitlementSpec( + post_conversion_mortality=post_conversion + ), + ) + q, cells = model.evaluate(frame, context) + assert q[0] == 0.2 + if post_conversion == "di_origin": + assert q[1] == 0.2 + assert cells["fixed_persons"].tolist() == [2] + else: + assert q[1] > 0.2 + assert cells["fixed_persons"].tolist() == [1] + assert not cells["infeasible"].any() + + +def test_composition_preserves_unclipped_ratios_and_baseline_origin( + monkeypatch, +): + frame = invented_frame(di=True) + context = invented_context(frame) + rates = invented_rates() + model = DifferentialMortality( + invented_baseline(), + invented_specification(3), + weight_column="weight", + di_configuration=rates, + ) + captured = {} + original = di_entitlement._net_of_di_origin + + def capture(q, population, fixed, sex, age, weight, bands): + out, cells = original(q, population, fixed, sex, age, weight, bands) + captured.update(population=population, fixed=fixed, out=out) + return out, cells + + monkeypatch.setattr(di_entitlement, "_net_of_di_origin", capture) + cells = {} + apply_di_aware_mortality( + frame, + context, + np.random.default_rng(6), + population_model=model, + rates=rates, + weight_column="weight", + cell_log=cells, + ) + ordered = frame.sort_values("person_id") + q_base = model.baseline(ordered, context) + fixed = captured["fixed"] + np.testing.assert_array_equal(captured["population"][fixed], q_base[fixed]) + other = ~fixed + factor = captured["out"][other] / captured["population"][other] + np.testing.assert_allclose(factor, factor[0], rtol=1e-13, atol=1e-13) + log = cells[2020] + np.testing.assert_allclose( + log["expected_deaths_applied"], + log["expected_deaths_population_model"], + rtol=1e-12, + atol=1e-12, + ) + + +def test_adapter_clipping_residual_is_derived_from_its_logged_totals( + monkeypatch, +): + """Intended: the unchanged adapter clips once after SES normalization.""" + frame = invented_frame(di=True) + context = invented_context(frame) + rates = invented_rates(death=0.1, q=0.6) + model = DifferentialMortality( + invented_baseline(q=0.6), + invented_specification(3), + weight_column="weight", + di_configuration=rates, + ) + captured = {} + original = di_entitlement._net_of_di_origin + + def capture(q, population, fixed, sex, age, weight, bands): + out, cells = original(q, population, fixed, sex, age, weight, bands) + captured.update(out=out, population=population, weight=weight) + return out, cells + + monkeypatch.setattr(di_entitlement, "_net_of_di_origin", capture) + cells = {} + apply_di_aware_mortality( + frame, + context, + np.random.default_rng(6), + population_model=model, + rates=rates, + weight_column="weight", + cell_log=cells, + ) + logged = cells[2020] + residual = ( + logged["expected_deaths_population_model"] + - logged["expected_deaths_applied"] + ) + assert not logged["infeasible"].any() + assert residual.gt(0).all() + assert (captured["out"] == 1).any() + np.testing.assert_allclose( + np.dot(captured["weight"], captured["out"]), + np.dot(captured["weight"], captured["population"]) - residual.sum(), + rtol=1e-12, + atol=1e-12, + ) + + +def test_unclassified_cell_counts_on_reference_and_refusal(): + frame = invented_frame() + frame.loc[[0, 1], "education"] = "invented_unclassified" + model = DifferentialMortality( + invented_baseline(), invented_specification(), weight_column="weight" + ) + with pytest.raises(UnclassifiedValueError) as raised: + model(frame, invented_context(frame)) + assert raised.value.cell_log["unclassified_refuse"].sum() == 2 + reference = replace( + model, + specification=replace( + model.specification, unclassified_rule="reference" + ), + ) + _, cells = reference.evaluate(frame, invented_context(frame)) + assert cells["unclassified_reference"].sum() == 2 + + +@pytest.mark.parametrize("age", [-1, 1.5, np.nan, np.inf]) +def test_invalid_age_refused(age): + frame = invented_frame().drop(columns="birth_year") + frame["age"] = frame["age"].astype(float) + frame.loc[0, "age"] = age + model = DifferentialMortality( + invented_baseline(), invented_specification(), weight_column=None + ) + with pytest.raises(ValueError): + model(frame, invented_context(frame)) + + +def test_missing_columns_and_ambiguous_di_refused(): + frame = invented_frame() + model = DifferentialMortality( + invented_baseline(), invented_specification(), weight_column="weight" + ) + for column in ("age", "sex", "weight", "education"): + with pytest.raises(ValueError, match="missing"): + model(frame.drop(columns=column), invented_context(frame)) + with pytest.raises(ValueError, match="DI columns require"): + model(invented_frame(di=True), invented_context(frame)) + configured = replace(model, di_configuration=DIEntitlementSpec()) + with pytest.raises(ValueError, match="missing columns"): + configured(frame, invented_context(frame)) + + +@pytest.mark.parametrize( + "column,value", + [ + ("di_entitled", None), + ("di_entitled", "yes"), + ("di_conversion_year", 2018.5), + ], +) +def test_invalid_di_state_refused(column, value): + frame = invented_frame(di=True) + frame[column] = frame[column].astype(object) + frame.loc[0, column] = value + model = DifferentialMortality( + invented_baseline(), + invented_specification(), + weight_column=None, + di_configuration=DIEntitlementSpec(), + ) + with pytest.raises(ValueError): + model(frame, invented_context(frame)) + + +def test_model_probability_attribute_and_wrong_interfaces_refused(): + class InventedYearlessBaseline: + bands = SINGLE_YEAR_AGE_BANDS + probabilities = None + + def __call__(self, frame, context): + return np.zeros(len(frame)) + + with pytest.raises(ValueError, match="probabilities attribute"): + DifferentialMortality( + InventedYearlessBaseline(), + invented_specification(), + weight_column=None, + ) + + class InventedBadSubclass(DifferentialMortality): + probabilities = None + + with pytest.raises(ValueError, match="probabilities attribute"): + InventedBadSubclass( + invented_baseline(), invented_specification(), weight_column=None + ) + with pytest.raises(TypeError): + DifferentialMortality(invented_baseline(), invented_specification()) + + class InventedWrongBands: + bands = ((0, 120),) + + def __call__(self, frame, context): + return np.zeros(len(frame)) + + with pytest.raises(ValueError, match="SINGLE_YEAR"): + DifferentialMortality( + InventedWrongBands(), + invented_specification(), + weight_column=None, + ) + + +def test_step_weight_must_match_model(): + frame = invented_frame() + model = DifferentialMortality( + invented_baseline(), invented_specification(), weight_column="weight" + ) + with pytest.raises(ValueError, match="weight_column must match"): + apply_differential_mortality( + frame, + invented_context(frame), + np.random.default_rng(1), + model=model, + weight_column=None, + ) diff --git a/tests/mortality_ses/test_normalizer.py b/tests/mortality_ses/test_normalizer.py new file mode 100644 index 00000000..daaee1ce --- /dev/null +++ b/tests/mortality_ses/test_normalizer.py @@ -0,0 +1,264 @@ +"""Cell accounting tests using invented probabilities, risks and weights.""" + +import numpy as np +import pandas as pd +import pytest + +from populace_dynamics.engine.di_entitlement import _net_of_di_origin +from populace_dynamics.mortality_ses.normalizer import ( + CELL_LOG_COLUMNS, + normalize_probabilities, +) + + +def normalize( + base, risks, fixed=None, prescribed=None, weights=None, cells=None +): + size = len(base) + return normalize_probabilities( + np.asarray(base), + np.asarray(risks), + np.zeros(size, dtype=bool) if fixed is None else np.asarray(fixed), + np.asarray(base if prescribed is None else prescribed), + [(2030, "male", 70)] * size if cells is None else cells, + np.ones(size) if weights is None else np.asarray(weights), + ) + + +def test_invented_uncapped_cell_conserves_and_logs_counts(): + q, log = normalize_probabilities( + np.array([0.2, 0.2, 0.2]), + np.array([1.0, 2.0, 3.0]), + np.array([True, False, False]), + np.array([0.2, 0.2, 0.2]), + [(2030, "female", 70)] * 3, + np.array([1.0, 2.0, 1.0]), + unclassified_counts={"unclassified_reference": np.array([0, 1, 2])}, + ) + np.testing.assert_allclose(q, [0.2, 0.2 * 2 * 3 / 7, 0.2 * 3 * 3 / 7]) + assert log.columns.tolist() == list(CELL_LOG_COLUMNS) + cell = log.iloc[0] + assert (cell.year, cell.sex, cell.age) == (2030, "female", 70) + assert cell.persons == 3 + assert cell.fixed_persons == 1 + assert cell.weight == 4 + assert cell.fixed_weight == 1 + assert cell.baseline_expected_deaths == pytest.approx(0.8) + assert cell.applied_expected_deaths == pytest.approx(0.8) + assert cell.factor == pytest.approx(3 / 7) + assert cell.clipped_persons == 0 + assert cell.clipping_residual == 0 + assert not cell.infeasible + assert cell.unclassified_reference == 3 + assert cell.unclassified_refuse == 0 + + +def test_iterative_capping_reallocates_invented_remaining_deaths(): + q, log = normalize([0.9, 0.9, 0.9], [1, 10, 100]) + np.testing.assert_allclose(q, [0.7, 1, 1], rtol=1e-12, atol=1e-12) + assert log.iloc[0].clipped_persons == 2 + assert log.iloc[0].factor == pytest.approx(7 / 9) + assert log.iloc[0].clipping_residual == pytest.approx(0, abs=1e-12) + assert not log.iloc[0].infeasible + + +def test_fixed_deaths_above_invented_target_are_infeasible(): + q, log = normalize( + [0.1, 0.1], [1, 3], fixed=[True, False], prescribed=[0.5, 0.1] + ) + np.testing.assert_array_equal(q, [0.5, 0]) + assert log.iloc[0].infeasible + assert log.iloc[0].clipping_residual == pytest.approx(-0.3) + + +def test_invented_target_beyond_free_positive_support_is_logged(): + q, log = normalize( + [1, 0.1, 0], + [1, 2, 10], + fixed=[True, False, False], + prescribed=[0, 0.1, 0], + ) + np.testing.assert_array_equal(q, [0, 1, 0]) + assert log.iloc[0].infeasible + assert log.iloc[0].clipping_residual == pytest.approx(0.1) + + +def test_equal_risks_preserve_baseline_bits_with_fixed_baseline_rows(): + baseline = np.array([0.0, 0.003, 0.25, 1.0]) + for risks in ([1] * 4, [1e-308] * 4, [1e308] * 4): + q, log = normalize( + baseline, + risks, + fixed=[True, False, False, True], + weights=[0.0, 2.0, 7.0, 0.25], + ) + np.testing.assert_array_equal( + q.view(np.uint64), baseline.view(np.uint64) + ) + assert not log.iloc[0].infeasible + assert log.iloc[0].clipping_residual == 0 + + +def test_finite_extreme_risks_do_not_underflow_positive_support(): + q, log = normalize([0.9, 0.9], [1e-308, 1e308]) + np.testing.assert_allclose(q, [0.8, 1], rtol=1e-12, atol=1e-12) + assert not log.iloc[0].infeasible + assert np.isfinite(log.iloc[0].factor) + assert log.iloc[0].clipping_residual == pytest.approx(0, abs=1e-12) + + +def test_extreme_invented_products_retain_positive_probability_support(): + # The first free row's q*RR underflows float64, but its probability is + # needed after the second free row caps. Fixed rows fund the target. + q, log = normalize( + [0.8, 1e-300, 0.8], + [1, 1e-300, 1e300], + fixed=[True, False, False], + prescribed=[0, 1e-300, 0.8], + ) + np.testing.assert_allclose(q, [0, 0.6, 1], rtol=1e-12, atol=1e-12) + assert not log.iloc[0].infeasible + assert log.iloc[0].clipping_residual == pytest.approx(0, abs=1e-12) + + +def test_zero_baselines_and_zero_weights_remain_bounded(): + q, log = normalize([0, 0.4, 0.4], [9, 1, 3], weights=[1, 0, 1]) + np.testing.assert_allclose(q, [0, 0.4 / 3, 0.4]) + assert log.iloc[0].baseline_expected_deaths == pytest.approx(0.4) + assert not log.iloc[0].infeasible + zero_q, zero_log = normalize([0.2, 0.6], [1, 2], weights=[0, 0]) + assert np.all((zero_q >= 0) & (zero_q <= 1)) + assert zero_log.iloc[0].weight == 0 + assert zero_log.iloc[0].clipping_residual == 0 + + +def test_cells_are_separate_and_sorted_and_inputs_are_unchanged(): + baseline = np.array([0.3, 0.1, 0.3, 0.1]) + risks = np.array([1.0, 1.0, 2.0, 2.0]) + fixed = np.zeros(4, dtype=bool) + weights = np.array([2.0, 1.0, 2.0, 1.0]) + cells = [(2031, "male", 71), (2030, "female", 70)] * 2 + before = [x.copy() for x in (baseline, risks, fixed, weights)] + q, log = normalize_probabilities( + baseline, risks, fixed, baseline, cells, weights + ) + np.testing.assert_allclose(q, [0.2, 1 / 15, 0.4, 2 / 15]) + assert log[["year", "sex", "age"]].values.tolist() == [ + [2030, "female", 70], + [2031, "male", 71], + ] + for actual, original in zip( + (baseline, risks, fixed, weights), before, strict=True + ): + np.testing.assert_array_equal(actual, original) + + +def test_empty_input_has_declared_log_schema(): + q, log = normalize([], []) + assert q.shape == (0,) + assert log.empty + assert list(log.columns) == list(CELL_LOG_COLUMNS) + + +@pytest.mark.parametrize("bad_risk", [0, -1, np.nan, np.inf, -np.inf]) +def test_refuses_nonpositive_or_nonfinite_invented_risks(bad_risk): + with pytest.raises(ValueError, match="relative risks"): + normalize([0.2], [bad_risk]) + + +@pytest.mark.parametrize( + ("field", "value", "message"), + [ + ("q_base", [1.1], "baseline probabilities"), + ("fixed_probabilities", [-0.1], "fixed probabilities"), + ("weights", [-1], "weights"), + ("weights", [np.inf], "weights"), + ("fixed_mask", [1], "boolean"), + ("cell_ids", [], "cell_ids"), + ("cell_ids", [(2030, "male")], "triples"), + ("cell_ids", [(2030, None, 70)], "missing"), + ("relative_risk", [[1]], "one-dimensional"), + ], +) +def test_refuses_invalid_accounting_inputs(field, value, message): + arguments = dict( + q_base=[0.2], + relative_risk=[1], + fixed_mask=np.array([False]), + fixed_probabilities=[0.2], + cell_ids=[(2030, "male", 70)], + weights=[1], + ) + arguments[field] = value + with pytest.raises(ValueError, match=message): + normalize_probabilities(**arguments) + + +@pytest.mark.parametrize( + "counts", [{"weight": [1]}, {"missing": [-1]}, {"missing": [0.5]}] +) +def test_refuses_invalid_unclassified_counts(counts): + with pytest.raises(ValueError, match="count"): + normalize_probabilities( + [0.2], + [1], + np.array([False]), + [0.2], + [(2030, "male", 70)], + [1], + unclassified_counts=counts, + ) + + +def test_matches_di_netting_arithmetic_without_capping(): + baseline = np.array([0.2, 0.2, 0.2, 0.1, 0.1]) + fixed = np.array([True, False, False, True, False]) + prescribed = np.array([0.3, 0.2, 0.2, 0.05, 0.1]) + weights = np.array([1, 2, 3, 2, 4]) + sexes = np.array([0, 0, 0, 1, 1]) + ages = np.array([70, 70, 70, 71, 71]) + actual, log = normalize( + baseline, + np.ones(5), + fixed=fixed, + prescribed=prescribed, + weights=weights, + cells=list(zip([2030] * 5, sexes, ages, strict=True)), + ) + expected, _ = _net_of_di_origin( + prescribed, baseline, fixed, sexes, ages, weights, [(70, 70), (71, 71)] + ) + np.testing.assert_allclose(actual, expected, rtol=1e-12, atol=1e-12) + assert not log.infeasible.any() + assert not log.clipped_persons.any() + + +def test_intended_iterative_capping_differs_from_adapters_single_clip(): + # Invented heterogeneous baselines make the adapter's single clip bind. + baseline = np.array([0.9, 0.9, 0.1]) + fixed = np.array([True, False, False]) + prescribed = np.array([0.1, 0.9, 0.1]) + actual, log = normalize( + baseline, [1, 1, 1], fixed=fixed, prescribed=prescribed + ) + expected, _ = _net_of_di_origin( + prescribed, + baseline, + fixed, + np.zeros(3, dtype=int), + np.full(3, 70), + np.ones(3), + [(70, 70)], + ) + np.testing.assert_allclose(actual, [0.1, 1, 0.8], atol=1e-12) + assert np.sum(actual) == pytest.approx(np.sum(baseline), abs=1e-12) + assert np.sum(expected) < np.sum(baseline) + assert log.iloc[0].clipped_persons == 1 + + +def test_repeat_calls_are_deterministic_in_probabilities_and_log(): + arguments = ([0.8, 0.8, 0.8], [1, 3, 20]) + first_q, first_log = normalize(*arguments) + second_q, second_log = normalize(*arguments) + np.testing.assert_array_equal(first_q, second_q) + pd.testing.assert_frame_equal(first_log, second_log) diff --git a/tests/mortality_ses/test_properties.py b/tests/mortality_ses/test_properties.py new file mode 100644 index 00000000..d3d93ea2 --- /dev/null +++ b/tests/mortality_ses/test_properties.py @@ -0,0 +1,762 @@ +"""Accounting properties using only invented populations and coefficients. + +All categories, coefficients, probability tables and DI rates below are +INVENTED. These tests establish mechanism identities, not an empirical +mortality gradient. Conservation uses absolute and relative tolerance 1e-12 +for finite double-precision products, sums and the adapter's cell arithmetic. +""" + +from __future__ import annotations + +from copy import deepcopy +from dataclasses import replace +from unittest.mock import patch + +import numpy as np +import pandas as pd +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from populace_dynamics.baselines.base import BaselineYearAwareMortality +from populace_dynamics.engine.di_entitlement import ( + _net_of_di_origin, + apply_di_aware_mortality, +) +from populace_dynamics.engine.di_entitlement_rates import ( + SINGLE_YEAR_AGE_BANDS, + DIEntitlementRates, + DIEntitlementSpec, +) +from populace_dynamics.engine.loop import PeriodContext +from populace_dynamics.engine.rng import ( + ProjectionModule, + ProjectionRNGRegistry, +) +from populace_dynamics.engine.steps import apply_mortality +from populace_dynamics.mortality_ses.model import DifferentialMortality +from populace_dynamics.mortality_ses.normalizer import normalize_probabilities +from populace_dynamics.mortality_ses.specification import ( + DIMENSIONS, + HoldLastCoefficientPath, + RelativeRiskSpec, + RiskDimension, + UnclassifiedValueError, +) +from populace_dynamics.mortality_ses.step import apply_differential_mortality + +RTOL = ATOL = 1e-12 +PROPERTY = settings(max_examples=35, deadline=None) +REF = "invented_ref" +OTHER = "invented_other" +UNCLASSIFIED = "invented_unclassified" +SEXES = ("female", "male") +POSITIVE = st.floats(0.05, 20, allow_nan=False, allow_infinity=False) +PROBABILITY = st.floats(0, 1, allow_nan=False, allow_infinity=False) + + +def invented_dimensions(): + return tuple( + RiskDimension(name, name, (REF, OTHER), REF, (UNCLASSIFIED,)) + for name in DIMENSIONS + ) + + +def invented_table(*, years=(2020,), risk=2.0): + """Invented full-coverage age bands, with zero reference coefficients.""" + return pd.DataFrame( + [ + dict( + year=year, + dimension=name, + category=category, + sex=sex, + age_lower=0, + age_upper=120, + log_relative_risk=(0.0 if category == REF else np.log(risk)), + ) + for year in years + for name in DIMENSIONS + for category in (REF, OTHER) + for sex in SEXES + ] + ) + + +def invented_spec(*, risk=2.0, rule="refuse", path=None): + return RelativeRiskSpec( + invented_dimensions(), + path or HoldLastCoefficientPath(invented_table(risk=risk)), + unclassified_rule=rule, + ) + + +def invented_baseline(q=0.1, *, years=(2021,)): + return BaselineYearAwareMortality( + "invented_sm1", + {year: {sex: np.full(121, q) for sex in SEXES} for year in years}, + ) + + +def invented_frame(size=8, *, year=2021, di=False): + frame = pd.DataFrame( + { + "person_id": np.arange(size, 0, -1), + "age": np.full(size, 50), + "birth_year": np.full(size, year - 1 - 50), + "year": np.full(size, year - 1), + "sex": [SEXES[i % 2] for i in range(size)], + "weight": np.ones(size), + **{ + name: [REF if i % 3 else OTHER for i in range(size)] + for name in DIMENSIONS + }, + } + ) + if di: + frame["di_entitled"] = [i % 6 == 0 for i in range(size)] + frame["di_conversion_year"] = pd.array( + [year - 2 if i % 6 == 1 else pd.NA for i in range(size)], + dtype="Int64", + ) + for column in ("di_award_year", "di_award_age", "di_recovery_year"): + frame[column] = pd.array([pd.NA] * size, dtype="Int64") + frame["di_converted"] = False + frame["di_event"] = "invented_opening" + return frame + + +def invented_context(frame, *, seed=0, registry=False, year=2021): + return PeriodContext( + 1, + year, + seed, + {}, + rng_registry=ProjectionRNGRegistry(seed, 2) if registry else None, + person_ordinals={ + person_id: index + for index, person_id in enumerate(sorted(frame["person_id"])) + }, + ) + + +def invented_rates(*, multiplier=1.5, post_conversion="di_origin"): + shape = (2, 121) + return DIEntitlementRates( + spec=DIEntitlementSpec(post_conversion_mortality=post_conversion), + incidence=np.zeros(shape), + recovery_attained=np.zeros(shape), + death_attained=np.full(shape, multiplier * 0.1), + population_reference_death=np.full(shape, 0.1), + recovery_select=None, + death_select=None, + recovery_level_factor=1.0, + death_level_factor=1.0, + ) + + +def assert_conserved(q, base, cells, weights, log): + for row in log.itertuples(index=False): + mask = np.asarray( + [cell == (row.year, row.sex, row.age) for cell in cells] + ) + target = np.sum(weights[mask] * base[mask]) + applied = np.sum(weights[mask] * q[mask]) + assert applied + row.clipping_residual == pytest.approx( + target, rel=RTOL, abs=ATOL + ) + assert row.baseline_expected_deaths == pytest.approx( + target, rel=RTOL, abs=ATOL + ) + assert row.applied_expected_deaths == pytest.approx( + applied, rel=RTOL, abs=ATOL + ) + if not row.clipped_persons and not row.infeasible: + assert row.clipping_residual == 0 + + +@st.composite +def invented_normalizer_inputs(draw): + size = draw(st.integers(1, 16)) + base = np.asarray( + draw(st.lists(PROBABILITY, min_size=size, max_size=size)) + ) + risk = np.asarray(draw(st.lists(POSITIVE, min_size=size, max_size=size))) + weights = np.asarray( + draw(st.lists(POSITIVE, min_size=size, max_size=size)) + ) + fixed = np.asarray( + draw(st.lists(st.booleans(), min_size=size, max_size=size)) + ) + groups = draw(st.lists(st.integers(0, 3), min_size=size, max_size=size)) + cells = [(2021, SEXES[group % 2], 50 + group // 2) for group in groups] + return base, risk, fixed, cells, weights + + +@PROPERTY +@given(invented_normalizer_inputs()) +def test_cell_conservation_bounds_and_baseline_fixed_feasibility(inputs): + base, risk, fixed, cells, weights = inputs + q, log = normalize_probabilities(base, risk, fixed, base, cells, weights) + assert ((0 <= q) & (q <= 1)).all() + np.testing.assert_array_equal(q[fixed], base[fixed]) + assert not log["infeasible"].any() + assert_conserved(q, base, cells, weights, log) + + +@PROPERTY +@given(invented_normalizer_inputs()) +def test_identity_is_bit_for_bit_and_normalizer_is_deterministic(inputs): + base, _, fixed, cells, weights = inputs + q, log = normalize_probabilities( + base, np.ones(len(base)), fixed, base, cells, weights + ) + repeated, repeated_log = normalize_probabilities( + base, np.ones(len(base)), fixed, base, cells, weights + ) + np.testing.assert_array_equal(q, base) + np.testing.assert_array_equal(q, repeated) + pd.testing.assert_frame_equal(log, repeated_log) + + +@PROPERTY +@given(invented_normalizer_inputs(), POSITIVE, st.integers(0, 10000)) +def test_permutation_and_common_weight_scaling_preserve_probabilities( + inputs, scale, seed +): + base, risk, fixed, cells, weights = inputs + q, _ = normalize_probabilities(base, risk, fixed, base, cells, weights) + permutation = np.random.default_rng(seed).permutation(len(base)) + permuted, _ = normalize_probabilities( + base[permutation], + risk[permutation], + fixed[permutation], + base[permutation], + [cells[i] for i in permutation], + weights[permutation], + ) + scaled, _ = normalize_probabilities( + base, risk, fixed, base, cells, weights * scale + ) + np.testing.assert_array_equal(permuted, q[permutation]) + np.testing.assert_allclose(scaled, q, rtol=RTOL, atol=ATOL) + + +@PROPERTY +@given(PROBABILITY, POSITIVE, st.integers(1, 15)) +def test_one_free_risk_group_returns_baseline_probability(q_base, risk, size): + base = np.full(size, q_base) + q, _ = normalize_probabilities( + base, + np.full(size, risk), + np.arange(size) % 3 == 0, + base, + [(2021, "female", 50)] * size, + np.arange(1, size + 1, dtype=float), + ) + np.testing.assert_array_equal(q, base) + + +@PROPERTY +@given(st.floats(0.01, 0.4), st.floats(0.01, 0.4), POSITIVE) +def test_fixed_deaths_above_target_are_flagged(base_q, excess, weight): + base = np.full(2, base_q) + prescribed = np.asarray([min(1.0, 2 * base_q + excess), base_q]) + q, log = normalize_probabilities( + base, + np.ones(2), + np.asarray([True, False]), + prescribed, + [(2021, "female", 50)] * 2, + np.full(2, weight), + ) + assert log.loc[0, "infeasible"] + assert q[0] == prescribed[0] + assert q[1] == 0 + assert log.loc[0, "clipping_residual"] < 0 + + +@PROPERTY +@given(st.floats(0.02, 0.15), st.floats(0.5, 1.75), POSITIVE) +def test_unit_risk_reproduces_adapter_arithmetic_without_clipping( + baseline_q, multiplier, weight_scale +): + base = np.full(8, baseline_q) + fixed = np.arange(8) % 4 == 0 + prescribed = base.copy() + prescribed[fixed] *= multiplier + sex = np.arange(8) % 2 + age = np.full(8, 50) + weights = np.arange(1, 9, dtype=float) * weight_scale + cells = [(2021, SEXES[s], 50) for s in sex] + q, log = normalize_probabilities( + base, np.ones(8), fixed, prescribed, cells, weights + ) + old, _ = _net_of_di_origin( + prescribed, base, fixed, sex, age, weights, SINGLE_YEAR_AGE_BANDS + ) + assert not log["clipped_persons"].any() + np.testing.assert_allclose(q, old, rtol=RTOL, atol=ATOL) + + +def test_iterative_capping_intentionally_differs_from_adapter_one_pass(): + """Intended difference: reallocate the one-pass clipping shortfall.""" + base = np.asarray([0.8, 0.8, 0.1]) + fixed = np.asarray([True, False, False]) + prescribed = np.asarray([0.0, 0.8, 0.1]) + weights = np.ones(3) + cells = [(2021, "female", 50)] * 3 + q, log = normalize_probabilities( + base, np.ones(3), fixed, prescribed, cells, weights + ) + old, _ = _net_of_di_origin( + prescribed, + base, + fixed, + np.zeros(3, dtype=int), + np.full(3, 50), + weights, + SINGLE_YEAR_AGE_BANDS, + ) + np.testing.assert_allclose(q, [0, 1, 0.7], rtol=RTOL, atol=ATOL) + assert old.sum() < base.sum() + assert not np.array_equal(q, old) + assert_conserved(q, base, cells, weights, log) + + +class InventedYearlessBaseline: + """Equivalent invented baseline for the engine's yearless interface.""" + + def __init__(self, baseline, context): + self.baseline = baseline + self.context = context + + def probabilities(self, frame): + return self.baseline(frame, self.context) + + +@PROPERTY +@given(st.integers(0, 10000), st.booleans(), PROBABILITY) +def test_step_unit_risks_draw_identical_deaths_with_or_without_registry( + seed, registry, q_base +): + frame = invented_frame(12) + context = invented_context(frame, seed=seed, registry=registry) + baseline = invented_baseline(q_base) + model = DifferentialMortality( + baseline, invented_spec(risk=1), weight_column=None + ) + np.testing.assert_array_equal( + model(frame, context), baseline(frame, context) + ) + expected_rng = np.random.default_rng(seed) + applied_rng = np.random.default_rng(seed) + expected = apply_mortality( + frame, + context, + expected_rng, + model=InventedYearlessBaseline(baseline, context), + ) + deaths, cells = {}, {} + actual = apply_differential_mortality( + frame, + context, + applied_rng, + model=model, + weight_column=None, + death_log=deaths, + cell_log=cells, + ) + pd.testing.assert_frame_equal(actual, expected) + assert actual["person_id"].is_monotonic_increasing + assert applied_rng.bit_generator.state == expected_rng.bit_generator.state + assert set(deaths[2021]["person_id"]) == set(frame["person_id"]) - set( + actual["person_id"] + ) + assert not cells[2021]["infeasible"].any() + + +@PROPERTY +@given( + st.integers(0, 10000), + st.booleans(), + st.sampled_from(("di_origin", "population")), + st.floats(0.02, 0.4), +) +def test_di_adapter_unit_risks_draw_identical_deaths( + seed, registry, post_conversion, q_base +): + frame = invented_frame(24, di=True) + context = invented_context(frame, seed=seed, registry=registry) + baseline = invented_baseline(q_base) + rates = invented_rates(post_conversion=post_conversion) + model = DifferentialMortality( + baseline, + invented_spec(risk=1), + weight_column="weight", + di_configuration=rates, + ) + np.testing.assert_array_equal( + model(frame, context), baseline(frame, context) + ) + expected_deaths, expected_cells = {}, {} + actual_deaths, actual_cells = {}, {} + expected_rng = np.random.default_rng(seed) + actual_rng = np.random.default_rng(seed) + kwargs = dict(rates=rates, weight_column="weight") + expected = apply_di_aware_mortality( + frame, + context, + expected_rng, + population_model=baseline, + death_log=expected_deaths, + cell_log=expected_cells, + **kwargs, + ) + actual = apply_di_aware_mortality( + frame, + context, + actual_rng, + population_model=model, + death_log=actual_deaths, + cell_log=actual_cells, + **kwargs, + ) + pd.testing.assert_frame_equal(actual, expected) + pd.testing.assert_frame_equal(actual_deaths[2021], expected_deaths[2021]) + pd.testing.assert_frame_equal(actual_cells[2021], expected_cells[2021]) + assert actual_rng.bit_generator.state == expected_rng.bit_generator.state + + +class InventedZeroUniforms: + """Expose every positive applied probability through the death log.""" + + def random(self, size): + return np.zeros(size) + + +@PROPERTY +@given( + st.floats(0.02, 0.15), + st.floats(1.0, 2.5), + st.floats(1.0, 1.75), + st.sampled_from(("di_origin", "population")), + POSITIVE, +) +def test_composition_conserves_cells_preserves_ratios_and_fixes_di_q( + q_base, risk, multiplier, post_conversion, weight_scale +): + # Each cell has one entitled, one converted, and four other persons. + frame = pd.concat( + [ + invented_frame(6, di=True).assign(sex=sex, age=age) + for sex in SEXES + for age in (50, 70) + ], + ignore_index=True, + ) + frame["person_id"] = np.arange(len(frame)) + frame["birth_year"] = 2020 - frame["age"] + frame["weight"] = weight_scale + frame["education"] = REF + frame["race_hispanic_origin"] = REF + context = invented_context(frame) + baseline = invented_baseline(q_base) + rates = invented_rates( + multiplier=multiplier, post_conversion=post_conversion + ) + model = DifferentialMortality( + baseline, + invented_spec(risk=risk), + weight_column="weight", + di_configuration=rates, + ) + q_population, mechanism_cells = model.evaluate(frame, context) + deaths, cells = {}, {} + apply_di_aware_mortality( + frame, + context, + InventedZeroUniforms(), + population_model=model, + rates=rates, + weight_column="weight", + death_log=deaths, + cell_log=cells, + ) + logged = deaths[2021].set_index("person_id").loc[frame["person_id"]] + q_applied = logged["q_applied"].to_numpy() + fixed = frame["di_entitled"].to_numpy() | ( + frame["di_conversion_year"].notna().to_numpy() + & (post_conversion == "di_origin") + ) + np.testing.assert_array_equal( + q_population[fixed], baseline(frame, context)[fixed] + ) + np.testing.assert_array_equal( + logged["q_population"].to_numpy()[fixed], + baseline(frame, context)[fixed], + ) + assert not cells[2021]["infeasible"].any() + for row in mechanism_cells.itertuples(index=False): + cell = (frame["sex"] == row.sex) & (frame["age"] == row.age) + assert np.sum( + frame.loc[cell, "weight"] * q_applied[cell] + ) == pytest.approx( + row.baseline_expected_deaths - row.clipping_residual, + rel=RTOL, + abs=ATOL, + ) + free = cell.to_numpy() & ~fixed + assert np.all(q_applied[free] < 1) + ratios_before = q_population[free] / q_population[free][0] + ratios_after = q_applied[free] / q_applied[free][0] + np.testing.assert_allclose( + ratios_after, ratios_before, rtol=RTOL, atol=ATOL + ) + + +@PROPERTY +@given(st.floats(1, 20), st.sampled_from(("lifetime_earnings", "education"))) +def test_ordered_earnings_and_education_risks_preserve_probability_order( + risk, dimension +): + frame = invented_frame(10).assign(sex="female") + for name in DIMENSIONS: + frame[name] = REF + frame[dimension] = [REF] * 5 + [OTHER] * 5 + model = DifferentialMortality( + invented_baseline(0.5), invented_spec(risk=risk), weight_column=None + ) + q, _ = model.evaluate(frame, invented_context(frame)) + assert q[:5].max() <= q[5:].min() + + +@PROPERTY +@given(POSITIVE, st.integers(0, 10000)) +def test_model_permutation_scaling_and_no_randomness(scale, seed): + frame = invented_frame(12) + frame["weight"] = np.arange(1, 13, dtype=float) + context = invented_context(frame, seed=seed, registry=True) + model = DifferentialMortality( + invented_baseline(), invented_spec(), weight_column="weight" + ) + rng = np.random.default_rng(seed) + rng_before = deepcopy(rng.bit_generator.state) + person_rng = context.person_generator(ProjectionModule.MORTALITY, 1) + person_before = deepcopy(person_rng.bit_generator.state) + with ( + patch.object( + PeriodContext, + "person_generator", + side_effect=AssertionError("model requested a person generator"), + ), + patch( + "numpy.random.default_rng", + side_effect=AssertionError("model constructed a generator"), + ), + ): + q, log = model.evaluate(frame, context) + repeated, repeated_log = model.evaluate(frame, context) + permutation = np.random.default_rng(seed).permutation(len(frame)) + permuted = model(frame.iloc[permutation], context) + scaled = model(frame.assign(weight=frame["weight"] * scale), context) + np.testing.assert_array_equal(q, repeated) + pd.testing.assert_frame_equal(log, repeated_log) + np.testing.assert_array_equal(permuted, q[permutation]) + np.testing.assert_allclose(scaled, q, rtol=RTOL, atol=ATOL) + assert rng.bit_generator.state == rng_before + assert person_rng.bit_generator.state == person_before + + +class InventedWideningPath: + """A separately registered invented path changing only named years.""" + + def __init__(self, base, changed_year): + self.base = base + self.changed_year = changed_year + + @property + def years(self): + return tuple(sorted(set(self.base.years) | {self.changed_year})) + + def coefficients(self, year): + table = self.base.coefficients(year) + if year == self.changed_year: + table["log_relative_risk"] *= 2 + return table + + +@PROPERTY +@given(st.integers(2021, 2100), st.floats(1.1, 4)) +def test_hold_last_refuse_before_first_and_pluggable_named_year(year, risk): + frame = invented_frame(6, year=year) + path = HoldLastCoefficientPath(invented_table(risk=risk)) + spec = invented_spec(path=path) + ages = frame["age"].to_numpy() + held, _ = spec.relative_risks(frame, year, ages) + observed, _ = spec.relative_risks(frame, 2020, ages) + np.testing.assert_array_equal(held, observed) + with pytest.raises(ValueError, match="precedes"): + spec.relative_risks(frame, 2019, ages) + alternative = invented_spec(path=InventedWideningPath(path, year)) + changed, _ = alternative.relative_risks(frame, year, ages) + unchanged, _ = alternative.relative_risks(frame, year + 1, ages) + np.testing.assert_allclose(changed, held**2, rtol=RTOL, atol=ATOL) + np.testing.assert_array_equal(unchanged, held) + + +@PROPERTY +@given(st.sampled_from((np.nan, np.inf, -np.inf, -1000.0, 1000.0))) +def test_invalid_coefficients_are_refused(value): + table = invented_table() + table.loc[table["category"] == OTHER, "log_relative_risk"] = value + with pytest.raises(ValueError, match="finite and positive"): + invented_spec(path=HoldLastCoefficientPath(table)) + + +@PROPERTY +@given(st.sampled_from((0.0, -1.0, np.nan, np.inf, -np.inf))) +def test_invalid_relative_risks_are_refused(value): + with pytest.raises(ValueError, match="relative risks"): + normalize_probabilities( + np.asarray([0.2]), + np.asarray([value]), + np.asarray([False]), + np.asarray([0.2]), + [(2021, "female", 50)], + np.ones(1), + ) + + +@PROPERTY +@given(st.sampled_from(("age", "sex", "weight", *DIMENSIONS))) +def test_missing_model_columns_are_refused(column): + frame = invented_frame().drop(columns=column) + model = DifferentialMortality( + invented_baseline(), invented_spec(), weight_column="weight" + ) + with pytest.raises(ValueError, match="missing"): + model(frame, invented_context(frame)) + + +@PROPERTY +@given(st.sampled_from(DIMENSIONS), st.sampled_from((None, UNCLASSIFIED))) +def test_unclassified_refusal_and_reference_counts_are_explicit( + dimension, value +): + frame = invented_frame() + frame.loc[0, dimension] = value + context = invented_context(frame) + refusing = DifferentialMortality( + invented_baseline(), invented_spec(), weight_column=None + ) + with pytest.raises(UnclassifiedValueError) as raised: + refusing(frame, context) + assert raised.value.cell_log["unclassified_refuse"].sum() == 1 + reference = DifferentialMortality( + invented_baseline(), + invented_spec(rule="reference"), + weight_column=None, + ) + q, log = reference.evaluate(frame, context) + substituted = frame.copy() + substituted.loc[0, dimension] = REF + np.testing.assert_array_equal(q, reference(substituted, context)) + assert log["unclassified_reference"].sum() == 1 + assert log[f"{dimension}_unclassified_reference"].sum() == 1 + assert log["unclassified_refuse"].sum() == 0 + + +@PROPERTY +@given(st.sampled_from(DIMENSIONS)) +def test_unknown_categories_are_refused_even_under_reference_rule(dimension): + frame = invented_frame() + frame.loc[0, dimension] = "invented_unknown" + model = DifferentialMortality( + invented_baseline(), + invented_spec(rule="reference"), + weight_column=None, + ) + with pytest.raises(ValueError, match="unknown category"): + model(frame, invented_context(frame)) + + +@PROPERTY +@given(st.integers(1, 10), st.booleans()) +def test_age_cell_mismatch_is_refused_even_when_adapter_validation_off( + offset, validate_age +): + frame = invented_frame(di=True) + frame.loc[0, "age"] += offset + spec = DIEntitlementSpec(validate_age_column=validate_age) + model = DifferentialMortality( + invented_baseline(), + invented_spec(), + weight_column=None, + di_configuration=spec, + ) + with pytest.raises(ValueError, match="before capping"): + model(frame, invented_context(frame)) + + +def test_probability_attribute_is_refused_and_model_has_context_only_api(): + frame = invented_frame() + context = invented_context(frame) + baseline = invented_baseline() + spec = invented_spec() + model = DifferentialMortality(baseline, spec, weight_column=None) + assert not hasattr(model, "probabilities") + with pytest.raises(ValueError, match="probabilities"): + DifferentialMortality( + InventedYearlessBaseline(baseline, context), + spec, + weight_column=None, + ) + + class InventedInvalidInterface(DifferentialMortality): + probabilities = True + + with pytest.raises(ValueError, match="probabilities"): + InventedInvalidInterface(baseline, spec, weight_column=None) + + +def test_capped_age_cells_match_the_single_year_di_adapter_cells(): + frame = invented_frame(6, di=True).assign(age=125, birth_year=1895) + context = invented_context(frame) + rates = invented_rates() + model = DifferentialMortality( + invented_baseline(), + invented_spec(), + weight_column=None, + di_configuration=rates, + ) + _, log = model.evaluate(frame, context) + cells = {} + apply_di_aware_mortality( + frame, + context, + np.random.default_rng(1), + population_model=model, + rates=rates, + weight_column=None, + cell_log=cells, + ) + assert set(log["age"]) == {120} + assert set(cells[2021]["band_lower"]) == {120} + assert set(cells[2021]["band_upper"]) == {120} + assert set(log["sex"]) == set(cells[2021]["sex"]) + + +def test_model_requires_explicit_weights_and_matching_di_configuration(): + with pytest.raises(TypeError): + DifferentialMortality(invented_baseline(), invented_spec()) + frame = invented_frame(di=True) + context = invented_context(frame) + model = DifferentialMortality( + invented_baseline(), invented_spec(), weight_column=None + ) + with pytest.raises(ValueError, match="di_configuration"): + model(frame, context) + configured = replace(model, di_configuration=DIEntitlementSpec()) + with pytest.raises(ValueError, match="missing"): + configured(frame.drop(columns="di_conversion_year"), context) diff --git a/tests/mortality_ses/test_specification.py b/tests/mortality_ses/test_specification.py new file mode 100644 index 00000000..10c96cbd --- /dev/null +++ b/tests/mortality_ses/test_specification.py @@ -0,0 +1,223 @@ +"""INVENTED categories and coefficients; no empirical gradient claims.""" + +from dataclasses import replace + +import numpy as np +import pandas as pd +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from populace_dynamics.mortality_ses.specification import ( + DIMENSIONS, + HoldLastCoefficientPath, + RelativeRiskSpec, + RiskDimension, + UnclassifiedValueError, +) + + +def invented_dimensions(): + return tuple( + RiskDimension( + name, + name, + ("invented_reference", "invented_higher_risk"), + "invented_reference", + ("invented_unclassified",), + ) + for name in DIMENSIONS + ) + + +def invented_table(years=(2020, 2022), risks=(2.0, 3.0)): + return pd.DataFrame( + [ + { + "year": year, + "dimension": dimension.name, + "category": category, + "sex": sex, + "age_lower": 0, + "age_upper": 120, + "log_relative_risk": ( + 0.0 if category == dimension.reference else np.log(risk) + ), + } + for year, risk in zip(years, risks, strict=True) + for dimension in invented_dimensions() + for category in dimension.categories + for sex in ("female", "male") + ] + ) + + +def invented_frame(): + return pd.DataFrame( + { + "sex": ["female", "male"], + **{ + name: ["invented_reference", "invented_higher_risk"] + for name in DIMENSIONS + }, + } + ) + + +def test_invented_bands_log_addition_and_hold_last(): + spec = RelativeRiskSpec( + invented_dimensions(), HoldLastCoefficientPath(invented_table()) + ) + frame = invented_frame() + risk, counts = spec.relative_risks(frame, 2020, np.array([0, 120])) + np.testing.assert_allclose(risk, [1, 8]) + assert not counts["unclassified_reference"].any() + np.testing.assert_array_equal( + spec.relative_risks(frame, 2022, np.array([0, 120]))[0], + spec.relative_risks(frame, 2100, np.array([0, 120]))[0], + ) + with pytest.raises(ValueError, match="precedes"): + spec.relative_risks(frame, 2019, np.array([0, 120])) + + +def test_pluggable_path_changes_only_named_year(): + primary = HoldLastCoefficientPath(invented_table()) + + class InventedAlternative: + years = primary.years + + def coefficients(self, year): + table = primary.coefficients(year) + if year == 2025: + table["log_relative_risk"] *= 2 + return table + + first = RelativeRiskSpec(invented_dimensions(), primary) + other = replace(first, path=InventedAlternative()) + for year in (2020, 2022, 2024, 2026, 2100): + np.testing.assert_array_equal( + first.relative_risks(invented_frame(), year, [20, 20])[0], + other.relative_risks(invented_frame(), year, [20, 20])[0], + ) + assert other.relative_risks(invented_frame(), 2025, [20, 20])[0][1] > ( + first.relative_risks(invented_frame(), 2025, [20, 20])[0][1] + ) + + +@given(st.sampled_from([None, "invented_unclassified"])) +@settings(deadline=None) +def test_unclassified_rules_are_counted(value): + frame = invented_frame() + frame.loc[0, DIMENSIONS[0]] = value + spec = RelativeRiskSpec( + invented_dimensions(), HoldLastCoefficientPath(invented_table()) + ) + with pytest.raises(UnclassifiedValueError) as raised: + spec.relative_risks(frame, 2020, [20, 20]) + assert raised.value.unclassified_counts["unclassified_refuse"].sum() == 1 + risk, counts = replace(spec, unclassified_rule="reference").relative_risks( + frame, 2020, [20, 20] + ) + assert risk[0] == 1 + assert counts["unclassified_reference"].sum() == 1 + + +@pytest.mark.parametrize("rule", ["refuse", "reference"]) +def test_unknown_category_always_refused(rule): + spec = RelativeRiskSpec( + invented_dimensions(), + HoldLastCoefficientPath(invented_table()), + rule, + ) + frame = invented_frame() + frame.loc[0, DIMENSIONS[0]] = "invented_unknown" + with pytest.raises(ValueError, match="unknown category"): + spec.relative_risks(frame, 2020, [20, 20]) + + +@pytest.mark.parametrize("log_risk", [np.nan, np.inf, -np.inf, 1000, -1000]) +def test_nonfinite_nonpositive_or_overflowed_risk_refused(log_risk): + table = invented_table() + table.loc[1, "log_relative_risk"] = log_risk + with pytest.raises(ValueError, match="finite and positive"): + RelativeRiskSpec(invented_dimensions(), HoldLastCoefficientPath(table)) + + +@pytest.mark.parametrize( + "mutation", ["missing", "extra", "duplicate", "reference"] +) +def test_table_coverage_and_reference_refused(mutation): + table = invented_table() + if mutation == "missing": + table = table.iloc[1:] + elif mutation == "extra": + table.loc[0, "category"] = "invented_unknown" + elif mutation == "duplicate": + table = pd.concat([table, table.iloc[[0]]]) + else: + table.loc[0, "log_relative_risk"] = 0.5 + with pytest.raises(ValueError): + RelativeRiskSpec(invented_dimensions(), HoldLastCoefficientPath(table)) + + +def test_missing_columns_and_malformed_specification_refused(): + path = HoldLastCoefficientPath(invented_table()) + spec = RelativeRiskSpec(invented_dimensions(), path) + with pytest.raises(ValueError, match="missing"): + spec.relative_risks( + invented_frame().drop(columns="education"), 2020, [20, 20] + ) + with pytest.raises(ValueError, match="three dimensions"): + RelativeRiskSpec(invented_dimensions()[:2], path) + with pytest.raises(ValueError, match="unclassified rule"): + replace(spec, unclassified_rule="silent") + with pytest.raises(ValueError, match="unique"): + replace( + invented_dimensions()[0], categories=("duplicate", "duplicate") + ) + + +def test_empty_custom_path_and_out_of_range_year_refused(): + class InventedEmptyPath: + years = (2020,) + + def coefficients(self, year): + return invented_table().drop(columns="year").iloc[:0] + + with pytest.raises(ValueError, match="coverage"): + RelativeRiskSpec(invented_dimensions(), InventedEmptyPath()) + table = invented_table() + table["year"] = 1e30 + with pytest.raises(ValueError, match="finite integers"): + HoldLastCoefficientPath(table) + + +@given( + st.floats(min_value=-100, max_value=100, allow_nan=False), + st.integers(min_value=0, max_value=120), +) +@settings(max_examples=30, deadline=None) +def test_single_age_and_bands_are_equivalent(log_risk, age): + table = invented_table(years=(2020,), risks=(1,)) + table.loc[ + table["category"] != "invented_reference", "log_relative_risk" + ] = log_risk + expanded = pd.DataFrame( + [ + {**row, "age": a} + for row in table.drop(columns=["age_lower", "age_upper"]).to_dict( + "records" + ) + for a in range(121) + ] + ) + band_spec = RelativeRiskSpec( + invented_dimensions(), HoldLastCoefficientPath(table) + ) + age_spec = RelativeRiskSpec( + invented_dimensions(), HoldLastCoefficientPath(expanded) + ) + np.testing.assert_array_equal( + band_spec.relative_risks(invented_frame(), 2020, [age, age])[0], + age_spec.relative_risks(invented_frame(), 2020, [age, age])[0], + ) From 0aad6e66f7425b50e7ff307db62a42317d245781 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 10:45:41 -0400 Subject: [PATCH 02/17] SM1: design note for the differential mortality mechanism Documents the interfaces, the cells and fixed persons, the unclassified rule, normalization and its tolerance, traceability to the lane design (docs/design/mortality_groups_financing.md, approved at 2c19f5b6 under d1922), the invariant evidence, and what stays outside SM1. Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 191 +++++++++++++++++++++++++++++++ 1 file changed, 191 insertions(+) create mode 100644 docs/design/mortality_ses_sm1.md diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md new file mode 100644 index 00000000..75a311b2 --- /dev/null +++ b/docs/design/mortality_ses_sm1.md @@ -0,0 +1,191 @@ +# SM1: differential mortality mechanism + +SM1 implements accounting conformance on invented inputs. It admits the +mechanism's identities only, with no empirical gradient claim. Authority: +the SM1 row of section 3.1, section 5.4, section 6 and section 7.1 of +design #558 at `2c19f5b6856b9da62924706067ef96f00f55012d`, approved with +its methodology by decision d1922 (23 October 2026). Q4 fixes baseline +expected deaths within age-sex-year cells, Q5 makes lifetime earnings the +primary dimension, Q6 holds the last observed gradient, and Q11 assigns +death alignment to SM. + +The package is opt-in. No runner, registered run, configuration, baseline, +engine module or gate is changed. The coefficients and categories in every +new test are explicitly invented. SM2 supplies empirical definitions and +coefficients. + +## Interfaces + +All public objects are exported from `populace_dynamics.mortality_ses`. + +- `RiskDimension(name, column, categories, reference, + unclassified_values=())` declares one of `lifetime_earnings`, `education` + and `race_hispanic_origin`. Categories are ordered, unique strings; the + reference is one of them. All three dimensions must appear in the spec. +- `RelativeRiskSpec(dimensions, path, unclassified_rule="refuse")` validates + complete category × sex × single-age coverage for every named table + year. Coefficient tables contain `dimension`, `category`, `sex`, `age` + and `log_relative_risk`. Inclusive `age_lower`/`age_upper` bands may + replace `age`; expansion must cover ages 0–120 exactly once. The reference + coefficient is zero. Each exponentiated coefficient and each combined + log-additive risk must be finite and strictly positive. +- `CoefficientPath` requires ordered integer `years` and + `coefficients(year) -> DataFrame`. `HoldLastCoefficientPath(table)` takes + an additional `year` column, holds between named observations and after + the final observation, and refuses years before the first. A custom + path can supply a registered widening alternative without mechanism + changes; every returned table is revalidated. +- `normalize_probabilities(q_base, relative_risk, fixed_mask, + fixed_probabilities, cell_ids, weights, *, unclassified_counts=None)` + returns `(probabilities, cell_log)`. Inputs are matching row vectors; + cell IDs are `(year, sex, age)` triples. Weights are finite and + nonnegative; zero-weight rows use the same factor but cannot determine + it. Inputs are never mutated. +- `DifferentialMortality(baseline, specification, *, weight_column, + di_configuration=None)` exposes single-year `bands` and + `__call__(frame, context) -> ndarray`. `evaluate(frame, context)` also + returns its cell log. `weight_column=None` explicitly means equal + person weights. The baseline must be callable with context, have + single-year bands, and have no `probabilities` attribute. Neither does + the new model: the DI adapter would otherwise call it without the year. +- `apply_differential_mortality(frame, context, rng, *, model, + weight_column, death_log=None, cell_log=None)` is the frame step for a + projection without the DI adapter. Its weight declaration must match + the model's. It draws in stable `person_id` order using exactly + `engine.steps.apply_mortality`'s batch or person-keyed mortality streams, + then returns survivors in that order. Logs are keyed by projection year. + +## Cells, fixed persons and unclassified attributes + +Age must be a finite, nonnegative integer. When `birth_year` is present, +SM1 always refuses a row unless its uncapped `age` is exactly +`context.year - 1 - birth_year`, even when the DI spec disables its own +optional check. The cell uses `min(age, 120)`, the frame's sex and +`context.year`. This agrees with the baseline's age cap and with the DI +adapter's single-year cells, including ages over 120. Without birth year, +the supplied age is authoritative. + +DI origin is exactly `di_entitled | (di_conversion_year.notna() & +(post_conversion_mortality == "di_origin"))`. The model accepts the +adapter's `DIEntitlementRates` or `DIEntitlementSpec`; callers must pass +the same setting to both. A frame with neither DI columns nor DI +configuration fixes no person. DI columns without configuration, or an +missing configured origin columns, refuse. The adapter separately requires +its complete DI state. Origin rows get the baseline q +bit for bit, so the adapter's multiplier sees the baseline probability. +Under `non_di_mortality="net_of_di_origin"`, the adapter then scales +non-origin probabilities by one cell factor, preserving their ratios +wherever clipping does not bind. + +Missing attributes and each dimension's declared `unclassified_values` +use the named rule. `refuse` is the default. Its +`UnclassifiedValueError.cell_log` records per-cell refusal counts, using +baseline probabilities for diagnostics; no applied gradient is returned. +`reference` contributes zero log risk for that dimension. Unknown ordinary +categories always refuse. Logs include aggregate `unclassified_refuse` +and `unclassified_reference` counts and separate counts for each +dimension. Counts are dimension contributions: a person unclassified in +two dimensions contributes two to the aggregate. Thus neither substitution +nor refusal is silent. + +## Normalization and tolerance + +For each cell, the target is `sum(weight * q_base)`. Fixed rows retain +their prescribed probabilities. Free rows receive +`min(1, q_base * relative_risk * factor)`: after a set reaches one, its +deaths are subtracted and the factor is recalculated on the remaining +rows. Fixed deaths above the target make the cell infeasible and give +all free rows zero. Arbitrary fixed probabilities can also leave a target +above free positive-q support; the normalizer saturates that support and +flags infeasibility. Neither case can arise with fixed rows at baseline q. + +The exact baseline path bypasses multiply/divide rounding for unit risks +and for a common free-row risk with fixed rows at baseline. Stable sums +reduce order sensitivity; Decimal arithmetic handles extreme magnitudes +that would underflow or overflow ordinary products. Probabilities returned +are float64; the log accounts for those returned values. + +The cell log contains year, sex, age, persons, fixed persons, total weight, +fixed weight, baseline expected deaths, applied expected deaths, factor, +clipped persons, signed clipping residual and infeasibility, plus the +unclassified counts. `clipping_residual = baseline - applied`. In a +feasible cell without capping it is explicitly zero; ordinary rounding +is handled by the tolerance, rather than labeled a clipping loss. + +Conservation checks use `rtol=1e-12, atol=1e-12` in expected-death units. +This allows float64 product, summation and multiply/divide roundoff and +the existing adapter's separate reductions; it is an accounting tolerance, +not an empirical admission threshold. Identity checks use exact array +equality. Permutation checks are exact where supported by stable sums; +weight scaling and adapter comparisons use the stated tolerance. + +Intended monotonicity applies among free rows with the same baseline q in +a cell, holding other dimensions constant. This is the setting supplied +by the age-sex-year baseline and exercised by the ordering properties. +Fixed DI-origin probabilities are prescribed baseline values and are +outside that gradient ordering check. + +The existing DI adapter caps once. Its feasible cells can consequently +fall short of the target; its residual is computed from +`expected_deaths_population_model - expected_deaths_applied` in its +existing cell log. The SM1 log describes normalization before that step. +Tests label the normalizer's iterative redistribution versus the +adapter's one-pass clipping as an intended difference. No shared-module +change is required. + +## Design and build traceability + +Paths below are relative to `src/populace_dynamics/mortality_ses/` for code +and `tests/mortality_ses/` for tests. + +| Design or build element | Code | Tests | +|---|---|---| +| Design 3.1 SM1: three-dimensional, normalized mechanism | `specification.py`, `normalizer.py`, `model.py` | `test_specification.py`, `test_properties.py` | +| Design 5.4 SM1: context interface, fixed DI origin, no RNG, new step | `model.py`, `step.py` | `test_model_step.py`, `test_properties.py` | +| Design 6 SM1: accounting invariants | All package modules | All four test modules; properties below | +| Design 7.1: no wiring; later change admission | Package-only implementation, historical exclusions | Exclusion tuple and transitive reachability in `tests/estimates/test_birth_evidence_artifact.py` | +| Brief 2.1: categories, coefficients, time path, refusals | `specification.py` | `test_specification.py`, path and refusal properties | +| Brief 2.2: pure normalizer and cell diagnostics | `normalizer.py` | `test_normalizer.py`, normalizer properties | +| Brief 2.3: baseline wrapper and age/DI/weight rules | `model.py` | `test_model_step.py`, model properties | +| Brief 2.4: non-DI projection step | `step.py` | Step death-identity examples and properties | +| Brief 2.5: no wiring | Diff review: only package, tests and housekeeping edits | Historical reducer reachability test | +| Brief 2.6: design note | This document | Traceability and invariant evidence below | + +## Invariant evidence + +Every generated property uses invented values, alongside explicit examples. +Hypothesis deadlines are disabled because these are accounting tests, not +performance benchmarks. Any failing generated value is rerun and minimized +before being characterized as a mechanism defect. + +| Invariant | Implementation | Example and property coverage | +|---|---|---| +| Cell conservation and zero uncapped residual | Iterative cap in `normalizer.py` | `test_invented_uncapped_cell_conserves_and_logs_counts`, `test_cell_conservation_bounds_and_baseline_fixed_feasibility` | +| Feasibility with baseline fixed q; arbitrary excess flagged | Fixed target calculation | `test_fixed_deaths_above_invented_target_are_infeasible`, `test_fixed_deaths_above_target_are_flagged`, baseline-fixed property | +| Unit-risk bit identity | Exact common-risk baseline path | `test_equal_risks_preserve_baseline_bits_with_fixed_baseline_rows`, `test_identity_is_bit_for_bit_and_normalizer_is_deterministic` | +| Same deaths with/without registry, ordinary and DI steps | `step.py` and context-only wrapper | `test_identity_probabilities_and_step_deaths`, `test_di_adapter_identity`, corresponding generated death-identity properties | +| Bounds | Capping and input validation | Conservation/bounds property and extreme-risk examples | +| DI composition, ratios and baseline origin q | Fixed origin in `model.py`, existing adapter | Composition examples and `test_composition_conserves_cells_preserves_ratios_and_fixes_di_q` | +| Intended monotonicity by earnings and education | Positive cell factor, log-additive risks | `test_ordered_earnings_and_education_risks_preserve_probability_order` | +| Row/weight invariance and one-risk cells | Stable reductions, exact common-risk path | Permutation/scaling examples and properties; `test_one_free_risk_group_returns_baseline_probability` | +| Differential adapter arithmetic; intended clipping difference | Same uncapped algebra, iterative cap | Private-adapter differential examples and properties; named intended capping divergence | +| Time path and alternative path isolation | Path protocol and hold-last default | `test_pluggable_path_changes_only_named_year`, generated hold-last/refusal property | +| Determinism and no RNG consumption | Pure `evaluate` | Repeated-call examples, generator state and patched-generator property | +| Refusals | Specification, frame and interface validation | Nonfinite/nonpositive risks, missing columns, unknown/unclassified values, age mismatch and `probabilities` attribute examples and properties | + +## What remains outside SM1 + +SM2 owes source-audited category definitions, lifetime-earnings construction, +empirical joint coefficients, input provenance and calibration, and registered +primary/alternative paths. Its admitted scope is external-target conformance, +labeled "calibrated, not validated". SM1 supplies none of those inputs. + +SM3 owes the authorized, prospectively registered mortality-gradient gate, +prediction-boundary attributes, death alignment, power checks, simultaneous +bounds and the S0.1 exposure audit. An inconclusive dimension admits nothing. +Before any Track B projection uses SM mortality, SM3 must pass and the +section 7.1 change-admission protocol must be frozen. Affected mortality, +marital, DI, joint-payment, recipient-denominator, survivor-pathway and timing +checks must be completed or rerun; earlier J/B7 evidence does not transfer. +This build changes no registered projection and admits no population, +gradient, lifetime-benefit or national financing claim. From 0c226e3f58944a5e7a7219ef575bca6aca1928cb Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 10:45:41 -0400 Subject: [PATCH 03/17] SM1: exclusion tuple and tier counts Adds the five mortality_ses sources to POST_REVIEW_SOURCE_EXCLUSIONS and its exact-tuple test, and asserts that they stay unreachable from the birth-evidence reducer. It also recounts the tiers: unit 18,088 -> 18,178, total 24,619 -> 24,709, from recount-tiers.py with the empty PSID root. Co-Authored-By: Claude Opus 5.5 --- scripts/first_estimates_birth_evidence.py | 6 ++++++ tests/README-tiers.md | 4 ++-- tests/estimates/test_birth_evidence_artifact.py | 15 +++++++++++++++ tests/tier_counts.json | 2 +- 4 files changed, 24 insertions(+), 3 deletions(-) diff --git a/scripts/first_estimates_birth_evidence.py b/scripts/first_estimates_birth_evidence.py index 83a0ca1f..c672e2bb 100644 --- a/scripts/first_estimates_birth_evidence.py +++ b/scripts/first_estimates_birth_evidence.py @@ -409,6 +409,12 @@ Path("src/populace_dynamics/group_breakdowns/fra68.py"), Path("src/populace_dynamics/group_breakdowns/uniform_cut.py"), Path("src/populace_dynamics/group_breakdowns/min_benefit.py"), + # SM1 is opt-in and absent from the historical projection call graph. + Path("src/populace_dynamics/mortality_ses/__init__.py"), + Path("src/populace_dynamics/mortality_ses/specification.py"), + Path("src/populace_dynamics/mortality_ses/normalizer.py"), + Path("src/populace_dynamics/mortality_ses/model.py"), + Path("src/populace_dynamics/mortality_ses/step.py"), ) POST_REVIEW_SHARED_SOURCE_BLOBS = { Path( diff --git a/tests/README-tiers.md b/tests/README-tiers.md index e2715c12..084d7202 100644 --- a/tests/README-tiers.md +++ b/tests/README-tiers.md @@ -39,9 +39,9 @@ pytest --collect-only -q -m oracle_policyengine | tail -1 | Tier | Tests at HEAD | |---|---:| -| `unit` | 18,088 | +| `unit` | 18,178 | | `artifact` | 4,253 | | `integration_psid` | 1,450 | | `reproduction_legacy` | 520 | | `oracle_policyengine` | 308 | -| **Total** | **24,619** | +| **Total** | **24,709** | diff --git a/tests/estimates/test_birth_evidence_artifact.py b/tests/estimates/test_birth_evidence_artifact.py index 3a0a2086..acd70aea 100644 --- a/tests/estimates/test_birth_evidence_artifact.py +++ b/tests/estimates/test_birth_evidence_artifact.py @@ -253,6 +253,11 @@ def test_post_review_sources_are_outside_historical_reducer_identity(): Path("src/populace_dynamics/group_breakdowns/fra68.py"), Path("src/populace_dynamics/group_breakdowns/uniform_cut.py"), Path("src/populace_dynamics/group_breakdowns/min_benefit.py"), + Path("src/populace_dynamics/mortality_ses/__init__.py"), + Path("src/populace_dynamics/mortality_ses/specification.py"), + Path("src/populace_dynamics/mortality_ses/normalizer.py"), + Path("src/populace_dynamics/mortality_ses/model.py"), + Path("src/populace_dynamics/mortality_ses/step.py"), ) assert reducer.POST_REVIEW_SHARED_SOURCE_BLOBS == { Path( @@ -598,6 +603,16 @@ def test_post_review_exclusions_are_unreachable_from_birth_evidence(): "birth-evidence reducer" ) assert "populace_dynamics.engine.steps" in reachable + mortality_ses_modules = { + name + for name in module_paths + if name == "populace_dynamics.mortality_ses" + or name.startswith("populace_dynamics.mortality_ses.") + } + assert len(mortality_ses_modules) == 5 + assert mortality_ses_modules.isdisjoint( + reachable + ), "SM1 became reachable from the historical birth-evidence reducer" scenario_module = "populace_dynamics.scenario_benefits" assert scenario_module in module_paths assert scenario_module not in reachable, ( diff --git a/tests/tier_counts.json b/tests/tier_counts.json index 0303dc6a..5f704c67 100644 --- a/tests/tier_counts.json +++ b/tests/tier_counts.json @@ -1,7 +1,7 @@ { "schema_version": 1, "counts": { - "unit": 18088, + "unit": 18178, "artifact": 4253, "integration_psid": 1450, "reproduction_legacy": 520, From 8f341b0f49ed1d655c7f49987e5d379f89b974aa Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 14:17:22 -0400 Subject: [PATCH 04/17] SM1: weightless free rows take the unit-weight factor, not zero When every free row with positive baseline q in a cell had zero weight, the weighted target was zero and the normalizer set the factor to zero, so those persons could never die. An invented cell with baseline 0.2 and relative risks 1 and 2 returned q = [0, 0], logged factor 0 and raised no flag; through the model, an all-zero-weight frame did the same. No factor changes such a cell's weighted deaths, so the weighted target says nothing about it. The cell now sets the factor as if each person had weight one: its free rows keep their relative risks, and their unweighted expected deaths equal their unweighted baseline total. Both arithmetic paths (float and Decimal) apply the rule. If fixed probabilities off baseline leave a nonzero weighted target, the cell stays infeasible and the residual records the gap. With unit risks this reproduces the DI adapter, which keeps the factor at one when the other rows carry no expected deaths. Tests: - examples on both arithmetic paths, mixed cells, infeasible cells of both signs, and the model; - a property that zero weights give exactly the unit-weight probabilities and factors; - a differential property against _net_of_di_origin on every unclipped cell, over zero weights and arbitrary fixed probabilities; - the normalizer properties now draw zero weights. An ad hoc search of 20,000 invented trials (53,481 unclipped cells) found 1,207 mismatches with the adapter before this change and none after. Found by the SM1 review r1, whose attempt was lost to an API error before it reported. Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 23 +++- .../mortality_ses/normalizer.py | 49 +++++++-- tests/mortality_ses/test_normalizer.py | 28 +++++ tests/mortality_ses/test_properties.py | 100 +++++++++++++++++- 4 files changed, 183 insertions(+), 17 deletions(-) diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md index 75a311b2..34b382dd 100644 --- a/docs/design/mortality_ses_sm1.md +++ b/docs/design/mortality_ses_sm1.md @@ -40,7 +40,9 @@ All public objects are exported from `populace_dynamics.mortality_ses`. returns `(probabilities, cell_log)`. Inputs are matching row vectors; cell IDs are `(year, sex, age)` triples. Weights are finite and nonnegative; zero-weight rows use the same factor but cannot determine - it. Inputs are never mutated. + it, and a cell whose positive-q free rows all weigh zero sets its + factor with unit weights (see "Normalization and tolerance"). Inputs + are never mutated. - `DifferentialMortality(baseline, specification, *, weight_column, di_configuration=None)` exposes single-year `bands` and `__call__(frame, context) -> ndarray`. `evaluate(frame, context)` also @@ -98,6 +100,22 @@ rows. Fixed deaths above the target make the cell infeasible and give all free rows zero. Arbitrary fixed probabilities can also leave a target above free positive-q support; the normalizer saturates that support and flags infeasibility. Neither case can arise with fixed rows at baseline q. +Both rules apply when that support carries weight; the weightless case +follows. + +Weights decide the factor only when they can. If every free row with +positive baseline q has zero weight, no factor changes the cell's weighted +deaths, so the weighted target says nothing about it. The cell then sets +the factor as if each person had weight one: its free rows' unweighted +expected deaths equal their unweighted baseline total, and they keep their +relative risks. A factor of zero would instead make those persons +immortal. With unit risks the rule reproduces the DI adapter, which leaves +the factor at one when the other rows carry no expected deaths. If the +weighted target is nonzero in such a cell, because fixed probabilities +differ from baseline, no allocation conserves it: the cell is still +flagged infeasible and the residual records the gap. On the model's path +the fixed rows sit at baseline, so the target is zero and the cell is +feasible. The exact baseline path bypasses multiply/divide rounding for unit risks and for a common free-row risk with fixed rows at baseline. Stable sums @@ -168,7 +186,8 @@ before being characterized as a mechanism defect. | DI composition, ratios and baseline origin q | Fixed origin in `model.py`, existing adapter | Composition examples and `test_composition_conserves_cells_preserves_ratios_and_fixes_di_q` | | Intended monotonicity by earnings and education | Positive cell factor, log-additive risks | `test_ordered_earnings_and_education_risks_preserve_probability_order` | | Row/weight invariance and one-risk cells | Stable reductions, exact common-risk path | Permutation/scaling examples and properties; `test_one_free_risk_group_returns_baseline_probability` | -| Differential adapter arithmetic; intended clipping difference | Same uncapped algebra, iterative cap | Private-adapter differential examples and properties; named intended capping divergence | +| Weightless free support: unit-weight factor, never factor zero | Unit-weight allocation in both arithmetic paths of `normalizer.py` | `test_weightless_support_counts_persons_equally`, `test_weightless_support_uses_unit_weights_on_the_decimal_path`, `test_weightless_free_rows_in_infeasible_cells_keep_unit_weights`, `test_weightless_cells_count_persons_equally`, `test_model_gives_weightless_cells_their_unit_weight_probabilities`; the normalizer properties draw zero weights | +| Differential adapter arithmetic; intended clipping difference | Same uncapped algebra, iterative cap | Private-adapter differential examples and properties, including `test_unit_risk_matches_adapter_on_every_unclipped_cell` over zero weights, infeasible cells and arbitrary fixed probabilities; named intended capping divergence | | Time path and alternative path isolation | Path protocol and hold-last default | `test_pluggable_path_changes_only_named_year`, generated hold-last/refusal property | | Determinism and no RNG consumption | Pure `evaluate` | Repeated-call examples, generator state and patched-generator property | | Refusals | Specification, frame and interface validation | Nonfinite/nonpositive risks, missing columns, unknown/unclassified values, age mismatch and `probabilities` attribute examples and properties | diff --git a/src/populace_dynamics/mortality_ses/normalizer.py b/src/populace_dynamics/mortality_ses/normalizer.py index e009e8ad..59e45fd9 100644 --- a/src/populace_dynamics/mortality_ses/normalizer.py +++ b/src/populace_dynamics/mortality_ses/normalizer.py @@ -86,8 +86,15 @@ def total(values): ) support = [i for i in free if base[i] > zero] capacity = total(w[i] for i in support) + allocation = w + unmet = False + if capacity == zero and support: + unmet = target != zero + allocation = [one] * len(w) + target = total(base[i] for i in support) + capacity = total(allocation[i] for i in support) capped = set() - infeasible = False + infeasible = unmet factor = one if fixed_at_baseline and ( not free or all(rr[i] == rr[free[0]] for i in free) @@ -99,7 +106,7 @@ def total(values): elif target == zero: factor = zero elif target >= capacity: - infeasible = target > capacity + infeasible = target > capacity or unmet for i in support: result[i] = one capped = set(support) @@ -108,16 +115,18 @@ def total(values): ) else: raw = [qi * ri for qi, ri in zip(base, rr, strict=True)] - active = {i for i in support if w[i] > zero} + active = {i for i in support if allocation[i] > zero} remaining = target while active: - factor = remaining / total(w[i] * raw[i] for i in active) + factor = remaining / total( + allocation[i] * raw[i] for i in active + ) next_capped = {i for i in active if raw[i] * factor > one} if not next_capped: break capped |= next_capped active -= next_capped - remaining = target - total(w[i] for i in capped) + remaining = target - total(allocation[i] for i in capped) for i in free: result[i] = min(one, raw[i] * factor) capped = {i for i in free if raw[i] * factor > one} @@ -161,8 +170,18 @@ def _normalize_cell(base, rr, w, is_fixed, prescribed): ) support = free & (base > 0) capacity = _stable_sum(w[support]) + # Weightless support leaves the weighted target silent on the factor; + # set it by counting persons equally rather than letting zero weight + # decide a probability. A nonzero weighted target stays infeasible. + allocation = w + unmet = False + if capacity == 0 and np.any(support): + unmet = bool(target != 0) + allocation = np.ones(len(base), dtype=np.longdouble) + target = _stable_sum(base[support]) + capacity = _stable_sum(allocation[support]) capped = np.zeros(len(base), dtype=bool) - infeasible = False + infeasible = unmet factor = np.longdouble(1) if fixed_at_baseline and ( not np.any(free) or np.all(rr[free] == rr[free][0]) @@ -178,7 +197,9 @@ def _normalize_cell(base, rr, w, is_fixed, prescribed): elif target >= capacity: # Baseline fixed rows are feasible by construction. Independent # floating reductions can differ by a rounding bit at capacity. - infeasible = bool(target > capacity) and not fixed_at_baseline + infeasible = ( + bool(target > capacity) and not fixed_at_baseline + ) or unmet result[support] = 1 capped[support] = True factor = ( @@ -188,17 +209,17 @@ def _normalize_cell(base, rr, w, is_fixed, prescribed): ) else: raw = base * rr - active = support & (w > 0) + active = support & (allocation > 0) remaining = target while np.any(active): - denominator = _stable_sum(w[active] * raw[active]) + denominator = _stable_sum(allocation[active] * raw[active]) factor = remaining / denominator next_capped = active & (raw * factor > 1) if not np.any(next_capped): break capped |= next_capped active &= ~next_capped - remaining = target - _stable_sum(w[capped]) + remaining = target - _stable_sum(allocation[capped]) result[free] = np.minimum(1, raw[free] * factor) capped[free] = raw[free] * factor > 1 @@ -244,11 +265,17 @@ def normalize_probabilities( factor is recomputed after each capped group until the remaining target can be allocated without further capping. Zero-weight rows receive the same factor as other rows but do not affect its value. + When every free row with positive baseline probability has zero + weight, no factor changes the cell's weighted deaths. The cell then + counts persons equally to set the factor, so those rows keep their + relative risks rather than a probability that zero weight decided. + It is still infeasible when its weighted target is nonzero. Fixed probabilities above the entire baseline target are infeasible; free probabilities are then zero. A target above the free rows' positive-probability support is also infeasible; those rows saturate - at one. The signed residual is baseline minus applied deaths. On + at one. Both rules apply only where that support carries weight. + The signed residual is baseline minus applied deaths. On feasible cells without capping it is explicitly zero: arithmetic roundoff is covered by the accounting tolerance, rather than being described as a clipping loss. diff --git a/tests/mortality_ses/test_normalizer.py b/tests/mortality_ses/test_normalizer.py index daaee1ce..dee860b3 100644 --- a/tests/mortality_ses/test_normalizer.py +++ b/tests/mortality_ses/test_normalizer.py @@ -132,6 +132,34 @@ def test_zero_baselines_and_zero_weights_remain_bounded(): assert zero_log.iloc[0].clipping_residual == 0 +def test_weightless_support_counts_persons_equally(): + # With no weight on the free rows that can die, conservation says + # nothing about the factor; zero would make these persons immortal. + unit_q, unit_log = normalize([0.2, 0.6], [1, 2], weights=[1, 1]) + zero_q, zero_log = normalize([0.2, 0.6], [1, 2], weights=[0, 0]) + np.testing.assert_array_equal(zero_q, unit_q) + np.testing.assert_allclose(zero_q, [0.8 / 7, 4.8 / 7]) + assert zero_log.iloc[0].factor == unit_log.iloc[0].factor + assert zero_log.iloc[0].applied_expected_deaths == 0 + assert not zero_log.iloc[0].infeasible + # Positive-weight rows with q_base = 0 leave the weightless row free. + mixed_q, mixed_log = normalize([0, 0.2], [1, 2], weights=[1, 0]) + np.testing.assert_array_equal(mixed_q, [0, 0.2]) + assert mixed_log.iloc[0].factor == pytest.approx(0.5) + # Positive-weight support still sets the factor alone. + q, _ = normalize([0.2, 0.2], [1, 2], weights=[1, 0]) + np.testing.assert_allclose(q, [0.2, 0.4]) + + +def test_weightless_support_uses_unit_weights_on_the_decimal_path(): + risks = [1e-90, 2e-90] + unit_q, _ = normalize([0.2, 0.2], risks, weights=[1, 1]) + zero_q, zero_log = normalize([0.2, 0.2], risks, weights=[0, 0]) + np.testing.assert_array_equal(zero_q, unit_q) + np.testing.assert_allclose(zero_q, [0.4 / 3, 0.8 / 3]) + assert not zero_log.iloc[0].infeasible + + def test_cells_are_separate_and_sorted_and_inputs_are_unchanged(): baseline = np.array([0.3, 0.1, 0.3, 0.1]) risks = np.array([1.0, 1.0, 2.0, 2.0]) diff --git a/tests/mortality_ses/test_properties.py b/tests/mortality_ses/test_properties.py index d3d93ea2..88c25baf 100644 --- a/tests/mortality_ses/test_properties.py +++ b/tests/mortality_ses/test_properties.py @@ -53,6 +53,7 @@ SEXES = ("female", "male") POSITIVE = st.floats(0.05, 20, allow_nan=False, allow_infinity=False) PROBABILITY = st.floats(0, 1, allow_nan=False, allow_infinity=False) +WEIGHT = st.one_of(st.just(0.0), POSITIVE) def invented_dimensions(): @@ -182,9 +183,7 @@ def invented_normalizer_inputs(draw): draw(st.lists(PROBABILITY, min_size=size, max_size=size)) ) risk = np.asarray(draw(st.lists(POSITIVE, min_size=size, max_size=size))) - weights = np.asarray( - draw(st.lists(POSITIVE, min_size=size, max_size=size)) - ) + weights = np.asarray(draw(st.lists(WEIGHT, min_size=size, max_size=size))) fixed = np.asarray( draw(st.lists(st.booleans(), min_size=size, max_size=size)) ) @@ -242,6 +241,23 @@ def test_permutation_and_common_weight_scaling_preserve_probabilities( np.testing.assert_allclose(scaled, q, rtol=RTOL, atol=ATOL) +@PROPERTY +@given(invented_normalizer_inputs()) +def test_weightless_cells_count_persons_equally(inputs): + """Zero weight leaves the factor to unit weights, never to zero.""" + base, risk, fixed, cells, _ = inputs + zero, zero_log = normalize_probabilities( + base, risk, fixed, base, cells, np.zeros(len(base)) + ) + unit, unit_log = normalize_probabilities( + base, risk, fixed, base, cells, np.ones(len(base)) + ) + np.testing.assert_array_equal(zero, unit) + np.testing.assert_array_equal(zero_log["factor"], unit_log["factor"]) + assert not zero_log["infeasible"].any() + assert (zero_log["clipping_residual"] == 0).all() + + @PROPERTY @given(PROBABILITY, POSITIVE, st.integers(1, 15)) def test_one_free_risk_group_returns_baseline_probability(q_base, risk, size): @@ -277,7 +293,7 @@ def test_fixed_deaths_above_target_are_flagged(base_q, excess, weight): @PROPERTY -@given(st.floats(0.02, 0.15), st.floats(0.5, 1.75), POSITIVE) +@given(st.floats(0.02, 0.15), st.floats(0.5, 1.75), WEIGHT) def test_unit_risk_reproduces_adapter_arithmetic_without_clipping( baseline_q, multiplier, weight_scale ): @@ -299,6 +315,69 @@ def test_unit_risk_reproduces_adapter_arithmetic_without_clipping( np.testing.assert_allclose(q, old, rtol=RTOL, atol=ATOL) +@st.composite +def invented_adapter_cells(draw): + size = draw(st.integers(1, 10)) + rows = st.lists(PROBABILITY, min_size=size, max_size=size) + base, prescribed = np.asarray(draw(rows)), np.asarray(draw(rows)) + fixed = np.asarray( + draw(st.lists(st.booleans(), min_size=size, max_size=size)) + ) + prescribed[~fixed] = base[~fixed] + weights = np.asarray(draw(st.lists(WEIGHT, min_size=size, max_size=size))) + groups = np.asarray( + draw(st.lists(st.integers(0, 3), min_size=size, max_size=size)) + ) + return base, fixed, prescribed, weights, groups % 2, 50 + groups // 2 + + +@PROPERTY +@given(invented_adapter_cells()) +def test_unit_risk_matches_adapter_on_every_unclipped_cell(inputs): + """Includes zero weights, infeasible cells and weightless free rows.""" + base, fixed, prescribed, weights, sex, age = inputs + cells = [(2021, SEXES[s], int(a)) for s, a in zip(sex, age, strict=True)] + q, log = normalize_probabilities( + base, np.ones(len(base)), fixed, prescribed, cells, weights + ) + old, _ = _net_of_di_origin( + prescribed, base, fixed, sex, age, weights, SINGLE_YEAR_AGE_BANDS + ) + for row in log.itertuples(index=False): + if row.clipped_persons: + continue # Intended: iterative capping versus one-pass clipping. + mask = np.asarray( + [cell == (row.year, row.sex, row.age) for cell in cells] + ) + np.testing.assert_allclose(q[mask], old[mask], rtol=RTOL, atol=ATOL) + + +def test_weightless_free_rows_in_infeasible_cells_keep_unit_weights(): + """Fixed rows carry the weight; no free factor can conserve the cell.""" + cells = [(2021, "female", 50)] * 3 + fixed = np.asarray([True, False, False]) + for prescribed_q, sign in ((0.6, -1), (0.1, 1)): + base = np.asarray([0.3, 0.2, 0.4]) + prescribed = np.asarray([prescribed_q, 0.2, 0.4]) + weights = np.asarray([2.0, 0.0, 0.0]) + q, log = normalize_probabilities( + base, + np.asarray([1.0, 1.0, 3.0]), + fixed, + prescribed, + cells, + weights, + ) + unit, _ = normalize_probabilities( + base, np.asarray([1.0, 1.0, 3.0]), fixed, base, cells, np.ones(3) + ) + assert log.loc[0, "infeasible"] + assert np.sign(log.loc[0, "clipping_residual"]) == sign + assert q[0] == prescribed_q + np.testing.assert_array_equal(q[1:], unit[1:]) + assert_conserved(q, base, cells, weights, log) + + def test_iterative_capping_intentionally_differs_from_adapter_one_pass(): """Intended difference: reallocate the one-pass clipping shortfall.""" base = np.asarray([0.8, 0.8, 0.1]) @@ -568,6 +647,19 @@ def test_model_permutation_scaling_and_no_randomness(scale, seed): assert person_rng.bit_generator.state == person_before +def test_model_gives_weightless_cells_their_unit_weight_probabilities(): + frame = invented_frame(12) + context = invented_context(frame) + model = DifferentialMortality( + invented_baseline(0.2), invented_spec(), weight_column="weight" + ) + unit = model(frame, context) + weightless = model(frame.assign(weight=0.0), context) + np.testing.assert_array_equal(weightless, unit) + assert (weightless > 0).all() + assert len(np.unique(weightless)) > 1 + + class InventedWideningPath: """A separately registered invented path changing only named years.""" From 6801a686df29eba5b3d3602c59a549b066c43583 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 14:17:22 -0400 Subject: [PATCH 05/17] SM1: tier counts for the six new unit tests tests/mortality_ses collects 96 unit tests, up from 90, and no others. unit 18,178 -> 18,184; total 24,709 -> 24,715. The full-collection tiers_match_policy_manifest check passes with the empty PSID root. Co-Authored-By: Claude Opus 5.5 --- tests/README-tiers.md | 4 ++-- tests/tier_counts.json | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/README-tiers.md b/tests/README-tiers.md index 084d7202..a069e786 100644 --- a/tests/README-tiers.md +++ b/tests/README-tiers.md @@ -39,9 +39,9 @@ pytest --collect-only -q -m oracle_policyengine | tail -1 | Tier | Tests at HEAD | |---|---:| -| `unit` | 18,178 | +| `unit` | 18,184 | | `artifact` | 4,253 | | `integration_psid` | 1,450 | | `reproduction_legacy` | 520 | | `oracle_policyengine` | 308 | -| **Total** | **24,709** | +| **Total** | **24,715** | diff --git a/tests/tier_counts.json b/tests/tier_counts.json index 5f704c67..a2679da6 100644 --- a/tests/tier_counts.json +++ b/tests/tier_counts.json @@ -1,7 +1,7 @@ { "schema_version": 1, "counts": { - "unit": 18178, + "unit": 18184, "artifact": 4253, "integration_psid": 1450, "reproduction_legacy": 520, From b99300cc644dc9e8583c9f00f1a7fbcbb80feed8 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 16:54:25 -0400 Subject: [PATCH 06/17] SM1: settle close float calls exactly and log the normalization basis Review round 1b2 found two float-path faults. Roundoff in the capping loop could zero a capped cell's probabilities and hide the loss, and a rounded tie between target and capacity saturated rows it should not. The float path now hands the Decimal path every call it cannot make reliably: a target within 1e-3 of capacity or of zero, a capping round that leaves less than 1e-3 of the target, and a feasible result that misses conservation by more than 1e-12. Capped rows stay at one. Infeasibility is now decided before the unit-weight rule. A negative weighted target zeroes the free rows and a positive one saturates the support, whatever their weights; these are the limits as the support's weight vanishes. Unit weights apply only where the weighted target is zero. A target that exactly fills the weighted support takes the smallest capping factor, so zero-weight rows are not pushed to one. The cell log gains normalization_basis and persons_target. DifferentialMortality takes an optional cell_log mapping, so a run under the DI adapter keeps that log. Grouping checks missing values once per distinct cell, and the stable sum converts its terms in bulk. The tests add an exact rational oracle, a vanishing-weight limit property, a regression for each finding, and a composed DI run with unequal weights and clipping in both steps. Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 124 ++++-- src/populace_dynamics/mortality_ses/model.py | 15 +- .../mortality_ses/normalizer.py | 174 +++++--- tests/mortality_ses/test_normalizer.py | 116 ++++- tests/mortality_ses/test_properties.py | 401 +++++++++++++++++- 5 files changed, 724 insertions(+), 106 deletions(-) diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md index 34b382dd..9bb154a3 100644 --- a/docs/design/mortality_ses_sm1.md +++ b/docs/design/mortality_ses_sm1.md @@ -40,13 +40,15 @@ All public objects are exported from `populace_dynamics.mortality_ses`. returns `(probabilities, cell_log)`. Inputs are matching row vectors; cell IDs are `(year, sex, age)` triples. Weights are finite and nonnegative; zero-weight rows use the same factor but cannot determine - it, and a cell whose positive-q free rows all weigh zero sets its - factor with unit weights (see "Normalization and tolerance"). Inputs - are never mutated. + it. A feasible cell whose positive-q free rows all weigh zero sets its + factor with unit weights and logs that basis (see "Normalization and + tolerance"). Inputs are never mutated. - `DifferentialMortality(baseline, specification, *, weight_column, - di_configuration=None)` exposes single-year `bands` and + di_configuration=None, cell_log=None)` exposes single-year `bands` and `__call__(frame, context) -> ndarray`. `evaluate(frame, context)` also - returns its cell log. `weight_column=None` explicitly means equal + returns its cell log. The DI adapter keeps only probabilities, so a + `cell_log` mapping, when given, receives each call's cell log keyed by + projection year. `weight_column=None` explicitly means equal person weights. The baseline must be callable with context, have single-year bands, and have no `probabilities` attribute. Neither does the new model: the DI adapter would otherwise call it without the year. @@ -100,42 +102,85 @@ rows. Fixed deaths above the target make the cell infeasible and give all free rows zero. Arbitrary fixed probabilities can also leave a target above free positive-q support; the normalizer saturates that support and flags infeasibility. Neither case can arise with fixed rows at baseline q. -Both rules apply when that support carries weight; the weightless case -follows. +Both rules hold whatever the free rows' weights. + +When the target exactly fills the weighted support, any factor that caps +every weighted row conserves the cell. The normalizer takes the smallest, +which is the limit of the capping loop's factor as the target rises to +capacity, so zero-weight rows are not pushed to one. With fixed rows at +baseline and one free risk, that factor returns the baseline q. Weights decide the factor only when they can. If every free row with -positive baseline q has zero weight, no factor changes the cell's weighted -deaths, so the weighted target says nothing about it. The cell then sets -the factor as if each person had weight one: its free rows' unweighted -expected deaths equal their unweighted baseline total, and they keep their -relative risks. A factor of zero would instead make those persons -immortal. With unit risks the rule reproduces the DI adapter, which leaves -the factor at one when the other rows carry no expected deaths. If the -weighted target is nonzero in such a cell, because fixed probabilities -differ from baseline, no allocation conserves it: the cell is still -flagged infeasible and the residual records the gap. On the model's path -the fixed rows sit at baseline, so the target is zero and the cell is -feasible. +positive baseline q has zero weight and the weighted target is zero, every +factor conserves the cell's weighted deaths, so the weights say nothing +about it. The cell then sets the factor as if each person had weight one: +its free rows' unweighted expected deaths equal their unweighted baseline +total, and they keep their relative risks. A factor of zero would instead +make those persons immortal. The log marks such a cell with +`normalization_basis = "persons"` and reports that unweighted total as +`persons_target`; every other cell has basis `"weights"` and a missing +persons target. On the model's path the fixed rows sit at baseline, so a +weightless cell's target is zero and the rule applies. With fixed rows at +baseline the cell gets the probabilities it would get if every weight +were one. + +Unit weights apply only where the weighted target is zero. When it is +not, because fixed probabilities differ from baseline, the cell is +infeasible and follows the rules above. Those rules are the unique limit +as the free support's weight shrinks to zero: a negative target gives +zero at every small weight, and a positive one saturates the support. A +zero target has no unique limit, since the factor then depends on how the +vanishing weights compare. The unit-weight rule chooses equal weights. + +Zero weights can reach a projection frame: +- persons take their household's weight (`data/deployment_frame.py`), so + a zero household weight reaches every member; +- children take their parent's weight (`engine/steps.py`); +- the projection loop passes the whole frame to the mortality module + (`engine/loop.py`); +- population accounting treats zero-weight rows as ordinary persons and + never drops them (`engine/accounting.py`). + +A weighted gate such as SM3 cannot see these rows' probabilities, so +passing it does not validate the rule; the `normalization_basis` column +shows where it applied. The exact baseline path bypasses multiply/divide rounding for unit risks and for a common free-row risk with fixed rows at baseline. Stable sums -reduce order sensitivity; Decimal arithmetic handles extreme magnitudes -that would underflow or overflow ordinary products. Probabilities returned -are float64; the log accounts for those returned values. +reduce order sensitivity. The Decimal path (precision 2,400 digits) +handles extreme magnitudes that would underflow or overflow ordinary +products. The float path also hands the Decimal path every decision it +cannot make reliably: +- a target within 1e-3 of the free capacity, relative to that capacity; +- a target formed by subtracting fixed deaths that is within 1e-3 of zero, + relative to the larger total; +- a capping round that leaves less than 1e-3 of the target; +- a feasible result that misses conservation by more than 1e-12, + relative to the cell's expected deaths. + +Float reductions err by a few units of 2^-52, so above that margin the +error moves the factor by less than the 1e-12 accounting tolerance. +Probabilities returned are float64; the log accounts for those returned +values. -The cell log contains year, sex, age, persons, fixed persons, total weight, -fixed weight, baseline expected deaths, applied expected deaths, factor, -clipped persons, signed clipping residual and infeasibility, plus the -unclassified counts. `clipping_residual = baseline - applied`. In a -feasible cell without capping it is explicitly zero; ordinary rounding -is handled by the tolerance, rather than labeled a clipping loss. +The cell log contains year, sex, age, persons, fixed persons, total +weight, fixed weight, baseline expected deaths, applied expected deaths, +factor, clipped persons, signed clipping residual, infeasibility, +normalization basis and persons target, plus the unclassified counts. +`clipping_residual = baseline - applied`. It is nonzero only on +infeasible cells: iterative capping conserves every feasible cell, and the +float path verifies that rather than assuming it. Conservation checks use `rtol=1e-12, atol=1e-12` in expected-death units. This allows float64 product, summation and multiply/divide roundoff and the existing adapter's separate reductions; it is an accounting tolerance, not an empirical admission threshold. Identity checks use exact array equality. Permutation checks are exact where supported by stable sums; -weight scaling and adapter comparisons use the stated tolerance. +weight scaling and adapter comparisons use the stated tolerance. An +independent exact oracle in the tests, in rational arithmetic, finds the +capped rows by searching over how many rows cap rather than by iterating. +The normalizer must match it to `rtol=1e-9, atol=1e-15` on cells across +float magnitudes, near-one q and arbitrary fixed probabilities. Intended monotonicity applies among free rows with the same baseline q in a cell, holding other dimensions constant. This is the setting supplied @@ -148,8 +193,14 @@ fall short of the target; its residual is computed from `expected_deaths_population_model - expected_deaths_applied` in its existing cell log. The SM1 log describes normalization before that step. Tests label the normalizer's iterative redistribution versus the -adapter's one-pass clipping as an intended difference. No shared-module -change is required. +adapter's one-pass clipping as an intended difference. A second intended +difference: in an infeasible cell whose free support weighs nothing, the +adapter leaves the factor at one, while the rules above zero or saturate +those rows. No shared-module change is required. + +At the scale of a projection frame (1,000,000 invented rows in 242 +cells), `evaluate` takes about 1.3 seconds on the build host, and the +normalizer alone takes 0.4–0.6 seconds. ## Design and build traceability @@ -178,16 +229,19 @@ before being characterized as a mechanism defect. | Invariant | Implementation | Example and property coverage | |---|---|---| -| Cell conservation and zero uncapped residual | Iterative cap in `normalizer.py` | `test_invented_uncapped_cell_conserves_and_logs_counts`, `test_cell_conservation_bounds_and_baseline_fixed_feasibility` | +| Cell conservation; residual nonzero only when infeasible | Iterative cap and the float path's conservation check in `normalizer.py` | `test_invented_uncapped_cell_conserves_and_logs_counts`, `test_cell_conservation_bounds_and_baseline_fixed_feasibility`, `test_normalizer_matches_an_exact_oracle_across_magnitudes`, `test_cascading_caps_match_the_oracle`; `assert_conserved` requires conservation outright on feasible cells | | Feasibility with baseline fixed q; arbitrary excess flagged | Fixed target calculation | `test_fixed_deaths_above_invented_target_are_infeasible`, `test_fixed_deaths_above_target_are_flagged`, baseline-fixed property | | Unit-risk bit identity | Exact common-risk baseline path | `test_equal_risks_preserve_baseline_bits_with_fixed_baseline_rows`, `test_identity_is_bit_for_bit_and_normalizer_is_deterministic` | | Same deaths with/without registry, ordinary and DI steps | `step.py` and context-only wrapper | `test_identity_probabilities_and_step_deaths`, `test_di_adapter_identity`, corresponding generated death-identity properties | | Bounds | Capping and input validation | Conservation/bounds property and extreme-risk examples | -| DI composition, ratios and baseline origin q | Fixed origin in `model.py`, existing adapter | Composition examples and `test_composition_conserves_cells_preserves_ratios_and_fixes_di_q` | +| DI composition, ratios and baseline origin q | Fixed origin in `model.py`, existing adapter | Composition examples, `test_composition_conserves_cells_preserves_ratios_and_fixes_di_q`, and, with unequal weights, clipping in both steps and the retained model log, `test_composition_with_unequal_weights_keeps_every_accounting_link` and `test_composition_accounts_for_clipping_in_both_steps` | | Intended monotonicity by earnings and education | Positive cell factor, log-additive risks | `test_ordered_earnings_and_education_risks_preserve_probability_order` | | Row/weight invariance and one-risk cells | Stable reductions, exact common-risk path | Permutation/scaling examples and properties; `test_one_free_risk_group_returns_baseline_probability` | -| Weightless free support: unit-weight factor, never factor zero | Unit-weight allocation in both arithmetic paths of `normalizer.py` | `test_weightless_support_counts_persons_equally`, `test_weightless_support_uses_unit_weights_on_the_decimal_path`, `test_weightless_free_rows_in_infeasible_cells_keep_unit_weights`, `test_weightless_cells_count_persons_equally`, `test_model_gives_weightless_cells_their_unit_weight_probabilities`; the normalizer properties draw zero weights | -| Differential adapter arithmetic; intended clipping difference | Same uncapped algebra, iterative cap | Private-adapter differential examples and properties, including `test_unit_risk_matches_adapter_on_every_unclipped_cell` over zero weights, infeasible cells and arbitrary fixed probabilities; named intended capping divergence | +| Weightless free support with a zero target: unit-weight factor, logged as basis `persons` | Unit-weight allocation in both arithmetic paths of `normalizer.py`; `normalization_basis` and `persons_target` log columns; the model's optional `cell_log` | `test_weightless_support_counts_persons_equally`, `test_weightless_support_uses_unit_weights_on_the_decimal_path`, `test_weightless_cells_count_persons_equally`, `test_weightless_cells_take_the_vanishing_weight_limit`, `test_model_gives_weightless_cells_their_unit_weight_probabilities`; the normalizer properties draw zero weights | +| Infeasible weightless cells follow the weighted rules (the vanishing-weight limit) | Infeasibility decided before the unit-weight rule | `test_fixed_excess_zeroes_weightless_free_rows_on_both_paths`, `test_fixed_shortfall_saturates_weightless_support_on_both_paths`, `test_weightless_infeasible_cells_follow_the_weighted_rules`, `test_weightless_cells_take_the_vanishing_weight_limit` | +| Float decisions too close to call go to the Decimal path | `_ambiguous`, the remaining-target check and the conservation check in `normalizer.py` | `test_roundoff_cannot_erase_a_capped_cells_deaths`, `test_rounded_target_and_capacity_tie_is_decided_exactly`, `test_cancelling_target_is_decided_exactly`; the oracle property's examples | +| A target that exactly fills the weighted support takes the smallest capping factor | The Decimal path's tie branch | The oracle property's two tie examples | +| Differential adapter arithmetic; intended clipping and weightless-infeasible differences | Same uncapped algebra, iterative cap | Private-adapter differential examples and properties, including `test_unit_risk_matches_adapter_on_every_unclipped_cell` over zero weights, infeasible cells and arbitrary fixed probabilities; named intended divergences | | Time path and alternative path isolation | Path protocol and hold-last default | `test_pluggable_path_changes_only_named_year`, generated hold-last/refusal property | | Determinism and no RNG consumption | Pure `evaluate` | Repeated-call examples, generator state and patched-generator property | | Refusals | Specification, frame and interface validation | Nonfinite/nonpositive risks, missing columns, unknown/unclassified values, age mismatch and `probabilities` attribute examples and properties | diff --git a/src/populace_dynamics/mortality_ses/model.py b/src/populace_dynamics/mortality_ses/model.py index ee41b389..7fa81613 100644 --- a/src/populace_dynamics/mortality_ses/model.py +++ b/src/populace_dynamics/mortality_ses/model.py @@ -2,6 +2,7 @@ from __future__ import annotations +from collections.abc import MutableMapping from dataclasses import dataclass, field from typing import ClassVar, Protocol @@ -77,6 +78,8 @@ class DifferentialMortality: This model deliberately exposes no ``probabilities`` attribute: the DI adapter dispatches that attribute without a context or a projection year. ``evaluate`` additionally returns the pure normalizer's cell accounting. + The adapter keeps only probabilities, so ``cell_log``, when given, + receives each ``__call__``'s cell accounting keyed by projection year. """ baseline: YearAwareMortality @@ -85,6 +88,9 @@ class DifferentialMortality: di_configuration: DIEntitlementRates | DIEntitlementSpec | None = field( default=None, kw_only=True ) + cell_log: MutableMapping[int, pd.DataFrame] | None = field( + default=None, kw_only=True, compare=False, repr=False + ) bands: ClassVar[tuple[tuple[int, int], ...]] = SINGLE_YEAR_AGE_BANDS def __post_init__(self) -> None: @@ -110,6 +116,10 @@ def __post_init__(self) -> None: self.di_configuration, (DIEntitlementRates, DIEntitlementSpec) ): raise TypeError("di_configuration must be DI rates or their spec") + if self.cell_log is not None and not isinstance( + self.cell_log, MutableMapping + ): + raise TypeError("cell_log must be a mutable mapping or None") def _fixed_rows(self, frame: pd.DataFrame) -> np.ndarray: columns = {"di_entitled", "di_conversion_year"} @@ -214,4 +224,7 @@ def evaluate( def __call__( self, frame: pd.DataFrame, context: PeriodContext ) -> np.ndarray: - return self.evaluate(frame, context)[0] + probabilities, cells = self.evaluate(frame, context) + if self.cell_log is not None: + self.cell_log[int(context.year)] = cells + return probabilities diff --git a/src/populace_dynamics/mortality_ses/normalizer.py b/src/populace_dynamics/mortality_ses/normalizer.py index 59e45fd9..8c744b92 100644 --- a/src/populace_dynamics/mortality_ses/normalizer.py +++ b/src/populace_dynamics/mortality_ses/normalizer.py @@ -28,10 +28,22 @@ "clipped_persons", "clipping_residual", "infeasible", + "normalization_basis", + "persons_target", "unclassified_refuse", "unclassified_reference", ) +# The float path leaves a decision to the Decimal path when it compares +# totals this close, relative to their size, or when capping leaves a +# remaining target this small. Float reductions here err by a few units +# of 2**-52; above this margin that error moves the factor by less than +# the 1e-12 accounting tolerance. +_FLOAT_MARGIN = 1e-3 +# A feasible float cell that misses conservation by more than this, +# relative to its expected deaths, is recomputed on the Decimal path. +_CONSERVATION_RTOL = 1e-12 + def _stable_sum(values: np.ndarray) -> np.longdouble: """Sum nonnegative values without a row-order-dependent reduction. @@ -45,16 +57,24 @@ def _stable_sum(values: np.ndarray) -> np.longdouble: maximum = np.max(values) if maximum == 0: return np.longdouble(0) - return maximum * np.longdouble(fsum(float(x / maximum) for x in values)) + scaled = (values / maximum).astype(np.float64).tolist() + return maximum * np.longdouble(fsum(scaled)) + + +def _ambiguous(difference, scale) -> bool: + """Whether a float difference is too close to zero to decide on.""" + return bool(scale > 0 and abs(difference) <= _FLOAT_MARGIN * scale) def _decimal_cell(base, rr, w, is_fixed, prescribed): - """Handle extreme finite float64 inputs without platform assumptions. + """Handle extreme or ambiguous cells without platform assumptions. On some supported platforms ``longdouble`` has float64's exponent range. Decimal's exponent range retains positive q*RR support even when that product would underflow, and avoids overflowing products. - This slower path is restricted to cells containing extreme magnitudes. + Its precision also settles the comparisons and cancellations that the + float path refers here. This slower path is used only for those + cells. """ with localcontext() as decimal_context: # Enough to retain exact products of three finite binary64 inputs, @@ -85,16 +105,14 @@ def total(values): else baseline_deaths - fixed_deaths ) support = [i for i in free if base[i] > zero] - capacity = total(w[i] for i in support) allocation = w - unmet = False - if capacity == zero and support: - unmet = target != zero + persons_target = None + if target == zero and support and total(w[i] for i in support) == zero: allocation = [one] * len(w) - target = total(base[i] for i in support) - capacity = total(allocation[i] for i in support) + target = persons_target = total(base[i] for i in support) + capacity = total(allocation[i] for i in support) capped = set() - infeasible = unmet + infeasible = False factor = one if fixed_at_baseline and ( not free or all(rr[i] == rr[free[0]] for i in free) @@ -105,14 +123,25 @@ def total(values): infeasible, factor = True, zero elif target == zero: factor = zero - elif target >= capacity: - infeasible = target > capacity or unmet + elif target > capacity: + infeasible = True for i in support: result[i] = one capped = set(support) factor = max( (one / (base[i] * rr[i]) for i in support), default=one ) + elif target == capacity: + # Any factor that caps the whole weighted support conserves the + # cell. The smallest is the capping loop's limit as the target + # rises to capacity; zero-allocation rows take it too. + weighted = [i for i in support if allocation[i] > zero] + factor = max(one / (base[i] * rr[i]) for i in weighted) + for i in free: + result[i] = min(one, base[i] * rr[i] * factor) + for i in weighted: + result[i] = one + capped = {i for i in free if base[i] * rr[i] * factor > one} else: raw = [qi * ri for qi, ri in zip(base, rr, strict=True)] active = {i for i in support if allocation[i] > zero} @@ -128,14 +157,16 @@ def total(values): active -= next_capped remaining = target - total(allocation[i] for i in capped) for i in free: - result[i] = min(one, raw[i] * factor) - capped = {i for i in free if raw[i] * factor > one} + result[i] = one if i in capped else min(one, raw[i] * factor) + capped |= {i for i in free if raw[i] * factor > one} probabilities = np.asarray([float(value) for value in result]) applied = total( wi * Decimal.from_float(float(qi)) for wi, qi in zip(w, probabilities, strict=True) ) - residual = baseline_deaths - applied if infeasible or capped else zero + # Feasible cells conserve exactly; only the conversion of the + # returned probabilities to float64 rounds. + residual = baseline_deaths - applied if infeasible else zero return probabilities, { "weight": float(total(w)), "fixed_weight": float(total(w[i] for i in fixed)), @@ -145,43 +176,53 @@ def total(values): "clipped_persons": len(capped), "clipping_residual": float(residual), "infeasible": infeasible, + "normalization_basis": ( + "weights" if persons_target is None else "persons" + ), + "persons_target": ( + float("nan") + if persons_target is None + else float(persons_target) + ), } def _normalize_cell(base, rr, w, is_fixed, prescribed): + def exact(): + return _decimal_cell(base, rr, w, is_fixed, prescribed) + # Within these magnitudes, coefficients, weighted products, and capped # scaling products remain representable even on float64-only platforms. if any( np.any((values > 0) & (values < 1e-80)) or np.any(values > 1e80) for values in (base, rr, w, prescribed) ): - return _decimal_cell(base, rr, w, is_fixed, prescribed) + return exact() free = ~is_fixed result = np.zeros(len(base), dtype=np.longdouble) result[is_fixed] = prescribed[is_fixed] baseline_deaths = _stable_sum(w * base) fixed_deaths = _stable_sum(w[is_fixed] * prescribed[is_fixed]) fixed_at_baseline = np.array_equal(prescribed[is_fixed], base[is_fixed]) - # Avoid cancellation of small free totals when fixed rows dominate. - target = ( - _stable_sum(w[free] * base[free]) - if fixed_at_baseline - else baseline_deaths - fixed_deaths - ) + if fixed_at_baseline: + # Avoid cancellation of small free totals when fixed rows dominate. + target = _stable_sum(w[free] * base[free]) + else: + target = baseline_deaths - fixed_deaths + if _ambiguous(target, max(baseline_deaths, fixed_deaths)): + return exact() support = free & (base > 0) - capacity = _stable_sum(w[support]) - # Weightless support leaves the weighted target silent on the factor; - # set it by counting persons equally rather than letting zero weight - # decide a probability. A nonzero weighted target stays infeasible. + # Weightless support leaves a zero weighted target silent on the + # factor: every factor conserves the cell. Count persons equally to + # set it, rather than letting zero weight decide a probability. allocation = w - unmet = False - if capacity == 0 and np.any(support): - unmet = bool(target != 0) + persons_target = None + if target == 0 and np.any(support) and _stable_sum(w[support]) == 0: allocation = np.ones(len(base), dtype=np.longdouble) - target = _stable_sum(base[support]) - capacity = _stable_sum(allocation[support]) + target = persons_target = _stable_sum(base[support]) + capacity = _stable_sum(allocation[support]) capped = np.zeros(len(base), dtype=bool) - infeasible = unmet + infeasible = False factor = np.longdouble(1) if fixed_at_baseline and ( not np.any(free) or np.all(rr[free] == rr[free][0]) @@ -194,12 +235,11 @@ def _normalize_cell(base, rr, w, is_fixed, prescribed): factor = np.longdouble(0) elif target == 0: factor = np.longdouble(0) - elif target >= capacity: - # Baseline fixed rows are feasible by construction. Independent - # floating reductions can differ by a rounding bit at capacity. - infeasible = ( - bool(target > capacity) and not fixed_at_baseline - ) or unmet + elif _ambiguous(target - capacity, capacity): + return exact() + elif target > capacity: + # Ties went to the Decimal path, so this cell is infeasible. + infeasible = True result[support] = 1 capped[support] = True factor = ( @@ -220,13 +260,22 @@ def _normalize_cell(base, rr, w, is_fixed, prescribed): capped |= next_capped active &= ~next_capped remaining = target - _stable_sum(allocation[capped]) + if remaining <= _FLOAT_MARGIN * target: + return exact() result[free] = np.minimum(1, raw[free] * factor) - capped[free] = raw[free] * factor > 1 + result[capped] = 1 + capped |= free & (raw * factor > 1) probabilities = result.astype(float) applied_deaths = _stable_sum(w * probabilities) residual = baseline_deaths - applied_deaths - if not infeasible and not np.any(capped): + if not infeasible: + # Iterative capping conserves every feasible cell, so a material + # float residual is an arithmetic failure, never a clipping loss. + if abs(residual) > _CONSERVATION_RTOL * max( + baseline_deaths, applied_deaths + ): + return exact() residual = np.longdouble(0) return probabilities, { "weight": _stable_sum(w), @@ -237,6 +286,12 @@ def _normalize_cell(base, rr, w, is_fixed, prescribed): "clipped_persons": int(np.count_nonzero(capped)), "clipping_residual": residual, "infeasible": infeasible, + "normalization_basis": ( + "weights" if persons_target is None else "persons" + ), + "persons_target": ( + float("nan") if persons_target is None else float(persons_target) + ), } @@ -265,20 +320,30 @@ def normalize_probabilities( factor is recomputed after each capped group until the remaining target can be allocated without further capping. Zero-weight rows receive the same factor as other rows but do not affect its value. - When every free row with positive baseline probability has zero - weight, no factor changes the cell's weighted deaths. The cell then - counts persons equally to set the factor, so those rows keep their - relative risks rather than a probability that zero weight decided. - It is still infeasible when its weighted target is nonzero. Fixed probabilities above the entire baseline target are infeasible; free probabilities are then zero. A target above the free rows' positive-probability support is also infeasible; those rows saturate - at one. Both rules apply only where that support carries weight. - The signed residual is baseline minus applied deaths. On - feasible cells without capping it is explicitly zero: arithmetic - roundoff is covered by the accounting tolerance, rather than being - described as a clipping loss. + at one. These rules hold whatever the free rows' weights. A target + that exactly fills the weighted support conserves with any factor + that caps it; the cell takes the smallest such factor, so zero-weight + rows keep ``min(1, q_base * relative_risk * factor)``. + + When every free row with positive baseline probability has zero + weight and the weighted target is zero, every factor conserves the + cell, so the weights cannot set it. The cell then counts persons + equally: the factor makes those rows' unweighted expected deaths equal + their unweighted baseline total, which the log reports as + ``persons_target`` with ``normalization_basis == "persons"``. Every + other cell has basis ``"weights"`` and a missing persons target. + With fixed rows at baseline, such a cell gets the probabilities it + would get if every weight were one. + + The signed residual is baseline minus applied deaths. It is nonzero + only on infeasible cells; iterative capping conserves feasible cells. + The float path verifies that within the accounting tolerance, and + recomputes on the Decimal path any cell it cannot decide or conserve + in float arithmetic. All relative risks must be finite and strictly positive, including those on fixed rows. Probabilities lie in [0, 1], weights are finite @@ -310,15 +375,16 @@ def normalize_probabilities( for index, cell in enumerate(ids): if not isinstance(cell, (tuple, list, np.ndarray)) or len(cell) != 3: raise ValueError("cell ids must be (year, sex, age) triples") - key = tuple(cell) - if any(pd.isna(value) for value in key): - raise ValueError("cell ids must not contain missing values") try: - groups.setdefault(key, []).append(index) + groups.setdefault(tuple(cell), []).append(index) except TypeError as error: raise ValueError( "cell ids must contain hashable values" ) from error + # Every row's key is one of the distinct keys, so checking each + # distinct key once checks every row. + if any(pd.isna(value) for key in groups for value in key): + raise ValueError("cell ids must not contain missing values") counts = { "unclassified_refuse": np.zeros(size, dtype=np.int64), diff --git a/tests/mortality_ses/test_normalizer.py b/tests/mortality_ses/test_normalizer.py index dee860b3..dfdf1b23 100644 --- a/tests/mortality_ses/test_normalizer.py +++ b/tests/mortality_ses/test_normalizer.py @@ -1,10 +1,13 @@ """Cell accounting tests using invented probabilities, risks and weights.""" +from unittest.mock import patch + import numpy as np import pandas as pd import pytest from populace_dynamics.engine.di_entitlement import _net_of_di_origin +from populace_dynamics.mortality_ses import normalizer from populace_dynamics.mortality_ses.normalizer import ( CELL_LOG_COLUMNS, normalize_probabilities, @@ -49,6 +52,8 @@ def test_invented_uncapped_cell_conserves_and_logs_counts(): assert cell.clipped_persons == 0 assert cell.clipping_residual == 0 assert not cell.infeasible + assert cell.normalization_basis == "weights" + assert np.isnan(cell.persons_target) assert cell.unclassified_reference == 3 assert cell.unclassified_refuse == 0 @@ -142,22 +147,118 @@ def test_weightless_support_counts_persons_equally(): assert zero_log.iloc[0].factor == unit_log.iloc[0].factor assert zero_log.iloc[0].applied_expected_deaths == 0 assert not zero_log.iloc[0].infeasible + # The log names the basis and the unweighted target it conserved. + assert zero_log.iloc[0].normalization_basis == "persons" + assert zero_log.iloc[0].persons_target == pytest.approx(0.8) + assert unit_log.iloc[0].normalization_basis == "weights" + assert np.isnan(unit_log.iloc[0].persons_target) # Positive-weight rows with q_base = 0 leave the weightless row free. mixed_q, mixed_log = normalize([0, 0.2], [1, 2], weights=[1, 0]) np.testing.assert_array_equal(mixed_q, [0, 0.2]) assert mixed_log.iloc[0].factor == pytest.approx(0.5) + assert mixed_log.iloc[0].normalization_basis == "persons" + assert mixed_log.iloc[0].persons_target == pytest.approx(0.2) # Positive-weight support still sets the factor alone. - q, _ = normalize([0.2, 0.2], [1, 2], weights=[1, 0]) + q, log = normalize([0.2, 0.2], [1, 2], weights=[1, 0]) np.testing.assert_allclose(q, [0.2, 0.4]) + assert log.iloc[0].normalization_basis == "weights" def test_weightless_support_uses_unit_weights_on_the_decimal_path(): risks = [1e-90, 2e-90] - unit_q, _ = normalize([0.2, 0.2], risks, weights=[1, 1]) + unit_q, unit_log = normalize([0.2, 0.2], risks, weights=[1, 1]) zero_q, zero_log = normalize([0.2, 0.2], risks, weights=[0, 0]) np.testing.assert_array_equal(zero_q, unit_q) np.testing.assert_allclose(zero_q, [0.4 / 3, 0.8 / 3]) assert not zero_log.iloc[0].infeasible + assert zero_log.iloc[0].normalization_basis == "persons" + assert zero_log.iloc[0].persons_target == pytest.approx(0.4) + assert unit_log.iloc[0].normalization_basis == "weights" + assert np.isnan(unit_log.iloc[0].persons_target) + + +@pytest.mark.parametrize("risks", [[1, 2], [1, 1e100]]) +def test_fixed_excess_zeroes_weightless_free_rows_on_both_paths(risks): + # Fixed deaths alone exceed the cell's baseline deaths, so free rows + # are zero whatever their weights; unit weights apply only when the + # weighted target is zero. + q, log = normalize( + [0.25, 0.25], + risks, + fixed=[True, False], + prescribed=[0.5, 0.25], + weights=[1, 0], + ) + np.testing.assert_array_equal(q, [0.5, 0]) + assert log.iloc[0].infeasible + assert log.iloc[0].clipping_residual == pytest.approx(-0.25) + assert log.iloc[0].normalization_basis == "weights" + + +@pytest.mark.parametrize("risks", [[1, 1, 3], [1, 1e-90, 3e-90]]) +def test_fixed_shortfall_saturates_weightless_support_on_both_paths(risks): + # The weighted target is positive but the free support carries no + # weight, so no allocation meets it; the support saturates, as it + # does for any target above the support's weight. + q, log = normalize( + [0.3, 0.2, 0.4], + risks, + fixed=[True, False, False], + prescribed=[0.1, 0.2, 0.4], + weights=[2, 0, 0], + ) + np.testing.assert_array_equal(q, [0.1, 1, 1]) + assert log.iloc[0].infeasible + assert log.iloc[0].clipping_residual == pytest.approx(0.4) + assert log.iloc[0].normalization_basis == "weights" + + +def test_roundoff_cannot_erase_a_capped_cells_deaths(): + # Float roundoff caps the second row and leaves no remaining target; + # the cell must still conserve its 1.1 expected deaths. + weights = np.array([1, 0.1, 1.1]) + with patch.object( + normalizer, "_decimal_cell", wraps=normalizer._decimal_cell + ) as exact: + q, log = normalize([0.5] * 3, [20, 1, 1e-25], weights=weights) + exact.assert_called_once() + assert q[0] == 1 + assert q[1] == pytest.approx(1, abs=1e-12) + assert 0 <= q[2] < 1e-15 + assert np.sum(weights * q) == pytest.approx(1.1, rel=1e-12) + assert log.iloc[0].clipped_persons >= 1 + assert log.iloc[0].clipping_residual == 0 + assert not log.iloc[0].infeasible + + +def test_rounded_target_and_capacity_tie_is_decided_exactly(): + # Rounded to float64, the target 1e16 + 0.5 and the capacity 1e16 + 1 + # are equal; exactly, the target is lower and nothing saturates. + with patch.object( + normalizer, "_decimal_cell", wraps=normalizer._decimal_cell + ) as exact: + q, log = normalize([1, 0.5], [1, 1e-20], weights=[1e16, 1]) + exact.assert_called_once() + np.testing.assert_allclose(q, [1, 0.5], rtol=1e-12) + assert log.iloc[0].clipped_persons == 1 + assert not log.iloc[0].infeasible + + +def test_cancelling_target_is_decided_exactly(): + # Fixed deaths equal baseline deaths up to float cancellation. + with patch.object( + normalizer, "_decimal_cell", wraps=normalizer._decimal_cell + ) as exact: + q, log = normalize( + [0.3, 0.1, 0.2], + [1, 1, 2], + fixed=[True, False, False], + prescribed=[0.6, 0.1, 0.2], + ) + exact.assert_called_once() + np.testing.assert_allclose(q, [0.6, 0, 0], atol=1e-15) + assert not log.iloc[0].infeasible + assert log.iloc[0].clipping_residual == 0 def test_cells_are_separate_and_sorted_and_inputs_are_unchanged(): @@ -181,6 +282,17 @@ def test_cells_are_separate_and_sorted_and_inputs_are_unchanged(): np.testing.assert_array_equal(actual, original) +def test_unhashable_and_missing_cell_ids_are_refused_on_every_row(): + base = [0.2, 0.2, 0.2] + for cells, message in ( + ([(2030, "male", 70)] * 2 + [(2030, float("nan"), 70)], "missing"), + ([(2030, "male", 70)] * 2 + [(2030, pd.NA, 70)], "missing"), + ([(2030, "male", 70)] * 2 + [(2030, ["male"], 70)], "hashable"), + ): + with pytest.raises(ValueError, match=message): + normalize(base, [1, 1, 1], cells=cells) + + def test_empty_input_has_declared_log_schema(): q, log = normalize([], []) assert q.shape == (0,) diff --git a/tests/mortality_ses/test_properties.py b/tests/mortality_ses/test_properties.py index 88c25baf..8b2665a1 100644 --- a/tests/mortality_ses/test_properties.py +++ b/tests/mortality_ses/test_properties.py @@ -10,12 +10,13 @@ from copy import deepcopy from dataclasses import replace +from fractions import Fraction from unittest.mock import patch import numpy as np import pandas as pd import pytest -from hypothesis import given, settings +from hypothesis import example, given, settings from hypothesis import strategies as st from populace_dynamics.baselines.base import BaselineYearAwareMortality @@ -47,6 +48,9 @@ RTOL = ATOL = 1e-12 PROPERTY = settings(max_examples=35, deadline=None) +# The exact oracle is cheap on these small cells and is the main check +# of the arithmetic, so it gets more examples. +ORACLE_PROPERTY = settings(max_examples=150, deadline=None) REF = "invented_ref" OTHER = "invented_other" UNCLASSIFIED = "invented_unclassified" @@ -157,6 +161,7 @@ def invented_rates(*, multiplier=1.5, post_conversion="di_origin"): def assert_conserved(q, base, cells, weights, log): + """Feasible cells conserve outright; only infeasible ones log a gap.""" for row in log.itertuples(index=False): mask = np.asarray( [cell == (row.year, row.sex, row.age) for cell in cells] @@ -172,8 +177,233 @@ def assert_conserved(q, base, cells, weights, log): assert row.applied_expected_deaths == pytest.approx( applied, rel=RTOL, abs=ATOL ) - if not row.clipped_persons and not row.infeasible: + if not row.infeasible: assert row.clipping_residual == 0 + assert applied == pytest.approx(target, rel=RTOL, abs=ATOL) + + +def exact_cell(base, risk, fixed, prescribed, weights): + """An independent exact reading of the normalization rule. + + It uses rational arithmetic and finds the capped rows by trying each + count of largest-q*RR rows in turn, rather than iterating the factor + as the normalizer does. Returns probabilities, the infeasible flag + and the normalization basis. + """ + b, r, p, w = ( + [Fraction(float(x)) for x in values] + for values in (base, risk, prescribed, weights) + ) + rows = range(len(b)) + free = [i for i in rows if not fixed[i]] + q = [p[i] if fixed[i] else Fraction(0) for i in rows] + target = sum((w[i] * b[i] for i in rows), Fraction(0)) - sum( + (w[i] * p[i] for i in rows if fixed[i]), Fraction(0) + ) + support = [i for i in free if b[i] > 0] + allocation, basis = w, "weights" + if target == 0 and support and not any(w[i] for i in support): + allocation, basis = [Fraction(1)] * len(b), "persons" + target = sum(b[i] for i in support) + capacity = sum((allocation[i] for i in support), Fraction(0)) + if target <= 0: + return [float(x) for x in q], target < 0, basis + if target > capacity: + for i in support: + q[i] = Fraction(1) + return [float(x) for x in q], True, basis + raw = [b[i] * r[i] for i in rows] + if target == capacity: + # Every weighted row must reach one; zero-allocation rows take the + # smallest factor that does that. + factor = max(1 / raw[i] for i in support if allocation[i] > 0) + for i in free: + q[i] = min(Fraction(1), raw[i] * factor) + return [float(x) for x in q], False, basis + ranked = sorted( + (i for i in support if allocation[i] > 0), + key=lambda i: raw[i], + reverse=True, + ) + for count in range(len(ranked)): + top, rest = ranked[:count], ranked[count:] + factor = (target - sum(allocation[i] for i in top)) / sum( + allocation[i] * raw[i] for i in rest + ) + if all(raw[i] * factor <= 1 for i in rest): + break + else: # pragma: no cover - a target below capacity always stops + raise AssertionError("no capped set satisfies the cell") + for i in free: + q[i] = min(Fraction(1), raw[i] * factor) + return [float(x) for x in q], False, basis + + +def wide(low, high): + """Invented magnitudes spread evenly in log10 between the bounds.""" + return st.floats(low, high).map(lambda exponent: 10.0**exponent) + + +@st.composite +def invented_wide_cells(draw): + """Cells across float magnitudes, near-one q and arbitrary fixed q.""" + size = draw(st.integers(1, 9)) + rows = dict(min_size=size, max_size=size) + base = draw( + st.lists( + st.one_of( + st.sampled_from((0.0, 1.0)), PROBABILITY, st.floats(0.9, 1) + ), + **rows, + ) + ) + risk = draw( + st.lists( + st.one_of(wide(-30, 30), st.sampled_from((1e-100, 1e100))), + **rows, + ) + ) + weights = draw(st.lists(st.one_of(st.just(0.0), wide(-8, 16)), **rows)) + fixed = draw(st.lists(st.booleans(), **rows)) + prescribed = list(base) + if draw(st.booleans()): + prescribed = [ + draw(PROBABILITY) if is_fixed else q + for q, is_fixed in zip(base, fixed, strict=True) + ] + return tuple( + np.asarray(values, dtype=dtype) + for values, dtype in ( + (base, float), + (risk, float), + (fixed, bool), + (prescribed, float), + (weights, float), + ) + ) + + +@ORACLE_PROPERTY +@given(invented_wide_cells()) +@example( # Float roundoff once emptied the remaining target here. + ( + np.array([0.5, 0.5, 0.5]), + np.array([20, 1, 1e-25]), + np.zeros(3, dtype=bool), + np.array([0.5, 0.5, 0.5]), + np.array([1, 0.1, 1.1]), + ) +) +@example( # A rounded target-capacity tie once saturated both rows. + ( + np.array([1, 0.5]), + np.array([1, 1e-20]), + np.zeros(2, dtype=bool), + np.array([1, 0.5]), + np.array([1e16, 1]), + ) +) +@example( # An exactly filled weighted support beside weightless rows. + ( + np.array([1, 2.2691957742722138e-77, 1, 1, 0.05303249588052529]), + np.array( + [1e30, 4.0259657179840314e17, 4.0259657179840314e17, 1e30, 1e30] + ), + np.array([False, True, True, False, False]), + np.array([1, 2.2691957742722138e-77, 1, 1, 0.05303249588052529]), + np.array([0, 0, 131.91017397114445, 1.000000137244776, 0]), + ) +) +@example( # The same tie with unequal risks, on the Decimal path. + ( + np.array([1, 0.5, 0.2]), + np.array([1, 3, 2]), + np.zeros(3, dtype=bool), + np.array([1, 0.5, 0.2]), + np.array([2, 0, 0]), + ) +) +@example( # Fixed excess with weightless free rows. + ( + np.array([0.25, 0.25]), + np.array([1, 2]), + np.array([True, False]), + np.array([0.5, 0.25]), + np.array([1, 0]), + ) +) +def test_normalizer_matches_an_exact_oracle_across_magnitudes(inputs): + base, risk, fixed, prescribed, weights = inputs + cells = [(2021, "female", 50)] * len(base) + q, log = normalize_probabilities( + base, risk, fixed, prescribed, cells, weights + ) + expected, infeasible, basis = exact_cell( + base, risk, fixed, prescribed, weights + ) + np.testing.assert_allclose(q, expected, rtol=1e-9, atol=1e-15) + assert bool(log.loc[0, "infeasible"]) == infeasible + assert log.loc[0, "normalization_basis"] == basis + assert ((0 <= q) & (q <= 1)).all() + if not infeasible: + applied = sum( + Fraction(x) * Fraction(y) for x, y in zip(weights, q, strict=True) + ) + target = sum( + Fraction(x) * Fraction(y) + for x, y in zip(weights, base, strict=True) + ) + assert float(abs(applied - target)) <= RTOL * max(float(target), ATOL) + assert log.loc[0, "clipping_residual"] == 0 + + +def test_cascading_caps_match_the_oracle(): + """Each capping round pushes the next row over one.""" + base = np.full(6, 0.9) + risk = np.asarray([1, 1.1, 1.3, 1.7, 2.5, 5.0]) + weights = np.asarray([3.0, 1.0, 2.0, 1.0, 0.5, 1.0]) + fixed = np.zeros(6, dtype=bool) + cells = [(2021, "male", 90)] * 6 + q, log = normalize_probabilities(base, risk, fixed, base, cells, weights) + expected, infeasible, _ = exact_cell(base, risk, fixed, base, weights) + np.testing.assert_allclose(q, expected, rtol=1e-12, atol=1e-15) + assert log.loc[0, "clipped_persons"] >= 3 + assert not infeasible + assert_conserved(q, base, cells, weights, log) + + +@PROPERTY +@given( + st.lists(st.floats(0.01, 1), min_size=3, max_size=6), + st.lists(POSITIVE, min_size=3, max_size=6), + st.floats(0, 1), + st.lists(POSITIVE, min_size=2, max_size=5), + st.booleans(), +) +def test_weightless_cells_take_the_vanishing_weight_limit( + base, risk, fixed_q, free_weights, equal +): + """Infeasible cells: the limit is unique. Feasible: unit weights.""" + size = min(len(base), len(risk), len(free_weights) + 1) + base, risk = np.asarray(base[:size]), np.asarray(risk[:size]) + fixed = np.arange(size) == 0 + prescribed = base.copy() + prescribed[0] = fixed_q + cells = [(2021, "female", 50)] * size + vanishing = np.asarray(free_weights[: size - 1]) + if equal: + vanishing = np.full(size - 1, vanishing[0]) + zero, zero_log = normalize_probabilities( + base, risk, fixed, prescribed, cells, np.r_[1.0, np.zeros(size - 1)] + ) + small, _ = normalize_probabilities( + base, risk, fixed, prescribed, cells, np.r_[1.0, 1e-60 * vanishing] + ) + if zero_log.loc[0, "infeasible"]: + np.testing.assert_array_equal(zero, small) + elif equal: + assert zero_log.loc[0, "normalization_basis"] == "persons" + np.testing.assert_allclose(zero, small, rtol=RTOL, atol=ATOL) @st.composite @@ -256,6 +486,19 @@ def test_weightless_cells_count_persons_equally(inputs): np.testing.assert_array_equal(zero_log["factor"], unit_log["factor"]) assert not zero_log["infeasible"].any() assert (zero_log["clipping_residual"] == 0).all() + assert (unit_log["normalization_basis"] == "weights").all() + for row in zero_log.itertuples(index=False): + mask = np.asarray( + [cell == (row.year, row.sex, row.age) for cell in cells] + ) + support = mask & ~fixed & (base > 0) + if support.any(): + assert row.normalization_basis == "persons" + assert row.persons_target == pytest.approx( + np.sum(base[support]), rel=RTOL, abs=ATOL + ) + else: + assert row.normalization_basis == "weights" @PROPERTY @@ -334,7 +577,11 @@ def invented_adapter_cells(draw): @PROPERTY @given(invented_adapter_cells()) def test_unit_risk_matches_adapter_on_every_unclipped_cell(inputs): - """Includes zero weights, infeasible cells and weightless free rows.""" + """Includes zero weights, infeasible cells and weightless free rows. + + The one intended exception besides clipping is an infeasible cell + whose free support weighs nothing. + """ base, fixed, prescribed, weights, sex, age = inputs cells = [(2021, SEXES[s], int(a)) for s, a in zip(sex, age, strict=True)] q, log = normalize_probabilities( @@ -349,17 +596,22 @@ def test_unit_risk_matches_adapter_on_every_unclipped_cell(inputs): mask = np.asarray( [cell == (row.year, row.sex, row.age) for cell in cells] ) + support = mask & ~fixed & (base > 0) + if row.infeasible and support.any() and not weights[support].any(): + # Intended: the adapter leaves weightless rows at factor one, + # while the brief's infeasible-cell rules hold for any weight. + continue np.testing.assert_allclose(q[mask], old[mask], rtol=RTOL, atol=ATOL) -def test_weightless_free_rows_in_infeasible_cells_keep_unit_weights(): - """Fixed rows carry the weight; no free factor can conserve the cell.""" +def test_weightless_infeasible_cells_follow_the_weighted_rules(): + """Unit weights apply only where the weighted target is zero.""" cells = [(2021, "female", 50)] * 3 fixed = np.asarray([True, False, False]) - for prescribed_q, sign in ((0.6, -1), (0.1, 1)): - base = np.asarray([0.3, 0.2, 0.4]) + base = np.asarray([0.3, 0.2, 0.4]) + weights = np.asarray([2.0, 0.0, 0.0]) + for prescribed_q, free_q, sign in ((0.6, 0, -1), (0.1, 1, 1)): prescribed = np.asarray([prescribed_q, 0.2, 0.4]) - weights = np.asarray([2.0, 0.0, 0.0]) q, log = normalize_probabilities( base, np.asarray([1.0, 1.0, 3.0]), @@ -368,13 +620,10 @@ def test_weightless_free_rows_in_infeasible_cells_keep_unit_weights(): cells, weights, ) - unit, _ = normalize_probabilities( - base, np.asarray([1.0, 1.0, 3.0]), fixed, base, cells, np.ones(3) - ) assert log.loc[0, "infeasible"] + assert log.loc[0, "normalization_basis"] == "weights" assert np.sign(log.loc[0, "clipping_residual"]) == sign - assert q[0] == prescribed_q - np.testing.assert_array_equal(q[1:], unit[1:]) + np.testing.assert_array_equal(q, [prescribed_q, free_q, free_q]) assert_conserved(q, base, cells, weights, log) @@ -594,6 +843,113 @@ def test_composition_conserves_cells_preserves_ratios_and_fixes_di_q( ) +def composed_accounting(q_base, risk, multiplier, weights): + """Run the model under the DI adapter and check each accounting link. + + A multiplier below one gives DI-origin persons less than their share + of deaths, so the adapter scales everyone else up and can clip. + """ + frame = pd.concat( + [ + invented_frame(6, di=True).assign(sex=sex, age=age) + for sex in SEXES + for age in (50, 70) + ], + ignore_index=True, + ) + frame["person_id"] = np.arange(len(frame)) + frame["birth_year"] = 2020 - frame["age"] + frame["weight"] = weights + context = invented_context(frame) + baseline = invented_baseline(q_base) + rates = invented_rates(multiplier=multiplier) + retained = {} + model = DifferentialMortality( + baseline, + invented_spec(risk=risk), + weight_column="weight", + di_configuration=rates, + cell_log=retained, + ) + deaths, cells = {}, {} + apply_di_aware_mortality( + frame, + context, + InventedZeroUniforms(), + population_model=model, + rates=rates, + weight_column="weight", + death_log=deaths, + cell_log=cells, + ) + # The adapter keeps only probabilities; the model's own log survives. + _, mechanism = model.evaluate(frame, context) + pd.testing.assert_frame_equal(retained[2021], mechanism) + logged = deaths[2021].set_index("person_id").loc[frame["person_id"]] + q_population = logged["q_population"].to_numpy() + q_applied = logged["q_applied"].to_numpy() + w = frame["weight"].to_numpy() + fixed = ( + frame["di_entitled"].to_numpy() + | frame["di_conversion_year"].notna().to_numpy() + ) + np.testing.assert_array_equal( + q_population[fixed], baseline(frame, context)[fixed] + ) + adapter = cells[2021].set_index(["sex", "band_lower"]) + for row in mechanism.itertuples(index=False): + cell = ((frame["sex"] == row.sex) & (frame["age"] == row.age)).values + # Fixed rows sit at baseline, so SM1 conserves the baseline cell. + assert not row.infeasible + assert row.clipping_residual == 0 + assert np.sum(w[cell] * q_population[cell]) == pytest.approx( + row.baseline_expected_deaths, rel=RTOL, abs=ATOL + ) + # The adapter starts from SM1's deaths and logs its own result. + netted = adapter.loc[(row.sex, row.age)] + assert netted.expected_deaths_population_model == pytest.approx( + row.applied_expected_deaths, rel=RTOL, abs=ATOL + ) + assert np.sum(w[cell] * q_applied[cell]) == pytest.approx( + netted.expected_deaths_applied, rel=RTOL, abs=ATOL + ) + free = cell & ~fixed + if not (q_applied[free] == 1).any(): + # Without the adapter's one-pass clip, the composition conserves. + assert netted.expected_deaths_applied == pytest.approx( + row.baseline_expected_deaths, rel=RTOL, abs=ATOL + ) + kept = free & (q_population < 1) & (q_applied < 1) + if kept.sum() >= 2: + np.testing.assert_allclose( + q_applied[kept] / q_applied[kept][0], + q_population[kept] / q_population[kept][0], + rtol=RTOL, + atol=ATOL, + ) + return mechanism, q_applied, fixed + + +@PROPERTY +@given( + st.floats(0.3, 0.9), + st.floats(1.0, 4.0), + st.floats(0.05, 1.0), + st.lists(st.floats(0.1, 10), min_size=24, max_size=24), +) +def test_composition_with_unequal_weights_keeps_every_accounting_link( + q_base, risk, multiplier, weights +): + composed_accounting(q_base, risk, multiplier, weights) + + +def test_composition_accounts_for_clipping_in_both_steps(): + weights = np.linspace(0.5, 6, 24) + mechanism, q_applied, fixed = composed_accounting(0.8, 3.0, 0.1, weights) + assert mechanism["clipped_persons"].sum() > 0 + assert (q_applied[~fixed] == 1).any() + + @PROPERTY @given(st.floats(1, 20), st.sampled_from(("lifetime_earnings", "education"))) def test_ordered_earnings_and_education_risks_preserve_probability_order( @@ -650,14 +1006,31 @@ def test_model_permutation_scaling_and_no_randomness(scale, seed): def test_model_gives_weightless_cells_their_unit_weight_probabilities(): frame = invented_frame(12) context = invented_context(frame) + logs = {} model = DifferentialMortality( - invented_baseline(0.2), invented_spec(), weight_column="weight" + invented_baseline(0.2), + invented_spec(), + weight_column="weight", + cell_log=logs, ) unit = model(frame, context) + assert (logs[2021]["normalization_basis"] == "weights").all() weightless = model(frame.assign(weight=0.0), context) np.testing.assert_array_equal(weightless, unit) assert (weightless > 0).all() assert len(np.unique(weightless)) > 1 + assert (logs[2021]["normalization_basis"] == "persons").all() + assert (logs[2021]["persons_target"] > 0).all() + + +def test_model_cell_log_must_be_a_mapping(): + with pytest.raises(TypeError, match="cell_log"): + DifferentialMortality( + invented_baseline(), + invented_spec(), + weight_column=None, + cell_log=[], + ) class InventedWideningPath: From fda767c2b2a3cb1742cae8b633d4a9c8709efb54 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 16:54:25 -0400 Subject: [PATCH 07/17] SM1: tier counts for the fourteen new unit tests tests/mortality_ses now collects 110 unit tests, up from 96. Co-Authored-By: Claude Opus 5.5 --- tests/README-tiers.md | 4 ++-- tests/tier_counts.json | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/README-tiers.md b/tests/README-tiers.md index a069e786..d101c323 100644 --- a/tests/README-tiers.md +++ b/tests/README-tiers.md @@ -39,9 +39,9 @@ pytest --collect-only -q -m oracle_policyengine | tail -1 | Tier | Tests at HEAD | |---|---:| -| `unit` | 18,184 | +| `unit` | 18,198 | | `artifact` | 4,253 | | `integration_psid` | 1,450 | | `reproduction_legacy` | 520 | | `oracle_policyengine` | 308 | -| **Total** | **24,715** | +| **Total** | **24,729** | diff --git a/tests/tier_counts.json b/tests/tier_counts.json index a2679da6..531ec14a 100644 --- a/tests/tier_counts.json +++ b/tests/tier_counts.json @@ -1,7 +1,7 @@ { "schema_version": 1, "counts": { - "unit": 18184, + "unit": 18198, "artifact": 4253, "integration_psid": 1450, "reproduction_legacy": 520, From bf49d60854d123fb294829bfc2f724b6f5a86606 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 18:06:14 -0400 Subject: [PATCH 08/17] SM1: decide cap calls near one exactly, and bound the float error Round 2 of the review found a float-path cell whose largest row's scaled probability rounded to exactly one. The loop capped nothing, so rows of tiny weight took a factor ten times too small while the weighted totals still conserved. - A row within 1e-9 of one at any capping round now sends the cell to the Decimal path. The float path's error bound is about 1.3e-12, so every other cap decision is the exact one. - The remaining-target guard is now measured against the totals the target is taken from, not the target. Fixed-death and capping cancellation can no longer compound; they once reached an error of 1.8e-11. - The comment and the design note state the error bound and its derivation. The separate 1e-12 conservation check covers weighted deaths only. - Tests: - two new unit regressions; - two new oracle properties, one for cap boundaries beside tiny weights and one for cancellation near the margin; - the oracle tolerance tightens from 1e-9 to 1e-11. - The composition test now checks that applied deaths equal baseline deaths minus SM1's residual and the adapter's own clipping gap. It also checks that the gap is positive where the adapter clips. On invented million-row frames, the normalizer's time and fallback counts are unchanged. Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 59 +++++-- .../mortality_ses/normalizer.py | 35 ++-- tests/mortality_ses/test_normalizer.py | 38 +++++ tests/mortality_ses/test_properties.py | 155 +++++++++++++++++- 4 files changed, 261 insertions(+), 26 deletions(-) diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md index 9bb154a3..8feae19d 100644 --- a/docs/design/mortality_ses_sm1.md +++ b/docs/design/mortality_ses_sm1.md @@ -151,15 +151,35 @@ reduce order sensitivity. The Decimal path (precision 2,400 digits) handles extreme magnitudes that would underflow or overflow ordinary products. The float path also hands the Decimal path every decision it cannot make reliably: -- a target within 1e-3 of the free capacity, relative to that capacity; -- a target formed by subtracting fixed deaths that is within 1e-3 of zero, +- a target within m = 1e-3 of the free capacity, relative to that + capacity; +- a target formed by subtracting fixed deaths that is within m of zero, relative to the larger total; -- a capping round that leaves less than 1e-3 of the target; +- a capping round that leaves a remaining target within m of zero, + relative to the totals the target is taken from (the larger of baseline + and fixed deaths when fixed rows are off baseline, otherwise the target + itself); +- a row whose scaled probability `q_base * relative_risk * factor` lies + within 1e-9 of one, at any capping round; - a feasible result that misses conservation by more than 1e-12, relative to the cell's expected deaths. -Float reductions err by a few units of 2^-52, so above that margin the -error moves the factor by less than the 1e-12 accounting tolerance. +The margins rest on an error bound. Write u = 2^-53. Each stable sum of +nonnegative float products errs by at most about 4u of its value. A +target formed by subtraction then errs by at most about 8u of the larger +total, and a remaining target by about 12u of it. Each guarded quantity +exceeds m times that total, so its relative error is at most about +12u/m, roughly 1.3e-12, and so are the factor's and every returned +probability's. That is far below the 1e-9 margin for cap decisions, so +every cap decision the float path makes is the exact one. A row within +1e-9 of one is different: rounding can put it on the wrong side of one, +and whether it caps decides the factor of every other row. Rows with tiny +weight can then end up with very different probabilities while the +weighted totals still conserve. +Measured against the exact oracle, cells built near these margins err by +at most about 1.3e-13 on the float path. The separate 1e-12 conservation +check concerns weighted deaths only. + Probabilities returned are float64; the log accounts for those returned values. @@ -179,8 +199,14 @@ equality. Permutation checks are exact where supported by stable sums; weight scaling and adapter comparisons use the stated tolerance. An independent exact oracle in the tests, in rational arithmetic, finds the capped rows by searching over how many rows cap rather than by iterating. -The normalizer must match it to `rtol=1e-9, atol=1e-15` on cells across -float magnitudes, near-one q and arbitrary fixed probabilities. +The normalizer must match it to `rtol=1e-11, atol=1e-15`, about seven +times the float path's error bound. Three generators feed that check: +- cells across float magnitudes, near-one q and arbitrary fixed + probabilities; +- cells whose dominant row sits at the cap, a few units of rounding + either way, beside rows of tiny weight; +- cells of ordinary magnitude whose fixed deaths cancel all but a share + near m of the cell's deaths, where the float path's error is largest. Intended monotonicity applies among free rows with the same baseline q in a cell, holding other dimensions constant. This is the setting supplied @@ -192,15 +218,24 @@ The existing DI adapter caps once. Its feasible cells can consequently fall short of the target; its residual is computed from `expected_deaths_population_model - expected_deaths_applied` in its existing cell log. The SM1 log describes normalization before that step. -Tests label the normalizer's iterative redistribution versus the +So under the adapter a cell's applied deaths equal its baseline deaths +minus two separate gaps: SM1's residual, zero in every feasible cell, and +the adapter's clipping gap, its population-model deaths minus its applied +deaths. The composition conserves only where the adapter's clip does not +bind. Tests label the normalizer's iterative redistribution versus the adapter's one-pass clipping as an intended difference. A second intended difference: in an infeasible cell whose free support weighs nothing, the adapter leaves the factor at one, while the rules above zero or saturate those rows. No shared-module change is required. At the scale of a projection frame (1,000,000 invented rows in 242 -cells), `evaluate` takes about 1.3 seconds on the build host, and the -normalizer alone takes 0.4–0.6 seconds. +cells), `evaluate` takes about 1.4 seconds on the build host, and the +normalizer alone takes 0.5–0.9 seconds. No cell of those invented frames +goes to the Decimal path. In an invented frame whose baseline q follows a +Gompertz curve to one, the 42 oldest-age cells, where every free row's q +is one, fill their capacity exactly and take the Decimal path. The +normalizer then takes about 3 seconds, because each such cell holds about +4,000 invented rows; real frames have few persons at those ages. ## Design and build traceability @@ -234,12 +269,12 @@ before being characterized as a mechanism defect. | Unit-risk bit identity | Exact common-risk baseline path | `test_equal_risks_preserve_baseline_bits_with_fixed_baseline_rows`, `test_identity_is_bit_for_bit_and_normalizer_is_deterministic` | | Same deaths with/without registry, ordinary and DI steps | `step.py` and context-only wrapper | `test_identity_probabilities_and_step_deaths`, `test_di_adapter_identity`, corresponding generated death-identity properties | | Bounds | Capping and input validation | Conservation/bounds property and extreme-risk examples | -| DI composition, ratios and baseline origin q | Fixed origin in `model.py`, existing adapter | Composition examples, `test_composition_conserves_cells_preserves_ratios_and_fixes_di_q`, and, with unequal weights, clipping in both steps and the retained model log, `test_composition_with_unequal_weights_keeps_every_accounting_link` and `test_composition_accounts_for_clipping_in_both_steps` | +| DI composition, ratios and baseline origin q; applied deaths equal baseline deaths minus SM1's residual and the adapter's clipping gap | Fixed origin in `model.py`, existing adapter | Composition examples, `test_composition_conserves_cells_preserves_ratios_and_fixes_di_q`, and, with unequal weights, clipping in both steps and the retained model log, `test_composition_with_unequal_weights_keeps_every_accounting_link` and `test_composition_accounts_for_clipping_in_both_steps` | | Intended monotonicity by earnings and education | Positive cell factor, log-additive risks | `test_ordered_earnings_and_education_risks_preserve_probability_order` | | Row/weight invariance and one-risk cells | Stable reductions, exact common-risk path | Permutation/scaling examples and properties; `test_one_free_risk_group_returns_baseline_probability` | | Weightless free support with a zero target: unit-weight factor, logged as basis `persons` | Unit-weight allocation in both arithmetic paths of `normalizer.py`; `normalization_basis` and `persons_target` log columns; the model's optional `cell_log` | `test_weightless_support_counts_persons_equally`, `test_weightless_support_uses_unit_weights_on_the_decimal_path`, `test_weightless_cells_count_persons_equally`, `test_weightless_cells_take_the_vanishing_weight_limit`, `test_model_gives_weightless_cells_their_unit_weight_probabilities`; the normalizer properties draw zero weights | | Infeasible weightless cells follow the weighted rules (the vanishing-weight limit) | Infeasibility decided before the unit-weight rule | `test_fixed_excess_zeroes_weightless_free_rows_on_both_paths`, `test_fixed_shortfall_saturates_weightless_support_on_both_paths`, `test_weightless_infeasible_cells_follow_the_weighted_rules`, `test_weightless_cells_take_the_vanishing_weight_limit` | -| Float decisions too close to call go to the Decimal path | `_ambiguous`, the remaining-target check and the conservation check in `normalizer.py` | `test_roundoff_cannot_erase_a_capped_cells_deaths`, `test_rounded_target_and_capacity_tie_is_decided_exactly`, `test_cancelling_target_is_decided_exactly`; the oracle property's examples | +| Float decisions too close to call go to the Decimal path | `_ambiguous`, the remaining-target check, the cap-decision check and the conservation check in `normalizer.py` | `test_roundoff_cannot_erase_a_capped_cells_deaths`, `test_rounded_target_and_capacity_tie_is_decided_exactly`, `test_cancelling_target_is_decided_exactly`, `test_a_cap_decision_at_one_is_decided_exactly`, `test_compounded_cancellation_is_decided_exactly`; the oracle property's examples, `test_cap_boundaries_beside_tiny_weights_match_the_oracle` and `test_cancellation_near_the_margin_matches_the_oracle` | | A target that exactly fills the weighted support takes the smallest capping factor | The Decimal path's tie branch | The oracle property's two tie examples | | Differential adapter arithmetic; intended clipping and weightless-infeasible differences | Same uncapped algebra, iterative cap | Private-adapter differential examples and properties, including `test_unit_risk_matches_adapter_on_every_unclipped_cell` over zero weights, infeasible cells and arbitrary fixed probabilities; named intended divergences | | Time path and alternative path isolation | Path protocol and hold-last default | `test_pluggable_path_changes_only_named_year`, generated hold-last/refusal property | diff --git a/src/populace_dynamics/mortality_ses/normalizer.py b/src/populace_dynamics/mortality_ses/normalizer.py index 8c744b92..91569d83 100644 --- a/src/populace_dynamics/mortality_ses/normalizer.py +++ b/src/populace_dynamics/mortality_ses/normalizer.py @@ -34,12 +34,19 @@ "unclassified_reference", ) -# The float path leaves a decision to the Decimal path when it compares -# totals this close, relative to their size, or when capping leaves a -# remaining target this small. Float reductions here err by a few units -# of 2**-52; above this margin that error moves the factor by less than -# the 1e-12 accounting tolerance. +# The float path leaves a cell to the Decimal path when cancellation +# leaves too little to compute reliably: a target within this relative +# margin m of zero or of capacity, or a capping round that leaves a +# remaining target within it of zero, relative to the totals the target +# is taken from. Each total errs by at most about 4u of the totals it +# sums (u = 2**-53), so outside the margin the factor and every returned +# probability err by at most about 12u/m, roughly 1.3e-12, relative. _FLOAT_MARGIN = 1e-3 +# Under that error a row whose scaled probability is this close to one +# may cap or not, and which it does decides the factor every other row +# gets. Such a cell also goes to the Decimal path. The margin is far +# above the error bound, so every other cap decision is the exact one. +_CAP_MARGIN = 1e-9 # A feasible float cell that misses conservation by more than this, # relative to its expected deaths, is recomputed on the Decimal path. _CONSERVATION_RTOL = 1e-12 @@ -204,12 +211,15 @@ def exact(): baseline_deaths = _stable_sum(w * base) fixed_deaths = _stable_sum(w[is_fixed] * prescribed[is_fixed]) fixed_at_baseline = np.array_equal(prescribed[is_fixed], base[is_fixed]) + # ``scale`` is the size of the totals the target is taken from, which + # bounds its rounding error. if fixed_at_baseline: # Avoid cancellation of small free totals when fixed rows dominate. - target = _stable_sum(w[free] * base[free]) + target = scale = _stable_sum(w[free] * base[free]) else: target = baseline_deaths - fixed_deaths - if _ambiguous(target, max(baseline_deaths, fixed_deaths)): + scale = max(baseline_deaths, fixed_deaths) + if _ambiguous(target, scale): return exact() support = free & (base > 0) # Weightless support leaves a zero weighted target silent on the @@ -219,7 +229,7 @@ def exact(): persons_target = None if target == 0 and np.any(support) and _stable_sum(w[support]) == 0: allocation = np.ones(len(base), dtype=np.longdouble) - target = persons_target = _stable_sum(base[support]) + target = persons_target = scale = _stable_sum(base[support]) capacity = _stable_sum(allocation[support]) capped = np.zeros(len(base), dtype=bool) infeasible = False @@ -254,13 +264,18 @@ def exact(): while np.any(active): denominator = _stable_sum(allocation[active] * raw[active]) factor = remaining / denominator - next_capped = active & (raw * factor > 1) + scaled = raw * factor + # A row this close to one may be capped or not; which it is + # decides the factor every other row gets. + if np.any(active & (np.abs(scaled - 1) <= _CAP_MARGIN)): + return exact() + next_capped = active & (scaled > 1) if not np.any(next_capped): break capped |= next_capped active &= ~next_capped remaining = target - _stable_sum(allocation[capped]) - if remaining <= _FLOAT_MARGIN * target: + if remaining <= _FLOAT_MARGIN * scale: return exact() result[free] = np.minimum(1, raw[free] * factor) result[capped] = 1 diff --git a/tests/mortality_ses/test_normalizer.py b/tests/mortality_ses/test_normalizer.py index dfdf1b23..9b38df66 100644 --- a/tests/mortality_ses/test_normalizer.py +++ b/tests/mortality_ses/test_normalizer.py @@ -261,6 +261,44 @@ def test_cancelling_target_is_decided_exactly(): assert log.iloc[0].clipping_residual == 0 +def test_a_cap_decision_at_one_is_decided_exactly(): + # Rounded, the first row's scaled probability is exactly one, so float + # arithmetic caps nothing. Exactly, it caps, and the other rows share + # the small target that remains. Expected values are the exact + # oracle's. + with patch.object( + normalizer, "_decimal_cell", wraps=normalizer._decimal_cell + ) as exact: + q, log = normalize([0.5] * 3, [20, 1e-26, 1], weights=[1, 1, 1e-20]) + exact.assert_called_once() + np.testing.assert_allclose( + q, [1, 4.999995000005e-27, 0.4999995000005], rtol=1e-12 + ) + assert log.iloc[0].clipped_persons == 1 + assert not log.iloc[0].infeasible + + +def test_compounded_cancellation_is_decided_exactly(): + # Fixed deaths cancel most of the cell's deaths, and capping cancels + # most of the rest. What remains is too small, relative to the totals + # it came from, to settle in float arithmetic. + with patch.object( + normalizer, "_decimal_cell", wraps=normalizer._decimal_cell + ) as exact: + q, log = normalize( + [0.5] * 3, + [1, 1e12, 1], + fixed=[True, False, False], + prescribed=[0.49950165, 0.5, 0.5], + weights=[1, 0.0009989, 0.0000044], + ) + exact.assert_called_once() + np.testing.assert_allclose( + q, [0.49950165, 1, 0.24999999999562017], rtol=1e-12 + ) + assert not log.iloc[0].infeasible + + def test_cells_are_separate_and_sorted_and_inputs_are_unchanged(): baseline = np.array([0.3, 0.1, 0.3, 0.1]) risks = np.array([1.0, 1.0, 2.0, 2.0]) diff --git a/tests/mortality_ses/test_properties.py b/tests/mortality_ses/test_properties.py index 8b2665a1..aaa80c4f 100644 --- a/tests/mortality_ses/test_properties.py +++ b/tests/mortality_ses/test_properties.py @@ -51,6 +51,7 @@ # The exact oracle is cheap on these small cells and is the main check # of the arithmetic, so it gets more examples. ORACLE_PROPERTY = settings(max_examples=150, deadline=None) +ORACLE_RTOL = 1e-11 REF = "invented_ref" OTHER = "invented_other" UNCLASSIFIED = "invented_unclassified" @@ -332,8 +333,36 @@ def invented_wide_cells(draw): np.array([1, 0]), ) ) +@example( # A scaled probability rounded to exactly one once capped nothing. + ( + np.array([0.5, 0.5, 0.5]), + np.array([20, 1e-26, 1]), + np.zeros(3, dtype=bool), + np.array([0.5, 0.5, 0.5]), + np.array([1, 1, 1e-20]), + ) +) +@example( # Fixed-death and capping cancellation once compounded. + ( + np.array([0.5, 0.5, 0.5]), + np.array([1, 1e12, 1]), + np.array([True, False, False]), + np.array([0.49950165, 0.5, 0.5]), + np.array([1, 0.0009989, 0.0000044]), + ) +) def test_normalizer_matches_an_exact_oracle_across_magnitudes(inputs): - base, risk, fixed, prescribed, weights = inputs + assert_matches_oracle(*inputs) + + +def assert_matches_oracle(base, risk, fixed, prescribed, weights): + """The normalizer agrees with ``exact_cell`` within its error bound. + + The float path's bound is about 1.3e-12 relative (see + ``normalizer._FLOAT_MARGIN``), and the Decimal path is exact up to the + final float64 rounding, so ``ORACLE_RTOL`` leaves headroom of about + seven times. + """ cells = [(2021, "female", 50)] * len(base) q, log = normalize_probabilities( base, risk, fixed, prescribed, cells, weights @@ -341,7 +370,7 @@ def test_normalizer_matches_an_exact_oracle_across_magnitudes(inputs): expected, infeasible, basis = exact_cell( base, risk, fixed, prescribed, weights ) - np.testing.assert_allclose(q, expected, rtol=1e-9, atol=1e-15) + np.testing.assert_allclose(q, expected, rtol=ORACLE_RTOL, atol=1e-15) assert bool(log.loc[0, "infeasible"]) == infeasible assert log.loc[0, "normalization_basis"] == basis assert ((0 <= q) & (q <= 1)).all() @@ -357,6 +386,101 @@ def test_normalizer_matches_an_exact_oracle_across_magnitudes(inputs): assert log.loc[0, "clipping_residual"] == 0 +@st.composite +def invented_cap_boundary_cells(draw): + """Cells whose dominant row sits at the cap, beside tiny weights. + + Row 0's weight equals the target that rows 0 and 1 (and an optional + fixed row) supply, so its scaled probability is one up to rounding, + nudged a few units either way. Row 1 adds target but almost no + scaled mass, and the other rows carry tiny positive weights, so + whether row 0 caps decides their probabilities. + """ + q0 = draw(st.floats(0.05, 0.95)) + q1 = draw(st.floats(0.05, 1)) + w0 = draw(wide(-3, 6)) + share = w0 * (1 - q0) + base, risk, fixed, prescribed, weights = [], [], [], [], [] + if draw(st.booleans()): + # A fixed row off baseline moves part of row 1's target into a + # cancellation of baseline and fixed deaths. + qf = draw(st.floats(0.05, 0.95)) + pf = draw(st.floats(0.05, 0.95).filter(lambda p: abs(p - qf) > 0.01)) + moved = draw(st.floats(-0.9, 0.9)) * share + wf = abs(moved / (qf - pf)) + share -= wf * (qf - pf) + base, risk, fixed = [qf], [draw(wide(-2, 2))], [True] + prescribed, weights = [pf], [wf] + w1 = share / q1 + nudge = draw(st.integers(-3, 3)) + for _ in range(abs(nudge)): + w0 = float(np.nextafter(w0, np.inf if nudge > 0 else 0.0)) + extra = draw(st.integers(1, 4)) + rows = dict(min_size=extra, max_size=extra) + base += [q0, q1] + draw(st.lists(st.floats(0.01, 1), **rows)) + risk += [draw(wide(0, 3)), draw(wide(-60, -12))] + draw( + st.lists(wide(-2, 2), **rows) + ) + # Relative to row 0, from below one unit in the last place upward. + weights += [w0, w1] + [ + w0 * tiny for tiny in draw(st.lists(wide(-40, -8), **rows)) + ] + fixed += [False] * (2 + extra) + prescribed = prescribed + base[len(prescribed) :] + order = draw(st.permutations(range(len(base)))) + return tuple( + np.asarray([values[i] for i in order], dtype=dtype) + for values, dtype in ( + (base, float), + (risk, float), + (fixed, bool), + (prescribed, float), + (weights, float), + ) + ) + + +@ORACLE_PROPERTY +@given(invented_cap_boundary_cells()) +def test_cap_boundaries_beside_tiny_weights_match_the_oracle(inputs): + assert_matches_oracle(*inputs) + + +@st.composite +def invented_near_margin_cells(draw): + """Ordinary magnitudes with cancellation near the float margin. + + A fixed row's deaths cancel all but a drawn share of the cell's + deaths, near the margin on either side, and large risks cap free rows, + so capping cancels much of the rest. Most such cells stay on the + float path, where the error bound is tightest. + """ + size = draw(st.integers(1, 7)) + rows = dict(min_size=size, max_size=size) + free_base = np.asarray(draw(st.lists(st.floats(0.01, 1), **rows))) + free_risk = draw(st.lists(st.one_of(wide(-3, 0), wide(0, 6)), **rows)) + free_weights = np.asarray(draw(st.lists(wide(-4, 4), **rows))) + q0 = draw(st.floats(0.01, 0.5)) + free_deaths = float(np.sum(free_weights * free_base)) + # Large enough that the prescribed probability below stays in [0, 1]. + w0 = free_deaths * draw(st.floats(1, 3)) / (1 - q0) + keep = draw(wide(-3.5, -1)) + p0 = min(1.0, (w0 * q0 + free_deaths) * (1 - keep) / w0) + return ( + np.r_[q0, free_base], + np.r_[draw(wide(-2, 2)), free_risk], + np.r_[True, np.zeros(size, dtype=bool)], + np.r_[p0, free_base], + np.r_[w0, free_weights], + ) + + +@ORACLE_PROPERTY +@given(invented_near_margin_cells()) +def test_cancellation_near_the_margin_matches_the_oracle(inputs): + assert_matches_oracle(*inputs) + + def test_cascading_caps_match_the_oracle(): """Each capping round pushes the next row over one.""" base = np.full(6, 0.9) @@ -913,6 +1037,18 @@ def composed_accounting(q_base, risk, multiplier, weights): assert np.sum(w[cell] * q_applied[cell]) == pytest.approx( netted.expected_deaths_applied, rel=RTOL, abs=ATOL ) + # Baseline deaths reach applied deaths through two separate gaps: + # SM1's residual and the adapter's own clipping gap, its + # population-model deaths minus its applied deaths. + adapter_gap = ( + netted.expected_deaths_population_model + - netted.expected_deaths_applied + ) + assert np.sum(w[cell] * q_applied[cell]) == pytest.approx( + row.baseline_expected_deaths - row.clipping_residual - adapter_gap, + rel=RTOL, + abs=ATOL, + ) free = cell & ~fixed if not (q_applied[free] == 1).any(): # Without the adapter's one-pass clip, the composition conserves. @@ -927,7 +1063,7 @@ def composed_accounting(q_base, risk, multiplier, weights): rtol=RTOL, atol=ATOL, ) - return mechanism, q_applied, fixed + return mechanism, q_applied, fixed, cells[2021] @PROPERTY @@ -945,9 +1081,20 @@ def test_composition_with_unequal_weights_keeps_every_accounting_link( def test_composition_accounts_for_clipping_in_both_steps(): weights = np.linspace(0.5, 6, 24) - mechanism, q_applied, fixed = composed_accounting(0.8, 3.0, 0.1, weights) + mechanism, q_applied, fixed, adapter = composed_accounting( + 0.8, 3.0, 0.1, weights + ) assert mechanism["clipped_persons"].sum() > 0 assert (q_applied[~fixed] == 1).any() + # SM1 conserves, but the adapter's one-pass clip loses deaths in a + # feasible cell: the composition does not conserve there. + gap = ( + adapter["expected_deaths_population_model"] + - adapter["expected_deaths_applied"] + ) + assert (gap > 1e-6).any() + assert not adapter["infeasible"].any() + assert (mechanism["clipping_residual"] == 0).all() @PROPERTY From cee73cb780f8fb64c50e8ac101cac7030e76451a Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 18:06:15 -0400 Subject: [PATCH 09/17] SM1: tier counts for the four new unit tests Unit 18,198 to 18,202; total 24,729 to 24,733. These counts are against the PR's base and will be recounted on fresh master at the tier-count slot. Co-Authored-By: Claude Opus 5.5 --- tests/README-tiers.md | 4 ++-- tests/tier_counts.json | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/README-tiers.md b/tests/README-tiers.md index d101c323..961de999 100644 --- a/tests/README-tiers.md +++ b/tests/README-tiers.md @@ -39,9 +39,9 @@ pytest --collect-only -q -m oracle_policyengine | tail -1 | Tier | Tests at HEAD | |---|---:| -| `unit` | 18,198 | +| `unit` | 18,202 | | `artifact` | 4,253 | | `integration_psid` | 1,450 | | `reproduction_legacy` | 520 | | `oracle_policyengine` | 308 | -| **Total** | **24,729** | +| **Total** | **24,733** | diff --git a/tests/tier_counts.json b/tests/tier_counts.json index 531ec14a..9264576c 100644 --- a/tests/tier_counts.json +++ b/tests/tier_counts.json @@ -1,7 +1,7 @@ { "schema_version": 1, "counts": { - "unit": 18198, + "unit": 18202, "artifact": 4253, "integration_psid": 1450, "reproduction_legacy": 520, From ad827f1ed3bfd2e2d5c374707c79f0ce71470c5f Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 18:08:32 -0400 Subject: [PATCH 10/17] SM1: state the measured near-margin error from the final search The design note now gives the largest float-path error from 10,000 invented near-margin cells at cee73cb7: 2.4e-13 relative, under the 1.3e-12 bound. The earlier figure came from a search run before the cap-decision margin was set. Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md index 8feae19d..1af2ab34 100644 --- a/docs/design/mortality_ses_sm1.md +++ b/docs/design/mortality_ses_sm1.md @@ -176,8 +176,8 @@ every cap decision the float path makes is the exact one. A row within and whether it caps decides the factor of every other row. Rows with tiny weight can then end up with very different probabilities while the weighted totals still conserve. -Measured against the exact oracle, cells built near these margins err by -at most about 1.3e-13 on the float path. The separate 1e-12 conservation +Measured against the exact oracle in 10,000 invented cells built near +these margins, the float path erred by at most 2.4e-13. The separate 1e-12 conservation check concerns weighted deaths only. Probabilities returned are float64; the log accounts for those returned From 8edc9d90ad3fd011c87d5a907d812db5d95fc076 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 21:18:26 -0400 Subject: [PATCH 11/17] SM1: count zero-weight clips exactly, and check clip counts against the oracle Round 3 found that the float path's cap-decision check examined only the rows that set the factor. A zero-weight row takes the factor without setting it, so one whose exact scaled probability lay a rounding step above one could go uncounted in the log's clipped count, though its probability stayed within the error bound. After the last capping round, any uncapped free row within the 1e-9 cap margin of one now sends the cell to the Decimal path. The exact oracle now returns the clipped count, and the oracle check compares it exactly, with the reviewer's cell as an example and a new generator that places zero-weight rows a few rounding steps either side of one, sometimes after a round that caps a dominant row. At ad827f1e that generator finds 605 miscounted cells in 2,000; here it finds none in 4,000, and the other three generators find no disagreement in 24,000 cells. On a million invented rows the fallback and clip counts are unchanged. Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 18 ++-- .../mortality_ses/normalizer.py | 10 ++- tests/mortality_ses/test_normalizer.py | 14 +++ tests/mortality_ses/test_properties.py | 87 +++++++++++++++++-- 4 files changed, 113 insertions(+), 16 deletions(-) diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md index 1af2ab34..3a7ca0ed 100644 --- a/docs/design/mortality_ses_sm1.md +++ b/docs/design/mortality_ses_sm1.md @@ -161,6 +161,9 @@ cannot make reliably: itself); - a row whose scaled probability `q_base * relative_risk * factor` lies within 1e-9 of one, at any capping round; +- after the last round, an uncapped free row within 1e-9 of one. This + catches zero-weight rows, which take the factor without setting it, so + their clipping changes only the log's clipped count; - a feasible result that misses conservation by more than 1e-12, relative to the cell's expected deaths. @@ -171,7 +174,8 @@ total, and a remaining target by about 12u of it. Each guarded quantity exceeds m times that total, so its relative error is at most about 12u/m, roughly 1.3e-12, and so are the factor's and every returned probability's. That is far below the 1e-9 margin for cap decisions, so -every cap decision the float path makes is the exact one. A row within +every cap decision and clip count the float path makes is the exact one. +A row within 1e-9 of one is different: rounding can put it on the wrong side of one, and whether it caps decides the factor of every other row. Rows with tiny weight can then end up with very different probabilities while the @@ -199,14 +203,18 @@ equality. Permutation checks are exact where supported by stable sums; weight scaling and adapter comparisons use the stated tolerance. An independent exact oracle in the tests, in rational arithmetic, finds the capped rows by searching over how many rows cap rather than by iterating. -The normalizer must match it to `rtol=1e-11, atol=1e-15`, about seven -times the float path's error bound. Three generators feed that check: +The normalizer must match its probabilities to `rtol=1e-11, atol=1e-15`, +about seven times the float path's error bound, and its clipped count, +infeasibility flag and normalization basis exactly. Four generators feed +that check: - cells across float magnitudes, near-one q and arbitrary fixed probabilities; - cells whose dominant row sits at the cap, a few units of rounding either way, beside rows of tiny weight; - cells of ordinary magnitude whose fixed deaths cancel all but a share - near m of the cell's deaths, where the float path's error is largest. + near m of the cell's deaths, where the float path's error is largest; +- cells whose zero-weight rows sit a few units of rounding either side of + one, sometimes after a first round that caps a dominant row. Intended monotonicity applies among free rows with the same baseline q in a cell, holding other dimensions constant. This is the setting supplied @@ -274,7 +282,7 @@ before being characterized as a mechanism defect. | Row/weight invariance and one-risk cells | Stable reductions, exact common-risk path | Permutation/scaling examples and properties; `test_one_free_risk_group_returns_baseline_probability` | | Weightless free support with a zero target: unit-weight factor, logged as basis `persons` | Unit-weight allocation in both arithmetic paths of `normalizer.py`; `normalization_basis` and `persons_target` log columns; the model's optional `cell_log` | `test_weightless_support_counts_persons_equally`, `test_weightless_support_uses_unit_weights_on_the_decimal_path`, `test_weightless_cells_count_persons_equally`, `test_weightless_cells_take_the_vanishing_weight_limit`, `test_model_gives_weightless_cells_their_unit_weight_probabilities`; the normalizer properties draw zero weights | | Infeasible weightless cells follow the weighted rules (the vanishing-weight limit) | Infeasibility decided before the unit-weight rule | `test_fixed_excess_zeroes_weightless_free_rows_on_both_paths`, `test_fixed_shortfall_saturates_weightless_support_on_both_paths`, `test_weightless_infeasible_cells_follow_the_weighted_rules`, `test_weightless_cells_take_the_vanishing_weight_limit` | -| Float decisions too close to call go to the Decimal path | `_ambiguous`, the remaining-target check, the cap-decision check and the conservation check in `normalizer.py` | `test_roundoff_cannot_erase_a_capped_cells_deaths`, `test_rounded_target_and_capacity_tie_is_decided_exactly`, `test_cancelling_target_is_decided_exactly`, `test_a_cap_decision_at_one_is_decided_exactly`, `test_compounded_cancellation_is_decided_exactly`; the oracle property's examples, `test_cap_boundaries_beside_tiny_weights_match_the_oracle` and `test_cancellation_near_the_margin_matches_the_oracle` | +| Float decisions too close to call go to the Decimal path | `_ambiguous`, the remaining-target check, the cap-decision check and the conservation check in `normalizer.py` | `test_roundoff_cannot_erase_a_capped_cells_deaths`, `test_rounded_target_and_capacity_tie_is_decided_exactly`, `test_cancelling_target_is_decided_exactly`, `test_a_cap_decision_at_one_is_decided_exactly`, `test_compounded_cancellation_is_decided_exactly`, `test_a_zero_weight_row_at_one_is_counted_exactly`; the oracle property's examples, `test_cap_boundaries_beside_tiny_weights_match_the_oracle`, `test_cancellation_near_the_margin_matches_the_oracle` and `test_weightless_rows_at_one_match_the_oracle` | | A target that exactly fills the weighted support takes the smallest capping factor | The Decimal path's tie branch | The oracle property's two tie examples | | Differential adapter arithmetic; intended clipping and weightless-infeasible differences | Same uncapped algebra, iterative cap | Private-adapter differential examples and properties, including `test_unit_risk_matches_adapter_on_every_unclipped_cell` over zero weights, infeasible cells and arbitrary fixed probabilities; named intended divergences | | Time path and alternative path isolation | Path protocol and hold-last default | `test_pluggable_path_changes_only_named_year`, generated hold-last/refusal property | diff --git a/src/populace_dynamics/mortality_ses/normalizer.py b/src/populace_dynamics/mortality_ses/normalizer.py index 91569d83..52ce49fe 100644 --- a/src/populace_dynamics/mortality_ses/normalizer.py +++ b/src/populace_dynamics/mortality_ses/normalizer.py @@ -44,8 +44,9 @@ _FLOAT_MARGIN = 1e-3 # Under that error a row whose scaled probability is this close to one # may cap or not, and which it does decides the factor every other row -# gets. Such a cell also goes to the Decimal path. The margin is far -# above the error bound, so every other cap decision is the exact one. +# gets, or, for a zero-allocation row, the log's clipped count. Such a +# cell also goes to the Decimal path. The margin is far above the error +# bound, so every other cap decision and clip count is the exact one. _CAP_MARGIN = 1e-9 # A feasible float cell that misses conservation by more than this, # relative to its expected deaths, is recomputed on the Decimal path. @@ -277,6 +278,11 @@ def exact(): remaining = target - _stable_sum(allocation[capped]) if remaining <= _FLOAT_MARGIN * scale: return exact() + # The loop checked every row that sets the factor. Zero-allocation + # rows take the factor without setting it, so whether one clips + # changes only the log's count, but that count must be exact too. + if np.any(free & ~capped & (np.abs(raw * factor - 1) <= _CAP_MARGIN)): + return exact() result[free] = np.minimum(1, raw[free] * factor) result[capped] = 1 capped |= free & (raw * factor > 1) diff --git a/tests/mortality_ses/test_normalizer.py b/tests/mortality_ses/test_normalizer.py index 9b38df66..593866d1 100644 --- a/tests/mortality_ses/test_normalizer.py +++ b/tests/mortality_ses/test_normalizer.py @@ -278,6 +278,20 @@ def test_a_cap_decision_at_one_is_decided_exactly(): assert not log.iloc[0].infeasible +def test_a_zero_weight_row_at_one_is_counted_exactly(): + # The weightless third row takes the factor without setting it. + # Exactly, its scaled probability exceeds one by about 6e-17, so it + # clips; rounded, it lands on one and would not be counted. + with patch.object( + normalizer, "_decimal_cell", wraps=normalizer._decimal_cell + ) as exact: + q, log = normalize([0.25, 0.5, 0.5], [1, 3, 14 / 3], weights=[1, 1, 0]) + exact.assert_called_once() + np.testing.assert_allclose(q, [3 / 28, 9 / 14, 1], rtol=1e-15) + assert log.iloc[0].clipped_persons == 1 + assert not log.iloc[0].infeasible + + def test_compounded_cancellation_is_decided_exactly(): # Fixed deaths cancel most of the cell's deaths, and capping cancels # most of the rest. What remains is too small, relative to the totals diff --git a/tests/mortality_ses/test_properties.py b/tests/mortality_ses/test_properties.py index aaa80c4f..bbea890e 100644 --- a/tests/mortality_ses/test_properties.py +++ b/tests/mortality_ses/test_properties.py @@ -188,8 +188,8 @@ def exact_cell(base, risk, fixed, prescribed, weights): It uses rational arithmetic and finds the capped rows by trying each count of largest-q*RR rows in turn, rather than iterating the factor - as the normalizer does. Returns probabilities, the infeasible flag - and the normalization basis. + as the normalizer does. Returns probabilities, the infeasible flag, + the normalization basis and the number of free rows clipped at one. """ b, r, p, w = ( [Fraction(float(x)) for x in values] @@ -208,19 +208,23 @@ def exact_cell(base, risk, fixed, prescribed, weights): target = sum(b[i] for i in support) capacity = sum((allocation[i] for i in support), Fraction(0)) if target <= 0: - return [float(x) for x in q], target < 0, basis + return [float(x) for x in q], target < 0, basis, 0 if target > capacity: for i in support: q[i] = Fraction(1) - return [float(x) for x in q], True, basis + return [float(x) for x in q], True, basis, len(support) raw = [b[i] * r[i] for i in rows] + + def clipped(factor): + return sum(raw[i] * factor > 1 for i in free) + if target == capacity: # Every weighted row must reach one; zero-allocation rows take the # smallest factor that does that. factor = max(1 / raw[i] for i in support if allocation[i] > 0) for i in free: q[i] = min(Fraction(1), raw[i] * factor) - return [float(x) for x in q], False, basis + return [float(x) for x in q], False, basis, clipped(factor) ranked = sorted( (i for i in support if allocation[i] > 0), key=lambda i: raw[i], @@ -237,7 +241,7 @@ def exact_cell(base, risk, fixed, prescribed, weights): raise AssertionError("no capped set satisfies the cell") for i in free: q[i] = min(Fraction(1), raw[i] * factor) - return [float(x) for x in q], False, basis + return [float(x) for x in q], False, basis, clipped(factor) def wide(low, high): @@ -351,6 +355,15 @@ def invented_wide_cells(draw): np.array([1, 0.0009989, 0.0000044]), ) ) +@example( # A weightless row at one once went uncounted as clipped. + ( + np.array([0.25, 0.5, 0.5]), + np.array([1, 3, 14 / 3]), + np.zeros(3, dtype=bool), + np.array([0.25, 0.5, 0.5]), + np.array([1, 1, 0]), + ) +) def test_normalizer_matches_an_exact_oracle_across_magnitudes(inputs): assert_matches_oracle(*inputs) @@ -367,12 +380,13 @@ def assert_matches_oracle(base, risk, fixed, prescribed, weights): q, log = normalize_probabilities( base, risk, fixed, prescribed, cells, weights ) - expected, infeasible, basis = exact_cell( + expected, infeasible, basis, clipped = exact_cell( base, risk, fixed, prescribed, weights ) np.testing.assert_allclose(q, expected, rtol=ORACLE_RTOL, atol=1e-15) assert bool(log.loc[0, "infeasible"]) == infeasible assert log.loc[0, "normalization_basis"] == basis + assert log.loc[0, "clipped_persons"] == clipped assert ((0 <= q) & (q <= 1)).all() if not infeasible: applied = sum( @@ -481,6 +495,59 @@ def test_cancellation_near_the_margin_matches_the_oracle(inputs): assert_matches_oracle(*inputs) +@st.composite +def invented_weightless_cap_cells(draw): + """Weightless rows within a few units of rounding of one. + + Weighted rows set the factor, optionally after a first round that caps + one dominant row; none of the others caps. Their exact factor is a + ratio of rational sums, so each weightless row's risk can be placed a + drawn number of float steps from where its scaled probability is + exactly one, on either side. The rows take that factor without + setting it, so they test only the clipped count and their own + probabilities. + """ + size = draw(st.integers(1, 5)) + rows = dict(min_size=size, max_size=size) + base = draw(st.lists(st.floats(0.01, 0.2), **rows)) + risk = draw(st.lists(st.floats(0.5, 2), **rows)) + weights = draw(st.lists(st.floats(0.1, 10), **rows)) + deaths = sum( + Fraction(w) * Fraction(b) for w, b in zip(weights, base, strict=True) + ) + scaled = sum( + Fraction(w) * Fraction(b) * Fraction(r) + for w, b, r in zip(weights, base, risk, strict=True) + ) + if draw(st.booleans()): + # This row caps in the first round; the rest share what remains. + dominant = float(deaths) * draw(st.floats(0.1, 1.5)) + base, risk = [0.5, *base], [1000.0, *risk] + weights = [dominant, *weights] + deaths -= Fraction(dominant) / 2 + factor = deaths / scaled + for _ in range(draw(st.integers(1, 3))): + q0 = draw(st.floats(0.05, 0.95)) + r0 = float(1 / (Fraction(q0) * factor)) + for _ in range(abs(steps := draw(st.integers(-4, 4)))): + r0 = float(np.nextafter(r0, np.inf if steps > 0 else 0)) + base, risk, weights = [*base, q0], [*risk, r0], [*weights, 0.0] + base = np.asarray(base) + return ( + base, + np.asarray(risk), + np.zeros(len(base), dtype=bool), + base, + np.asarray(weights), + ) + + +@ORACLE_PROPERTY +@given(invented_weightless_cap_cells()) +def test_weightless_rows_at_one_match_the_oracle(inputs): + assert_matches_oracle(*inputs) + + def test_cascading_caps_match_the_oracle(): """Each capping round pushes the next row over one.""" base = np.full(6, 0.9) @@ -489,9 +556,11 @@ def test_cascading_caps_match_the_oracle(): fixed = np.zeros(6, dtype=bool) cells = [(2021, "male", 90)] * 6 q, log = normalize_probabilities(base, risk, fixed, base, cells, weights) - expected, infeasible, _ = exact_cell(base, risk, fixed, base, weights) + expected, infeasible, _, clipped = exact_cell( + base, risk, fixed, base, weights + ) np.testing.assert_allclose(q, expected, rtol=1e-12, atol=1e-15) - assert log.loc[0, "clipped_persons"] >= 3 + assert log.loc[0, "clipped_persons"] == clipped >= 3 assert not infeasible assert_conserved(q, base, cells, weights, log) From 8ca39c26d4f125d1083b8292797e1db3aa91808e Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 21:18:26 -0400 Subject: [PATCH 12/17] SM1: tier counts for the two new unit tests Unit 18,202 to 18,204; total 24,733 to 24,735. These counts are against the PR's base and will be recounted on fresh master at the tier-count slot. Co-Authored-By: Claude Opus 5.5 --- tests/README-tiers.md | 4 ++-- tests/tier_counts.json | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/README-tiers.md b/tests/README-tiers.md index 961de999..bf3e39c6 100644 --- a/tests/README-tiers.md +++ b/tests/README-tiers.md @@ -39,9 +39,9 @@ pytest --collect-only -q -m oracle_policyengine | tail -1 | Tier | Tests at HEAD | |---|---:| -| `unit` | 18,202 | +| `unit` | 18,204 | | `artifact` | 4,253 | | `integration_psid` | 1,450 | | `reproduction_legacy` | 520 | | `oracle_policyengine` | 308 | -| **Total** | **24,733** | +| **Total** | **24,735** | diff --git a/tests/tier_counts.json b/tests/tier_counts.json index 9264576c..bdf7245c 100644 --- a/tests/tier_counts.json +++ b/tests/tier_counts.json @@ -1,7 +1,7 @@ { "schema_version": 1, "counts": { - "unit": 18202, + "unit": 18204, "artifact": 4253, "integration_psid": 1450, "reproduction_legacy": 520, From e1155b4bbb92b43cc78a76acde9b28c6e8075fda Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 21:49:32 -0400 Subject: [PATCH 13/17] SM1: define the clipped count, and state reproducible search figures The delta review of 8edc9d90..8ca39c26 approved it, with two wording findings and a nit. The design note and the normalizer's docstring now define `clipped_persons`: free rows whose scaled probability exceeds one at the cell's final factor, not counting a row exactly at one; in a cell whose target exceeds its free capacity, every free row with positive baseline q. The note's measured error now comes from a derandomized search (1.8e-13 relative, near the margin), with the unseeded range beside it, and a paragraph is rewrapped. The search figures in 8edc9d90's message came from unseeded runs. In derandomized runs, ad827f1e miscounts 1,174 of the same 4,000 zero-weight cells this head gets right, and the other three generators give identical results at both heads. No code path changes here. Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 21 ++++++++++++------- .../mortality_ses/normalizer.py | 6 ++++++ 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md index 3a7ca0ed..c8495ba6 100644 --- a/docs/design/mortality_ses_sm1.md +++ b/docs/design/mortality_ses_sm1.md @@ -175,14 +175,14 @@ exceeds m times that total, so its relative error is at most about 12u/m, roughly 1.3e-12, and so are the factor's and every returned probability's. That is far below the 1e-9 margin for cap decisions, so every cap decision and clip count the float path makes is the exact one. -A row within -1e-9 of one is different: rounding can put it on the wrong side of one, -and whether it caps decides the factor of every other row. Rows with tiny -weight can then end up with very different probabilities while the -weighted totals still conserve. -Measured against the exact oracle in 10,000 invented cells built near -these margins, the float path erred by at most 2.4e-13. The separate 1e-12 conservation -check concerns weighted deaths only. +A row within 1e-9 of one is different: rounding can put it on the wrong +side of one, and whether it caps decides the factor of every other row. +Rows with tiny weight can then end up with very different probabilities +while the weighted totals still conserve. +Measured against the exact oracle in a derandomized search of 10,000 +invented cells built near these margins, the float path erred by at most +1.8e-13; unseeded searches of the same generator gave up to 2.4e-13. The +separate 1e-12 conservation check concerns weighted deaths only. Probabilities returned are float64; the log accounts for those returned values. @@ -191,6 +191,11 @@ The cell log contains year, sex, age, persons, fixed persons, total weight, fixed weight, baseline expected deaths, applied expected deaths, factor, clipped persons, signed clipping residual, infeasibility, normalization basis and persons target, plus the unclassified counts. +`clipped_persons` counts free rows whose scaled probability +`q_base * relative_risk * factor` exceeds one at the cell's final factor; +a row exactly at one is not counted. In a cell whose target exceeds its +free capacity, it counts every free row with positive baseline q, since +all of them saturate. `clipping_residual = baseline - applied`. It is nonzero only on infeasible cells: iterative capping conserves every feasible cell, and the float path verifies that rather than assuming it. diff --git a/src/populace_dynamics/mortality_ses/normalizer.py b/src/populace_dynamics/mortality_ses/normalizer.py index 52ce49fe..83bdba83 100644 --- a/src/populace_dynamics/mortality_ses/normalizer.py +++ b/src/populace_dynamics/mortality_ses/normalizer.py @@ -360,6 +360,12 @@ def normalize_probabilities( With fixed rows at baseline, such a cell gets the probabilities it would get if every weight were one. + The log's ``clipped_persons`` counts free rows whose scaled probability + ``q_base * relative_risk * factor`` exceeds one at the cell's final + factor; a row exactly at one is not counted. In a cell whose target + exceeds its free capacity, it counts every free row with positive + baseline probability, since all of them saturate. + The signed residual is baseline minus applied deaths. It is nonzero only on infeasible cells; iterative capping conserves feasible cells. The float path verifies that within the accounting tolerance, and From 5bda0fed1d1a3f041e27dac030253d10ffafc024 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 21:59:08 -0400 Subject: [PATCH 14/17] SM1: qualify the clipped-count definition by the exact factor The confirmation of e1155b4b found that the definition holds only at the cell's exact final factor. The logged factor is its float64 rounding, so a row exactly at one, such as a q = 1 row at a terminal age, can appear to exceed one when recomputed from the log. The design note and the normalize_probabilities docstring now say so. No code changes. Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 9 +++++---- src/populace_dynamics/mortality_ses/normalizer.py | 9 +++++---- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md index c8495ba6..c918b825 100644 --- a/docs/design/mortality_ses_sm1.md +++ b/docs/design/mortality_ses_sm1.md @@ -192,10 +192,11 @@ weight, fixed weight, baseline expected deaths, applied expected deaths, factor, clipped persons, signed clipping residual, infeasibility, normalization basis and persons target, plus the unclassified counts. `clipped_persons` counts free rows whose scaled probability -`q_base * relative_risk * factor` exceeds one at the cell's final factor; -a row exactly at one is not counted. In a cell whose target exceeds its -free capacity, it counts every free row with positive baseline q, since -all of them saturate. +`q_base * relative_risk * factor` exceeds one at the cell's exact final +factor (the logged `factor` is its float64 rounding, so a row at one can +appear to exceed it); a row exactly at one is not counted. In a cell +whose target exceeds its free capacity, it counts every free row with +positive baseline q, since all of them saturate. `clipping_residual = baseline - applied`. It is nonzero only on infeasible cells: iterative capping conserves every feasible cell, and the float path verifies that rather than assuming it. diff --git a/src/populace_dynamics/mortality_ses/normalizer.py b/src/populace_dynamics/mortality_ses/normalizer.py index 83bdba83..a865c16c 100644 --- a/src/populace_dynamics/mortality_ses/normalizer.py +++ b/src/populace_dynamics/mortality_ses/normalizer.py @@ -361,10 +361,11 @@ def normalize_probabilities( would get if every weight were one. The log's ``clipped_persons`` counts free rows whose scaled probability - ``q_base * relative_risk * factor`` exceeds one at the cell's final - factor; a row exactly at one is not counted. In a cell whose target - exceeds its free capacity, it counts every free row with positive - baseline probability, since all of them saturate. + ``q_base * relative_risk * factor`` exceeds one at the cell's exact + final factor (the logged ``factor`` is its float64 rounding, so a row + at one can appear to exceed it); a row exactly at one is not counted. + In a cell whose target exceeds its free capacity, it counts every free + row with positive baseline probability, since all of them saturate. The signed residual is baseline minus applied deaths. It is nonzero only on infeasible cells; iterative capping conserves feasible cells. From 3078d18be5f7a4049dfed530d33ab2b662599cc2 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 22:32:49 -0400 Subject: [PATCH 15/17] SM1: hold the adapter differential to the exact formula The confirmation of 5bda0fed found the differential test against the DI adapter's _net_of_di_origin flaky. On 1,000 fresh CI-sized runs at 5bda0fed, 4 failed. In every failing cell SM1 equals the adapter's formula evaluated exactly; the adapter forms max(P - D, 0) / E in float64 with no guard, so its factor can cancel, absorb the free deaths or overflow. The test now holds every unclipped cell to the adapter's formula in rational arithmetic (exact_adapter_cell), at the oracle's tolerance. It compares the adapter's float64 result only where neither P - D nor E lies within a thousandth of max(P, D), with the sums taken exactly so a subnormal E is not lost. That is a third named intended divergence. Four examples cover overflow, absorption, subnormal absorption beside a weightless row, and cancellation. The adapter's intended overflow warning is silenced around its call. Over 140,000 fresh draws the test passes; at 1e-3 the worst float64 comparison it keeps uses under 3% of its tolerance. A relative error of 5e-11 in SM1 fails all four examples through the exact check alone. The clipped-count definition now says the logged factor approximates the exact one in float64, within the float path's error bound, as the same review proposed. No change to SM1's code. Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 20 ++-- .../mortality_ses/normalizer.py | 9 +- tests/mortality_ses/test_properties.py | 100 +++++++++++++++++- 3 files changed, 114 insertions(+), 15 deletions(-) diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md index c918b825..627c95ba 100644 --- a/docs/design/mortality_ses_sm1.md +++ b/docs/design/mortality_ses_sm1.md @@ -193,10 +193,11 @@ factor, clipped persons, signed clipping residual, infeasibility, normalization basis and persons target, plus the unclassified counts. `clipped_persons` counts free rows whose scaled probability `q_base * relative_risk * factor` exceeds one at the cell's exact final -factor (the logged `factor` is its float64 rounding, so a row at one can -appear to exceed it); a row exactly at one is not counted. In a cell -whose target exceeds its free capacity, it counts every free row with -positive baseline q, since all of them saturate. +factor (the logged `factor` approximates it in float64, within the stated +error bound on the float path, so a row at one can appear to exceed it); +a row exactly at one is not counted. In a cell whose target exceeds its +free capacity, it counts every free row with positive baseline q, since +all of them saturate. `clipping_residual = baseline - applied`. It is nonzero only on infeasible cells: iterative capping conserves every feasible cell, and the float path verifies that rather than assuming it. @@ -240,7 +241,14 @@ bind. Tests label the normalizer's iterative redistribution versus the adapter's one-pass clipping as an intended difference. A second intended difference: in an infeasible cell whose free support weighs nothing, the adapter leaves the factor at one, while the rules above zero or saturate -those rows. No shared-module change is required. +those rows. A third: the adapter forms its factor `max(P - D, 0) / E` in +float64 with no guard on the subtraction or the division. Where `P - D` or +the free deaths `E` lies within a thousandth of `max(P, D)`, its factor +can lose precision, absorb the free deaths or overflow. The differential +test therefore holds every unclipped cell to the adapter's formula +evaluated exactly, at the oracle's tolerance, and compares the adapter's +float64 result only outside those cells. No shared-module change is +required. At the scale of a projection frame (1,000,000 invented rows in 242 cells), `evaluate` takes about 1.4 seconds on the build host, and the @@ -290,7 +298,7 @@ before being characterized as a mechanism defect. | Infeasible weightless cells follow the weighted rules (the vanishing-weight limit) | Infeasibility decided before the unit-weight rule | `test_fixed_excess_zeroes_weightless_free_rows_on_both_paths`, `test_fixed_shortfall_saturates_weightless_support_on_both_paths`, `test_weightless_infeasible_cells_follow_the_weighted_rules`, `test_weightless_cells_take_the_vanishing_weight_limit` | | Float decisions too close to call go to the Decimal path | `_ambiguous`, the remaining-target check, the cap-decision check and the conservation check in `normalizer.py` | `test_roundoff_cannot_erase_a_capped_cells_deaths`, `test_rounded_target_and_capacity_tie_is_decided_exactly`, `test_cancelling_target_is_decided_exactly`, `test_a_cap_decision_at_one_is_decided_exactly`, `test_compounded_cancellation_is_decided_exactly`, `test_a_zero_weight_row_at_one_is_counted_exactly`; the oracle property's examples, `test_cap_boundaries_beside_tiny_weights_match_the_oracle`, `test_cancellation_near_the_margin_matches_the_oracle` and `test_weightless_rows_at_one_match_the_oracle` | | A target that exactly fills the weighted support takes the smallest capping factor | The Decimal path's tie branch | The oracle property's two tie examples | -| Differential adapter arithmetic; intended clipping and weightless-infeasible differences | Same uncapped algebra, iterative cap | Private-adapter differential examples and properties, including `test_unit_risk_matches_adapter_on_every_unclipped_cell` over zero weights, infeasible cells and arbitrary fixed probabilities; named intended divergences | +| Differential adapter arithmetic; intended clipping, weightless-infeasible and adapter float64 differences | Same uncapped algebra, iterative cap | Private-adapter differential examples and properties, including `test_unit_risk_matches_adapter_on_every_unclipped_cell` over zero weights, infeasible cells and arbitrary fixed probabilities, which also holds every unclipped cell to the adapter's formula evaluated exactly (`exact_adapter_cell`); named intended divergences | | Time path and alternative path isolation | Path protocol and hold-last default | `test_pluggable_path_changes_only_named_year`, generated hold-last/refusal property | | Determinism and no RNG consumption | Pure `evaluate` | Repeated-call examples, generator state and patched-generator property | | Refusals | Specification, frame and interface validation | Nonfinite/nonpositive risks, missing columns, unknown/unclassified values, age mismatch and `probabilities` attribute examples and properties | diff --git a/src/populace_dynamics/mortality_ses/normalizer.py b/src/populace_dynamics/mortality_ses/normalizer.py index a865c16c..0410b94d 100644 --- a/src/populace_dynamics/mortality_ses/normalizer.py +++ b/src/populace_dynamics/mortality_ses/normalizer.py @@ -362,10 +362,11 @@ def normalize_probabilities( The log's ``clipped_persons`` counts free rows whose scaled probability ``q_base * relative_risk * factor`` exceeds one at the cell's exact - final factor (the logged ``factor`` is its float64 rounding, so a row - at one can appear to exceed it); a row exactly at one is not counted. - In a cell whose target exceeds its free capacity, it counts every free - row with positive baseline probability, since all of them saturate. + final factor (the logged ``factor`` approximates it in float64, within + the stated error bound on the float path, so a row at one can appear + to exceed it); a row exactly at one is not counted. In a cell whose + target exceeds its free capacity, it counts every free row with + positive baseline probability, since all of them saturate. The signed residual is baseline minus applied deaths. It is nonzero only on infeasible cells; iterative capping conserves feasible cells. diff --git a/tests/mortality_ses/test_properties.py b/tests/mortality_ses/test_properties.py index bbea890e..b1ecd11c 100644 --- a/tests/mortality_ses/test_properties.py +++ b/tests/mortality_ses/test_properties.py @@ -767,22 +767,102 @@ def invented_adapter_cells(draw): return base, fixed, prescribed, weights, groups % 2, 50 + groups // 2 +def exact_adapter_cell(base, fixed, prescribed, weights): + """The adapter's formula for one cell, in exact arithmetic. + + ``_net_of_di_origin`` multiplies a cell's free rows by + ``max(P - D, 0) / E`` and clips at one, where ``P`` is the cell's + baseline deaths, ``D`` its fixed rows' deaths and ``E`` its free rows' + baseline deaths; with no fixed row or no free deaths the factor is one. + Also returns whether the adapter's float64 factor is unreliable: when + ``P - D`` or ``E`` is within a thousandth of ``max(P, D)``, the + subtraction can cancel, absorb ``E`` or overflow the division. Exact + sums keep an ``E`` that would underflow in float64. + """ + weight = [Fraction(w) for w in weights] + population = sum( + w * Fraction(b) for w, b in zip(weight, base, strict=True) + ) + fixed_deaths = sum( + weight[i] * Fraction(prescribed[i]) for i in np.flatnonzero(fixed) + ) + free_deaths = sum( + weight[i] * Fraction(base[i]) for i in np.flatnonzero(~fixed) + ) + factor, unreliable = Fraction(1), False + if fixed.any(): + scale = max(population, fixed_deaths) / 1000 + gap = abs(population - fixed_deaths) + unreliable = 0 < free_deaths <= scale or 0 < gap <= scale + if free_deaths > 0: + factor = max(population - fixed_deaths, 0) / free_deaths + q = [ + p if f else float(min(1, Fraction(b) * factor)) + for b, f, p in zip(base, fixed, prescribed, strict=True) + ] + return np.asarray(q, dtype=float), unreliable + + @PROPERTY @given(invented_adapter_cells()) +@example( # The adapter's factor overflows on a subnormal free row. + ( + np.array([2.225073858507e-311, 0.5]), + np.array([False, True]), + np.array([2.225073858507e-311, 0.0]), + np.array([1.0, 1.0]), + np.zeros(2, dtype=int), + np.full(2, 50), + ) +) +@example( # Its baseline deaths absorb the free deaths: factor zero. + ( + np.array([0.5, 0.8, 3e-301]), + np.array([False, True, False]), + np.array([0.5, 0.8, 3e-301]), + np.array([0.0, 4.0, 16.0]), + np.zeros(3, dtype=int), + np.full(3, 50), + ) +) +@example( # The same with subnormal free deaths, beside a weightless row. + ( + np.array([5e-324, 0.5074563290805952, 0.5537272658366658]), + np.array([False, False, True]), + np.array([5e-324, 0.5074563290805952, 0.5537272658366658]), + np.array([1.6873477097608718, 0.0, 18.517829151607195]), + np.zeros(3, dtype=int), + np.full(3, 50), + ) +) +@example( # Its subtraction cancels to an error above the tolerance. + ( + np.array([0.9, 1.192092896e-07, 0.5543192935511442]), + np.array([False, False, True]), + np.array([0.9, 1.192092896e-07, 0.5543192935511442]), + np.array([0.0, 2.00001, 0.05000000000000001]), + np.zeros(3, dtype=int), + np.full(3, 50), + ) +) def test_unit_risk_matches_adapter_on_every_unclipped_cell(inputs): """Includes zero weights, infeasible cells and weightless free rows. - The one intended exception besides clipping is an infeasible cell - whose free support weighs nothing. + Besides clipping, one intended exception skips the cell: an + infeasible cell whose free support weighs nothing. Every other cell + matches the adapter's formula evaluated exactly, at the oracle's + tolerance, and the adapter's float64 result except where that is + unreliable (``exact_adapter_cell``), a second intended exception. """ base, fixed, prescribed, weights, sex, age = inputs cells = [(2021, SEXES[s], int(a)) for s, a in zip(sex, age, strict=True)] q, log = normalize_probabilities( base, np.ones(len(base)), fixed, prescribed, cells, weights ) - old, _ = _net_of_di_origin( - prescribed, base, fixed, sex, age, weights, SINGLE_YEAR_AGE_BANDS - ) + with np.errstate(over="ignore"): # The adapter's intended overflow. + old, _ = _net_of_di_origin( + prescribed, base, fixed, sex, age, weights, SINGLE_YEAR_AGE_BANDS + ) for row in log.itertuples(index=False): if row.clipped_persons: continue # Intended: iterative capping versus one-pass clipping. @@ -794,6 +874,16 @@ def test_unit_risk_matches_adapter_on_every_unclipped_cell(inputs): # Intended: the adapter leaves weightless rows at factor one, # while the brief's infeasible-cell rules hold for any weight. continue + exact, unreliable = exact_adapter_cell( + base[mask], fixed[mask], prescribed[mask], weights[mask] + ) + np.testing.assert_allclose( + q[mask], exact, rtol=ORACLE_RTOL, atol=1e-15 + ) + if unreliable: + # Intended: the adapter's float64 factor cancels, absorbs the + # free deaths or overflows here; SM1 meets the exact formula. + continue np.testing.assert_allclose(q[mask], old[mask], rtol=RTOL, atol=ATOL) From 3384364a5b8cadce726421041b7a61c45995da24 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sat, 24 Oct 2026 23:38:28 -0400 Subject: [PATCH 16/17] SM1: guard subnormal adapter cells and bound the float tolerance The confirmation of 3078d18b found that the adapter guard missed cells whose sums are all subnormal, where float64 rounds products absolutely, and that at a thousandth the float64 comparison held only empirically. - exact_adapter_cell flags a cell when P - D or E lies within a hundredth of max(P, D), or when max(P, D) is below 2**-1000. At a hundredth, a first-order bound keeps the adapter within RTOL and ATOL for every probability up to one; the design note states it. - The reviewer's three minimized subnormal cells become examples. - np.errstate also ignores invalid, for the 0 * inf the overflow forms. - The design note no longer says every unclipped cell meets the exact formula: weightless infeasible cells are skipped. Correction to 3078d18b's message, which this one supersedes: its "4 failed" of 1,000 fresh CI-sized runs at 5bda0fed was a measurement error by the old runner. That runner caught only AssertionError and stopped at the first ExceptionGroup, after 46 runs. A runner that counts every failure measures 70 of 1,000 at 5bda0fed (seeds 1000-1999), as the reviewer found, and 0 of 3,000 at this commit (seeds 1000-3999). Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 19 ++++++--- tests/mortality_ses/test_properties.py | 56 ++++++++++++++++++++++---- 2 files changed, 63 insertions(+), 12 deletions(-) diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md index 627c95ba..2b3b2790 100644 --- a/docs/design/mortality_ses_sm1.md +++ b/docs/design/mortality_ses_sm1.md @@ -243,12 +243,21 @@ difference: in an infeasible cell whose free support weighs nothing, the adapter leaves the factor at one, while the rules above zero or saturate those rows. A third: the adapter forms its factor `max(P - D, 0) / E` in float64 with no guard on the subtraction or the division. Where `P - D` or -the free deaths `E` lies within a thousandth of `max(P, D)`, its factor -can lose precision, absorb the free deaths or overflow. The differential -test therefore holds every unclipped cell to the adapter's formula +the free deaths `E` lies within a hundredth of `max(P, D)`, its factor +can lose precision, absorb the free deaths or overflow. Where `max(P, D)` +is below 2^-1000, float64 rounds the cell's products absolutely rather +than relatively, and the factor can be off by percents. Elsewhere, for +the at most ten rows the test draws per cell, first-order bounds on the +three sums give the adapter's factor a relative error of at most about +(19/0.01 + 12)u, roughly 2.1e-13 (u = 2^-53). With SM1's own float-path +bound of about 1.3e-12, the two results then differ by at most about +1.55e-12 of the probability, inside the test's +tolerance of 1e-12 plus 1e-12 of the probability for every probability +up to one. The differential test therefore holds every unclipped cell, +apart from the weightless infeasible ones above, to the adapter's formula evaluated exactly, at the oracle's tolerance, and compares the adapter's -float64 result only outside those cells. No shared-module change is -required. +float64 result only outside the unreliable cells. No shared-module change +is required. At the scale of a projection frame (1,000,000 invented rows in 242 cells), `evaluate` takes about 1.4 seconds on the build host, and the diff --git a/tests/mortality_ses/test_properties.py b/tests/mortality_ses/test_properties.py index b1ecd11c..ce2470b0 100644 --- a/tests/mortality_ses/test_properties.py +++ b/tests/mortality_ses/test_properties.py @@ -775,9 +775,12 @@ def exact_adapter_cell(base, fixed, prescribed, weights): baseline deaths, ``D`` its fixed rows' deaths and ``E`` its free rows' baseline deaths; with no fixed row or no free deaths the factor is one. Also returns whether the adapter's float64 factor is unreliable: when - ``P - D`` or ``E`` is within a thousandth of ``max(P, D)``, the - subtraction can cancel, absorb ``E`` or overflow the division. Exact - sums keep an ``E`` that would underflow in float64. + ``P - D`` or ``E`` is within a hundredth of ``max(P, D)``, the + subtraction can cancel, absorb ``E`` or overflow the division; when + ``max(P, D)`` is below ``2**-1000``, float64 rounds its products + absolutely, not relatively. Elsewhere a first-order error bound keeps + the adapter within ``RTOL`` and ``ATOL`` (the design note's intended + differences). Exact sums keep an ``E`` that would underflow in float64. """ weight = [Fraction(w) for w in weights] population = sum( @@ -791,9 +794,13 @@ def exact_adapter_cell(base, fixed, prescribed, weights): ) factor, unreliable = Fraction(1), False if fixed.any(): - scale = max(population, fixed_deaths) / 1000 + largest = max(population, fixed_deaths) gap = abs(population - fixed_deaths) - unreliable = 0 < free_deaths <= scale or 0 < gap <= scale + unreliable = ( + 0 < free_deaths <= largest / 100 + or 0 < gap <= largest / 100 + or 0 < largest < Fraction(2) ** -1000 + ) if free_deaths > 0: factor = max(population - fixed_deaths, 0) / free_deaths q = [ @@ -803,8 +810,41 @@ def exact_adapter_cell(base, fixed, prescribed, weights): return np.asarray(q, dtype=float), unreliable +TINY = 5e-324 # The smallest subnormal float64. + + @PROPERTY @given(invented_adapter_cells()) +@example( # Subnormal sums round absolutely: the factor errs by percents. + ( + np.array([12 * TINY, 0.5, 41 * TINY]), + np.array([True, False, False]), + np.array([0.0, 0.5, 41 * TINY]), + np.array([1.0, 0.0, 0.3]), + np.zeros(3, dtype=int), + np.full(3, 50), + ) +) +@example( # The same with single subnormal rows. + ( + np.array([TINY, 0.5, TINY]), + np.array([True, False, False]), + np.array([0.0, 0.5, TINY]), + np.array([0.3, 0.0, 1.7]), + np.zeros(3, dtype=int), + np.full(3, 50), + ) +) +@example( # Equal subnormal P and D, which leave no gap to flag. + ( + np.array([0.0, 0.0, 0.5, 2 * TINY]), + np.array([True, True, False, False]), + np.array([TINY, TINY, 0.5, 2 * TINY]), + np.array([0.3, 0.3, 0.0, 0.3]), + np.zeros(4, dtype=int), + np.full(4, 50), + ) +) @example( # The adapter's factor overflows on a subnormal free row. ( np.array([2.225073858507e-311, 0.5]), @@ -859,7 +899,8 @@ def test_unit_risk_matches_adapter_on_every_unclipped_cell(inputs): q, log = normalize_probabilities( base, np.ones(len(base)), fixed, prescribed, cells, weights ) - with np.errstate(over="ignore"): # The adapter's intended overflow. + # The adapter's intended overflow, and the 0 * inf it can then form. + with np.errstate(over="ignore", invalid="ignore"): old, _ = _net_of_di_origin( prescribed, base, fixed, sex, age, weights, SINGLE_YEAR_AGE_BANDS ) @@ -882,7 +923,8 @@ def test_unit_risk_matches_adapter_on_every_unclipped_cell(inputs): ) if unreliable: # Intended: the adapter's float64 factor cancels, absorbs the - # free deaths or overflows here; SM1 meets the exact formula. + # free deaths, overflows or rounds subnormal sums here; SM1 + # meets the exact formula. continue np.testing.assert_allclose(q[mask], old[mask], rtol=RTOL, atol=ATOL) From cd12003ace3ecd82d3a32b990bd241e35910eb09 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Sun, 25 Oct 2026 00:08:55 -0400 Subject: [PATCH 17/17] SM1: qualify the subnormal clause and state the zero-gap bound Wording only, from the confirmation of 3384364a: below 2**-1000 float64 can round a cell's products absolutely, but only those below 2**-1022 (subnormal) do; and where P = D exactly the adapter's factor is at most about 19u/0.01, under the test's absolute tolerance. Co-Authored-By: Claude Opus 5.5 --- docs/design/mortality_ses_sm1.md | 28 ++++++++++++++------------ tests/mortality_ses/test_properties.py | 9 +++++---- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/docs/design/mortality_ses_sm1.md b/docs/design/mortality_ses_sm1.md index 2b3b2790..d509a948 100644 --- a/docs/design/mortality_ses_sm1.md +++ b/docs/design/mortality_ses_sm1.md @@ -245,19 +245,21 @@ those rows. A third: the adapter forms its factor `max(P - D, 0) / E` in float64 with no guard on the subtraction or the division. Where `P - D` or the free deaths `E` lies within a hundredth of `max(P, D)`, its factor can lose precision, absorb the free deaths or overflow. Where `max(P, D)` -is below 2^-1000, float64 rounds the cell's products absolutely rather -than relatively, and the factor can be off by percents. Elsewhere, for -the at most ten rows the test draws per cell, first-order bounds on the -three sums give the adapter's factor a relative error of at most about -(19/0.01 + 12)u, roughly 2.1e-13 (u = 2^-53). With SM1's own float-path -bound of about 1.3e-12, the two results then differ by at most about -1.55e-12 of the probability, inside the test's -tolerance of 1e-12 plus 1e-12 of the probability for every probability -up to one. The differential test therefore holds every unclipped cell, -apart from the weightless infeasible ones above, to the adapter's formula -evaluated exactly, at the oracle's tolerance, and compares the adapter's -float64 result only outside the unreliable cells. No shared-module change -is required. +is below 2^-1000, float64 can round the cell's products absolutely rather +than relatively (a product below 2^-1022 is subnormal), and the factor +can be off by percents. Elsewhere, for the at most ten rows the test +draws per cell, first-order bounds on the three sums give the adapter's +factor a relative error of at most about (19/0.01 + 12)u, roughly 2.1e-13 +(u = 2^-53). With SM1's own float-path bound of about 1.3e-12, the two +results then differ by at most about 1.55e-12 of the probability, inside +the test's tolerance of 1e-12 plus 1e-12 of the probability for every +probability up to one. Where `P = D` exactly, the exact factor is zero and +the adapter's is at most about 19u/0.01, roughly 2.1e-13, under the 1e-12 +absolute tolerance. The differential test therefore holds every unclipped +cell, apart from the weightless infeasible ones above, to the adapter's +formula evaluated exactly, at the oracle's tolerance, and compares the +adapter's float64 result only outside the unreliable cells. No +shared-module change is required. At the scale of a projection frame (1,000,000 invented rows in 242 cells), `evaluate` takes about 1.4 seconds on the build host, and the diff --git a/tests/mortality_ses/test_properties.py b/tests/mortality_ses/test_properties.py index ce2470b0..26daa4ba 100644 --- a/tests/mortality_ses/test_properties.py +++ b/tests/mortality_ses/test_properties.py @@ -777,10 +777,11 @@ def exact_adapter_cell(base, fixed, prescribed, weights): Also returns whether the adapter's float64 factor is unreliable: when ``P - D`` or ``E`` is within a hundredth of ``max(P, D)``, the subtraction can cancel, absorb ``E`` or overflow the division; when - ``max(P, D)`` is below ``2**-1000``, float64 rounds its products - absolutely, not relatively. Elsewhere a first-order error bound keeps - the adapter within ``RTOL`` and ``ATOL`` (the design note's intended - differences). Exact sums keep an ``E`` that would underflow in float64. + ``max(P, D)`` is below ``2**-1000``, float64 can round its products + absolutely, not relatively (below ``2**-1022`` a product is subnormal). + Elsewhere a first-order error bound keeps the adapter within ``RTOL`` + and ``ATOL`` (the design note's intended differences). Exact sums keep + an ``E`` that would underflow in float64. """ weight = [Fraction(w) for w in weights] population = sum(