diff --git a/SYSTEM_REQUIREMENTS.md b/SYSTEM_REQUIREMENTS.md index 4424d09..937bc59 100644 --- a/SYSTEM_REQUIREMENTS.md +++ b/SYSTEM_REQUIREMENTS.md @@ -43,8 +43,8 @@ A US build is two heavy phases over a sampling Frame, plus loaders: which is what the imputation benchmark below measures directly. 2. **Calibration** (`populace.calibrate`) — compile targets (national + county control totals) into a sparse constraint matrix over the pool's weight - vector, then optimize log-weights with torch Adam (the bounded relative-error - loss), optionally with L0 generate-big-then-prune. See + vector, then optimize log-weights with torch Adam (capped weighted MAPE), + optionally with L0 generate-big-then-prune. See `packages/populace-calibrate/src/populace/calibrate/solve.py`. 3. **Loaders** (`populace.data`) — pull a published population artifact from the Hugging Face Hub and return it as a policyengine engine dataset. diff --git a/packages/populace-build/README.md b/packages/populace-build/README.md index c5a2a47..78ca155 100644 --- a/packages/populace-build/README.md +++ b/packages/populace-build/README.md @@ -26,9 +26,9 @@ names its donor survey and fails loudly — no silent fallbacks), and the - **rotated holdout** — deterministic target folds so *every* target is held out exactly once across rotations, instead of one lucky split. -All gate losses use the calibrator's bounded relative-error loss -`mean(((est − target)/(target + 1))²)` — scorers consume the same functions, -so there is no calibrator-vs-scorer objective mismatch. +All gate losses use the calibrator's capped weighted-MAPE helper +`weighted_mean(min(abs((estimate − target) / scale), cap))` — scorers consume +the same functions, so there is no calibrator-vs-scorer objective mismatch. The `us` extra adds the rules engine for formula/export checks. Country source loaders are not Python dependencies: source stages are declared in packaged diff --git a/packages/populace-build/src/populace/build/us/puf_aggregate_records.py b/packages/populace-build/src/populace/build/us/puf_aggregate_records.py index 120948c..c57dc4d 100644 --- a/packages/populace-build/src/populace/build/us/puf_aggregate_records.py +++ b/packages/populace-build/src/populace/build/us/puf_aggregate_records.py @@ -380,7 +380,9 @@ def _source_loss( if not errors: return { "loss": 0.0, - "loss_formula": "mean(((estimate - target) / (target + 1)) ** 2)", + "loss_formula": ( + "mean(min(abs((estimate - target) / max(abs(target), 1)), 10))" + ), "n_columns": 0, "within_10pct": 1.0, "max_abs_relative_error": 0.0, @@ -401,7 +403,9 @@ def _source_loss( np.asarray([error["target_total"] for error in errors]), ) ), - "loss_formula": "mean(((estimate - target) / (target + 1)) ** 2)", + "loss_formula": ( + "mean(min(abs((estimate - target) / max(abs(target), 1)), 10))" + ), "n_columns": int(len(errors)), "within_10pct": _json_float( float(np.mean([abs(error["relative_error"]) <= 0.10 for error in errors])) diff --git a/packages/populace-build/tests/test_gates.py b/packages/populace-build/tests/test_gates.py index 2f89b65..d33dac9 100644 --- a/packages/populace-build/tests/test_gates.py +++ b/packages/populace-build/tests/test_gates.py @@ -210,14 +210,14 @@ class TestRelativeErrorLoss: def test_matches_the_calibrator_formula(self) -> None: est = np.asarray([110.0, 90.0]) tgt = np.asarray([100.0, 100.0]) - expected = (((est - tgt) / (tgt + 1.0)) ** 2).mean() + expected = np.abs((est - tgt) / np.maximum(np.abs(tgt), 1.0)).mean() assert relative_error_loss(est, tgt) == pytest.approx(expected) def test_accepts_target_loss_weights(self) -> None: est = np.asarray([110.0, 90.0]) tgt = np.asarray([100.0, 100.0]) weights = np.asarray([10.0, 1.0]) - residual = ((est - tgt) / (tgt + 1.0)) ** 2 + residual = np.abs((est - tgt) / np.maximum(np.abs(tgt), 1.0)) expected = np.average(residual, weights=weights) assert relative_error_loss( @@ -226,6 +226,19 @@ def test_accepts_target_loss_weights(self) -> None: target_loss_weights=weights, ) == pytest.approx(expected) + def test_accepts_target_loss_scales_and_caps_each_row(self) -> None: + est = np.asarray([1_000.0, 50.0]) + tgt = np.asarray([0.0, 100.0]) + scales = np.asarray([100.0, 100.0]) + expected = np.asarray([10.0, 0.5]).mean() + + assert relative_error_loss( + est, + tgt, + target_loss_scales=scales, + target_loss_cap=10.0, + ) == pytest.approx(expected) + def test_shape_mismatch_is_refused(self) -> None: with pytest.raises(ValueError, match="must align"): relative_error_loss(np.zeros(2), np.zeros(3)) diff --git a/packages/populace-build/tests/test_us_fiscal_refresh_builder.py b/packages/populace-build/tests/test_us_fiscal_refresh_builder.py index 84d21de..d91d765 100644 --- a/packages/populace-build/tests/test_us_fiscal_refresh_builder.py +++ b/packages/populace-build/tests/test_us_fiscal_refresh_builder.py @@ -403,18 +403,19 @@ def test_release_gate_failures_reject_missing_critical_targets() -> None: ] -def test_target_value_loss_weights_prioritize_large_targets() -> None: +def test_fiscal_target_loss_weights_prioritize_national_totals() -> None: builder = _load_builder_module() registry = TargetRegistry( ( TargetSpec( - name="small_count", + name="national_income_tax_total", entity="household", value=10.0, source="fixture", + metadata={"target_role": "federal_income_tax_total"}, ), TargetSpec( - name="large_amount", + name="distribution_row", entity="household", value=1_000_000.0, source="fixture", @@ -423,33 +424,26 @@ def test_target_value_loss_weights_prioritize_large_targets() -> None: country="us", ) - weights = builder._target_value_loss_weights(registry) + weights = builder._fiscal_target_loss_weights(registry) assert weights.shape == (2,) assert weights.mean() == 1.0 - assert weights[1] > weights[0] * 10_000 + assert weights[0] == weights[1] * builder.US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER -def test_target_value_loss_weights_boost_ctc_total_role() -> None: +def test_fiscal_target_loss_weights_downweight_state_rows() -> None: builder = _load_builder_module() registry = TargetRegistry( ( TargetSpec( - name="ctc_total", + name="state_role_row", entity="tax_unit", value=100.0, source="fixture", - metadata={"target_role": "ctc_total"}, + metadata={"state_fips": "06", "target_role": "tanf_total"}, ), TargetSpec( - name="other_same_value", - entity="tax_unit", - value=100.0, - source="fixture", - metadata={"target_role": "eitc_total"}, - ), - TargetSpec( - name="legacy_named_ctc_without_role", + name="national_row", entity="tax_unit", value=100.0, source="fixture", @@ -458,11 +452,10 @@ def test_target_value_loss_weights_boost_ctc_total_role() -> None: country="us", ) - weights = builder._target_value_loss_weights(registry) + weights = builder._fiscal_target_loss_weights(registry) assert weights.mean() == 1.0 - assert weights[0] == weights[1] * builder.US_CTC_TARGET_LOSS_WEIGHT_MULTIPLIER - assert weights[2] == weights[1] + assert weights[0] == weights[1] * builder.US_STATE_TARGET_LOSS_MULTIPLIER def test_release_gate_failures_reject_bad_ctc_fit() -> None: diff --git a/packages/populace-build/tests/test_us_puf_aggregate_records.py b/packages/populace-build/tests/test_us_puf_aggregate_records.py index aa001f6..f16c70c 100644 --- a/packages/populace-build/tests/test_us_puf_aggregate_records.py +++ b/packages/populace-build/tests/test_us_puf_aggregate_records.py @@ -346,7 +346,9 @@ def test_audit_reports_source_reconstruction_recovery( old_loss = audit["source_reconstruction_loss"]["old_drop_aggregate"] disaggregated_loss = audit["source_reconstruction_loss"]["disaggregated"] assert old_loss["loss"] > 0 - assert old_loss["loss_formula"] == "mean(((estimate - target) / (target + 1)) ** 2)" + assert old_loss["loss_formula"] == ( + "mean(min(abs((estimate - target) / max(abs(target), 1)), 10))" + ) subset_audit = audit_puf_aggregate_disaggregation( mini_puf, seed=42, @@ -365,7 +367,7 @@ def test_audit_reports_source_reconstruction_recovery( ) ) assert old_loss["within_10pct"] < 1.0 - assert disaggregated_loss["loss"] < 1e-18 + assert disaggregated_loss["loss"] < 1e-12 assert disaggregated_loss["within_10pct"] == pytest.approx(1.0) agi = next( diff --git a/packages/populace-calibrate/README.md b/packages/populace-calibrate/README.md index 072e5cb..08ce278 100644 --- a/packages/populace-calibrate/README.md +++ b/packages/populace-calibrate/README.md @@ -21,10 +21,11 @@ reproduces them. Uncompilable targets (missing column, zero `mean` denominator) are **skipped and reported**, never dropped silently. 3. **Solve for calibrated weights.** `calibrate(frame, targets, ...)` optimizes - the log-weights with torch Adam to minimize the **bounded relative-error - loss** `mean(((A @ w - b)/(b + 1))**2)`. Weights stay strictly positive by - construction (`w = exp(log_w)`). The - result carries a new `Frame` with `CALIBRATED` weights, per-target + the log-weights with torch Adam to minimize **capped weighted MAPE**: + `weighted_mean(min(abs((A @ w - b) / scale), cap))`. By default + `scale = max(abs(target), abs(initial_estimate), 1)` and `cap = 10` + (1000%). Weights stay strictly positive by construction (`w = exp(log_w)`). + The result carries a new `Frame` with `CALIBRATED` weights, per-target diagnostics, and the loss trajectory. ## Load-bearing options diff --git a/packages/populace-calibrate/pyproject.toml b/packages/populace-calibrate/pyproject.toml index 249d466..e27ccbf 100644 --- a/packages/populace-calibrate/pyproject.toml +++ b/packages/populace-calibrate/pyproject.toml @@ -1,7 +1,7 @@ [project] name = "populace-calibrate" version = "0.1.0" -description = "The populace representation operator: compile targets into a sparse constraint system over a Frame and solve for calibrated weights (bounded relative-error loss, torch on log-weights, hard max-weight-ratio guard, optional L0 hard-concrete pruning)" +description = "The populace representation operator: compile targets into a sparse constraint system over a Frame and solve for calibrated weights (capped weighted-MAPE loss, torch on log-weights, hard max-weight-ratio guard, optional L0 hard-concrete pruning)" readme = "README.md" requires-python = ">=3.13" dependencies = [ diff --git a/packages/populace-calibrate/src/populace/calibrate/__init__.py b/packages/populace-calibrate/src/populace/calibrate/__init__.py index 0d36917..049da8c 100644 --- a/packages/populace-calibrate/src/populace/calibrate/__init__.py +++ b/packages/populace-calibrate/src/populace/calibrate/__init__.py @@ -4,9 +4,9 @@ produced. Compiles declared facts — population control totals, counts, averages with standard-error-style tolerances — into a sparse linear constraint system over a :class:`~populace.frame.Frame`, then solves for the weight vector that -best reproduces them under the bounded relative-error loss -``mean(((A @ w - b)/(b + 1))**2)``, optimized with torch's Adam over the -log-weights (positivity by construction). Multi-period targets stack as +best reproduces them under capped weighted MAPE, +``weighted_mean(min(abs((A @ w - b) / scale), cap))``, optimized with torch's +Adam over the log-weights (positivity by construction). Multi-period targets stack as ``(target, period)`` rows over the *same* weight vector — the charter's "one weight per trajectory". @@ -95,6 +95,7 @@ def _assert_frame_compatible(version: str, required: tuple[int, int]) -> None: CalibrationResult, TargetDiagnostic, calibrate, + default_target_loss_scales, relative_error_loss, ) from populace.calibrate.target import ( # noqa: E402 - after the compat gate @@ -120,6 +121,7 @@ def _assert_frame_compatible(version: str, required: tuple[int, int]) -> None: "TargetSpec", "build_constraint_matrix", "calibrate", + "default_target_loss_scales", "diagnostics_payload", "relative_error_loss", "specs_from_pe_surface", diff --git a/packages/populace-calibrate/src/populace/calibrate/matrix.py b/packages/populace-calibrate/src/populace/calibrate/matrix.py index d4006d1..5f2fb55 100644 --- a/packages/populace-calibrate/src/populace/calibrate/matrix.py +++ b/packages/populace-calibrate/src/populace/calibrate/matrix.py @@ -30,15 +30,6 @@ __all__ = ["CalibrationProblem", "SkippedTarget", "build_constraint_matrix"] -#: Half-width of the forbidden band around ``-1`` for a compiled target value. -#: The bounded relative-error loss divides each residual by -#: ``(target_value + 1)``; a compiled -#: value within ``_DENOM_EPS`` of ``-1`` drives that denominator to ~0, so the -#: loss, its gradients, and every weight go NaN — surfacing only as the kernel's -#: opaque "Weights must be finite" error. We reject such targets at compile time -#: instead, naming the culprit and the cause. -_DENOM_EPS = 1e-8 - @dataclass(frozen=True) class SkippedTarget: @@ -314,30 +305,10 @@ def build_constraint_matrix( f"({len(skipped)} skipped): {detail}." ) - target_vector = np.asarray(values, dtype=np.float64) - # Guard the loss denominator: the bounded relative-error loss divides each - # residual by ``(target_value + 1)``. A compiled value at -1 (a raw - # ``value=-1``, or a ``mean`` whose ``value`` is exactly 1 below the current - # mean, since its compiled RHS is ``value - current_mean``) makes that - # denominator ~0 -> NaN loss -> NaN gradients -> all-NaN weights, which would - # otherwise surface only as the kernel's opaque "Weights must be finite". - near_minus_one = np.abs(target_vector + 1.0) < _DENOM_EPS - if near_minus_one.any(): - culprits = "; ".join( - f"{names[i]} (compiled target value {target_vector[i]:.6g})" - for i in np.flatnonzero(near_minus_one) - ) - raise ValueError( - "Target(s) compile to a value at -1, which zeroes the " - "relative-error loss denominator (target_value + 1) and makes every " - f"weight NaN: {culprits}. Shift the target value away from -1 (for a " - "'mean', away from exactly 1 below the current mean)." - ) - matrix = sparse.csr_array(np.vstack(rows)) return CalibrationProblem( matrix=matrix, - target_vector=target_vector, + target_vector=np.asarray(values, dtype=np.float64), names=tuple(names), initial_weights=initial, weight_entity=weight_entity, diff --git a/packages/populace-calibrate/src/populace/calibrate/solve.py b/packages/populace-calibrate/src/populace/calibrate/solve.py index 36f9fe1..0c8581e 100644 --- a/packages/populace-calibrate/src/populace/calibrate/solve.py +++ b/packages/populace-calibrate/src/populace/calibrate/solve.py @@ -3,9 +3,14 @@ :func:`calibrate` is the representation operator. It compiles a :class:`~populace.calibrate.target.TargetSet` against a :class:`~populace.frame.Frame` (:mod:`populace.calibrate.matrix`), then optimizes -the weight vector of ``weight_entity`` to minimize the bounded relative-error loss +the weight vector of ``weight_entity`` to minimize fixed-scale weighted MAPE - ``mean(((A @ w - b) / (b + 1))**2)`` + ``weighted_mean(abs((A @ w - b) / s))`` + +where ``s`` is a per-target scale fixed before optimization. By default, +``s = max(abs(b), abs(A @ w0), 1)``: zero-valued or tiny targets can no longer +dominate the objective merely because their denominator is near zero, while +large national aggregates remain measured as proportional misses. with torch's Adam over the **log-weights** (so weights stay strictly positive by construction). It returns a @@ -55,6 +60,7 @@ __all__ = [ "calibrate", + "default_target_loss_scales", "relative_error_loss", "CalibrationResult", "TargetDiagnostic", @@ -86,6 +92,10 @@ #: the ``log10`` bracket by 2^10, far finer than the count is resolvable. _DEFAULT_BUDGET_ITERS = 10 +# Per-target contribution cap for weighted MAPE. A target can contribute at most +# a 1000% scaled miss to the objective. +_DEFAULT_TARGET_LOSS_CAP = 10.0 + @dataclass(frozen=True) class TargetDiagnostic: @@ -140,7 +150,7 @@ class CalibrationResult: search settled on, not the value passed in. n_nonzero: Number of calibrated weights above the prune threshold. With a ``target_records`` budget, the quantity the search drives toward it. - closing_loss: The bounded relative-error loss evaluated once on the + closing_loss: The capped weighted-MAPE calibration loss evaluated once on the *returned* weights (after the closing mass/cap projections). Exposed as :attr:`final_loss`; recorded separately from the trajectory, whose tail is a pre-step/pre-projection value. @@ -170,12 +180,12 @@ class CalibrationResult: @property def initial_loss(self) -> float: - """The bounded relative-error loss under the input weights.""" + """The capped weighted-MAPE calibration loss under the input weights.""" return float(self.loss_trajectory[0]) @property def final_loss(self) -> float: - """The bounded relative-error loss of the *returned* weights. + """The capped weighted-MAPE calibration loss of the *returned* weights. This is a single eval-mode evaluation on the weights actually returned — after the closing mass/cap projections — so it describes the calibrated @@ -242,17 +252,23 @@ def relative_error_loss( targets: np.ndarray, *, target_loss_weights: np.ndarray | None = None, + target_loss_scales: np.ndarray | None = None, + target_loss_cap: float = _DEFAULT_TARGET_LOSS_CAP, ) -> float: - """THE loss, in numpy: weighted ``((est - tgt)/(tgt + 1))**2``. + """THE loss, in numpy: capped weighted MAPE on fixed row scales. The single canonical definition every measurement imports — the solver's closing loss, the acceptance gates, and scorers all call this function (the torch twin below is the autograd path of the same formula). Refuses non-finite inputs: a NaN estimate is a harness bug, not a large miss. - When ``target_loss_weights`` is omitted, this is the historical unweighted - mean. When supplied, weights must align to targets and are normalized by - their own sum, so multiplying all weights by a constant does not change the - objective. + + ``target_loss_scales`` is the row denominator ``s`` in + ``abs((estimate - target) / s)``. If omitted, the diagnostic helper uses + ``max(abs(target), 1)``; :func:`calibrate` supplies the stronger production + default ``max(abs(target), abs(initial_estimate), 1)`` once it has compiled + the target matrix. Target weights are normalized by their own sum, so + multiplying all weights by a constant does not change the objective. Each + row's scaled absolute miss is capped by ``target_loss_cap``. """ estimates = np.asarray(estimates, dtype=np.float64) targets = np.asarray(targets, dtype=np.float64) @@ -266,8 +282,13 @@ def relative_error_loss( "relative_error_loss requires finite inputs; got non-finite " "estimate or target values." ) - rel = (estimates - targets) / (targets + 1.0) - loss = rel**2 + scales = _validate_target_loss_scales( + target_loss_scales, + targets.shape, + targets=targets, + ) + cap = _validate_target_loss_cap(target_loss_cap) + loss = np.minimum(np.abs((estimates - targets) / scales), cap) weights = _validate_target_loss_weights(target_loss_weights, targets.shape) if weights is None: return float(loss.mean()) @@ -278,26 +299,61 @@ def _relative_error_loss( estimate: torch.Tensor, targets: torch.Tensor, target_loss_weights: torch.Tensor | None, + target_loss_scales: torch.Tensor, + target_loss_cap: float, ) -> torch.Tensor: - """The relative-error loss, optionally averaged with target row weights. - - The ``+1`` in the *denominator* is the regularizer: it keeps the loss finite - and well-scaled for targets near zero (a zero-valued count target then - contributes ``est**2`` rather than dividing by zero). The numerator is the - raw residual ``est - tgt``, so the loss is minimized exactly at ``est = tgt``. - - Some older reweighting formulas carried a ``+1`` in the numerator too. That - biases the optimum to ``est = tgt - 1`` and is fatal for small-valued - targets (a count of 5 converges to 4). We use the raw residual — this is - also the loss this docstring has always described. - """ - rel_error = (estimate - targets) / (targets + 1.0) - loss = rel_error**2 + """The capped weighted-MAPE target loss, optionally averaged with row weights.""" + scaled_error = (estimate - targets) / target_loss_scales + loss = torch.clamp( + torch.abs(scaled_error), + max=_validate_target_loss_cap(target_loss_cap), + ) if target_loss_weights is None: return loss.mean() return (loss * target_loss_weights).sum() / target_loss_weights.sum() +def _validate_target_loss_cap(target_loss_cap: float) -> float: + cap = float(target_loss_cap) + if not np.isfinite(cap) or cap <= 0.0: + raise ValueError( + f"target_loss_cap must be positive and finite, got {target_loss_cap!r}." + ) + return cap + + +def default_target_loss_scales( + targets: np.ndarray, + initial_estimates: np.ndarray, +) -> np.ndarray: + """Default fixed row scales for calibration. + + The old objective used ``target + 1`` as the denominator. That makes a + zero-valued target with a large starting estimate dominate the loss by many + orders of magnitude. The production scale is instead the largest meaningful + size already known before optimization: the absolute target, the absolute + starting estimate, or one unit. + """ + targets = np.asarray(targets, dtype=np.float64) + initial_estimates = np.asarray(initial_estimates, dtype=np.float64) + if targets.shape != initial_estimates.shape: + raise ValueError( + "targets and initial_estimates must align, got shapes " + f"{targets.shape} vs {initial_estimates.shape}." + ) + if not (np.isfinite(targets).all() and np.isfinite(initial_estimates).all()): + raise ValueError( + "default_target_loss_scales requires finite targets and estimates." + ) + return np.maximum.reduce( + [ + np.abs(targets), + np.abs(initial_estimates), + np.ones_like(targets, dtype=np.float64), + ] + ) + + def _validate_target_loss_weights( target_loss_weights: np.ndarray | None, shape: tuple[int, ...], @@ -319,6 +375,28 @@ def _validate_target_loss_weights( return weights +def _validate_target_loss_scales( + target_loss_scales: np.ndarray | None, + shape: tuple[int, ...], + *, + targets: np.ndarray, +) -> np.ndarray: + if target_loss_scales is None: + scales = np.maximum(np.abs(np.asarray(targets, dtype=np.float64)), 1.0) + else: + scales = np.asarray(target_loss_scales, dtype=np.float64) + if scales.shape != shape: + raise ValueError( + "target_loss_scales must align with targets, got shapes " + f"{scales.shape} vs {shape}." + ) + if not np.isfinite(scales).all(): + raise ValueError("target_loss_scales must be finite.") + if (scales <= 0).any(): + raise ValueError("target_loss_scales must be positive.") + return scales + + def _target_loss_weight_options( target_loss_weights: np.ndarray | None, ) -> Mapping[str, object]: @@ -334,6 +412,24 @@ def _target_loss_weight_options( } +def _target_loss_scale_options( + target_loss_scales: np.ndarray, + *, + kind: str, + target_loss_cap: float, +) -> Mapping[str, object]: + scales = np.asarray(target_loss_scales, dtype=np.float64) + return { + "kind": kind, + "formula": "weighted_mean(min(abs((estimate - target) / scale), cap))", + "cap": float(target_loss_cap), + "n": int(scales.shape[0]), + "min": float(scales.min()), + "median": float(np.median(scales)), + "max": float(scales.max()), + } + + def _build_diagnostics( problem: CalibrationProblem, frame: Frame, @@ -398,6 +494,8 @@ def _optimize( matrix: torch.Tensor, targets: torch.Tensor, target_loss_weights: torch.Tensor | None, + target_loss_scales: torch.Tensor, + target_loss_cap: float, initial_weights: np.ndarray, *, epochs: int, @@ -411,7 +509,7 @@ def _optimize( ) -> tuple[np.ndarray, np.ndarray]: """Run the torch optimization and return ``(final_weights, loss_trajectory)``. - Optimizes the log-weights with Adam against the bounded relative-error loss. + Optimizes the log-weights with Adam against capped weighted MAPE. Positivity is by construction (``w = exp(log_w)`` times optional gates). The hard constraints — mass conservation and ``max_weight_ratio`` — are applied by projecting the realized weights after each step, so they hold on the @@ -444,7 +542,13 @@ def _optimize( if gates is not None: weights = weights * gates() estimate = _apply_constraint(matrix, weights) - loss = _relative_error_loss(estimate, targets, target_loss_weights) + loss = _relative_error_loss( + estimate, + targets, + target_loss_weights, + target_loss_scales, + target_loss_cap, + ) penalty = ( l0_lambda * gates.get_penalty() if (gates is not None and l0_lambda > 0.0) @@ -509,6 +613,8 @@ def _search_l0_lambda_for_budget( matrix: torch.Tensor, targets: torch.Tensor, target_loss_weights: torch.Tensor | None, + target_loss_scales: torch.Tensor, + target_loss_cap: float, initial_weights: np.ndarray, *, target_records: int, @@ -567,6 +673,8 @@ def evaluate(lam: float) -> tuple[np.ndarray, np.ndarray, int]: matrix, targets, target_loss_weights, + target_loss_scales, + target_loss_cap, initial_weights, epochs=epochs, learning_rate=learning_rate, @@ -726,7 +834,7 @@ def calibrate( weight_entity: str = "household", method: str = "apg", epochs: int = 256, - learning_rate: float = 0.1, + learning_rate: float = 0.02, mass: str = FREE_MASS, max_weight_ratio: float | None = None, target_records: int | None = None, @@ -736,12 +844,14 @@ def calibrate( budget_iters: int = _DEFAULT_BUDGET_ITERS, seed: int = 0, target_loss_weights: np.ndarray | None = None, + target_loss_scales: np.ndarray | None = None, + target_loss_cap: float = _DEFAULT_TARGET_LOSS_CAP, ) -> CalibrationResult: """Calibrate ``weight_entity``'s weights to ``targets`` over ``frame``. Compiles the targets into a sparse system and optimizes the log-weights with - Adam to minimize the bounded relative-error loss - ``mean(((A @ w - b)/(b + 1))**2)``. Returns a new frame whose + Adam to minimize capped fixed-scale weighted MAPE + ``weighted_mean(min(abs((A @ w - b) / s), cap))``. Returns a new frame whose ``weight_entity`` weights are :class:`~populace.frame.WeightKind.CALIBRATED`. Args: @@ -755,9 +865,10 @@ def calibrate( first-order method here; the label is carried for the manifest and future solver swaps. Any other value is rejected. epochs: Number of optimization steps. - learning_rate: Adam learning rate on the log-weights. The reference - ``reweight`` is sensitive to this — it does well above ~0.1 and - breaks down near ~0.005 — so the default matches that regime. + learning_rate: Adam learning rate on the log-weights. Capped MAPE has a + nearly constant gradient away from zero, so the default is lower + than the old squared-error objective to avoid oscillating around + feasible targets. mass: :data:`FREE_MASS` (default) to let the total move, or :data:`CONSERVE_MASS` to hold it to the input total. max_weight_ratio: If given, a hard per-record cap: no calibrated weight @@ -783,10 +894,16 @@ def calibrate( seed: Seed for torch's RNG (the gate sampling), for reproducibility. target_loss_weights: Optional non-negative row weights aligned to the supplied :class:`TargetSet`. When omitted, every compiled target row - contributes equally (historical behavior). When supplied, the weights - for skipped targets are dropped with those targets, and the squared - bounded relative errors for compiled rows are averaged with the - remaining weights, normalized by their sum. + contributes equally. When supplied, the weights for skipped targets + are dropped with those targets, and the capped scaled misses for + compiled rows are averaged with the remaining weights, normalized by + their sum. + target_loss_scales: Optional positive row scales aligned to the supplied + :class:`TargetSet`. When omitted, compiled rows use + :func:`default_target_loss_scales`, fixed from the input weights. + Supplying scales is mainly for harnesses and specialized releases. + target_loss_cap: Positive per-row cap on scaled absolute misses. The + default caps each target's objective contribution at 1000%. Returns: A :class:`CalibrationResult` with the calibrated frame, per-target @@ -833,11 +950,21 @@ def calibrate( ) if budget_iters <= 0: raise ValueError(f"budget_iters must be positive, got {budget_iters!r}.") + target_loss_cap = _validate_target_loss_cap(target_loss_cap) target_loss_weights_input = _validate_target_loss_weights( target_loss_weights, (len(targets),), ) + target_loss_scales_input = ( + None + if target_loss_scales is None + else _validate_target_loss_scales( + target_loss_scales, + (len(targets),), + targets=np.asarray([target.value for target in targets], dtype=np.float64), + ) + ) problem = build_constraint_matrix(frame, targets, weight_entity) initial = problem.initial_weights w0 = initial.values @@ -847,6 +974,8 @@ def calibrate( matrix_t = _torch_constraint_matrix(problem.matrix) targets_t = torch.tensor(problem.target_vector, dtype=torch.float32) target_loss_weights_np: np.ndarray | None = None + target_loss_scales_np: np.ndarray + target_loss_scale_kind = "default_initial_or_target" if target_loss_weights_input is not None: weights_by_key = { target.key: weight @@ -860,11 +989,32 @@ def calibrate( target_loss_weights_np, problem.target_vector.shape, ) + if target_loss_scales_input is not None: + scales_by_key = { + target.key: scale + for target, scale in zip(targets, target_loss_scales_input, strict=True) + } + target_loss_scales_np = np.asarray( + [scales_by_key[target.key] for target in problem.targets], + dtype=np.float64, + ) + target_loss_scales_np = _validate_target_loss_scales( + target_loss_scales_np, + problem.target_vector.shape, + targets=problem.target_vector, + ) + target_loss_scale_kind = "provided" + else: + target_loss_scales_np = default_target_loss_scales( + problem.target_vector, + problem.estimates(w0), + ) target_loss_weights_t = ( torch.tensor(target_loss_weights_np, dtype=torch.float32) if target_loss_weights_np is not None else None ) + target_loss_scales_t = torch.tensor(target_loss_scales_np, dtype=torch.float32) if target_records is not None: # Budget control (Finding 3): search l0_lambda so the achieved non-zero @@ -875,6 +1025,8 @@ def calibrate( matrix_t, targets_t, target_loss_weights_t, + target_loss_scales_t, + target_loss_cap, w0, target_records=target_records, epochs=epochs, @@ -895,6 +1047,8 @@ def calibrate( matrix_t, targets_t, target_loss_weights_t, + target_loss_scales_t, + target_loss_cap, w0, epochs=epochs, learning_rate=learning_rate, @@ -924,14 +1078,15 @@ def calibrate( diagnostics = _build_diagnostics(problem, frame, w0, final_weights) # One closing eval-mode loss on the RETURNED weights (Finding 8): the same - # bounded relative-error loss the optimizer minimizes, - # mean(((A@w - b)/(b+1))**2), - # evaluated after the closing mass/cap projections — so final_loss describes + # capped weighted-MAPE loss the optimizer minimizes, evaluated after the closing + # mass/cap projections — so final_loss describes # what calibrate returns, not the trajectory's pre-projection tail. closing_loss = relative_error_loss( problem.estimates(final_weights), problem.target_vector, target_loss_weights=target_loss_weights_np, + target_loss_scales=target_loss_scales_np, + target_loss_cap=target_loss_cap, ) return CalibrationResult( @@ -955,6 +1110,11 @@ def calibrate( "target_records": target_records, "seed": seed, "target_loss_weights": _target_loss_weight_options(target_loss_weights_np), + "target_loss_scales": _target_loss_scale_options( + target_loss_scales_np, + kind=target_loss_scale_kind, + target_loss_cap=target_loss_cap, + ), "matrix_format": ( "sparse_csr" if matrix_t.layout == torch.sparse_csr else "dense" ), @@ -985,7 +1145,7 @@ def _apply_weights( factor = calibrated.total / initial.total if initial.total != 0 else None reason = ( f"calibrated {weight_entity!r} weights to {len(targets)} target(s) " - "(bounded relative-error loss); total mass free to move" + "(capped weighted-MAPE loss); total mass free to move" ) return frame.with_weights( weight_entity, diff --git a/packages/populace-calibrate/src/populace/calibrate/target.py b/packages/populace-calibrate/src/populace/calibrate/target.py index 9f2e959..7760d40 100644 --- a/packages/populace-calibrate/src/populace/calibrate/target.py +++ b/packages/populace-calibrate/src/populace/calibrate/target.py @@ -78,8 +78,7 @@ class Target: The frame is responsible for carrying period-specific columns; ``period`` is metadata the compiler stacks on, not a column lookup. tolerance: Optional absolute tolerance for "is this target hit". Used - by diagnostics, not by the loss (the loss is the bounded relative - error over all targets jointly). + by diagnostics, not by the capped weighted-MAPE loss. filter: Optional column name or callable producing a boolean (or 0/1) per-record mask. For ``count`` it selects which records to count; for ``mean`` it selects the denominator population; for diff --git a/packages/populace-calibrate/tests/test_matrix.py b/packages/populace-calibrate/tests/test_matrix.py index 24bb6e0..0748780 100644 --- a/packages/populace-calibrate/tests/test_matrix.py +++ b/packages/populace-calibrate/tests/test_matrix.py @@ -88,17 +88,10 @@ def test_multi_period_targets_become_distinct_rows(multiperiod_frame) -> None: assert any("2030" in n for n in problem.names) -def test_target_value_near_minus_one_is_rejected_with_a_named_error( +def test_target_value_near_minus_one_compiles_under_scaled_mape( feasible_frame, ) -> None: - """A compiled RHS of exactly -1 zeroes the loss denominator ``(b + 1)``. - - The bounded relative-error loss divides each residual by ``(target + 1)``. - A target value of exactly -1 makes that denominator zero, so the loss is NaN, the gradients are - NaN, every weight goes NaN, and the kernel raises an opaque "Weights must be - finite" error naming neither the target nor the cause. The compiler must - reject the offending target up front, naming it. - """ + """A target at -1 is no longer special: scaled MAPE has no ``b + 1`` denominator.""" frame, _ = feasible_frame(n=50) targets = TargetSet( ( @@ -110,17 +103,13 @@ def test_target_value_near_minus_one_is_rejected_with_a_named_error( ), ) ) - with pytest.raises(ValueError, match="counts_to_neg_one"): - build_constraint_matrix(frame, targets) - + problem = build_constraint_matrix(frame, targets) + assert problem.names == ("counts_to_neg_one@0",) + assert problem.target_vector.tolist() == [-1.0] -def test_mean_target_one_below_current_mean_is_rejected(feasible_frame) -> None: - """A ``mean`` whose compiled RHS lands at -1 is rejected too. - A ``mean`` target compiles to a right-hand side of ``value - current_mean``; - a value exactly 1 below the current mean lands the RHS at -1 — the same - zero-denominator landmine — and must be named, not silently NaN'd. - """ +def test_mean_target_one_below_current_mean_compiles(feasible_frame) -> None: + """Mean rows whose compiled RHS lands at -1 are valid under scaled MAPE.""" frame, _ = feasible_frame(n=50) current_mean = float(frame.table("household")["income"].to_numpy().mean()) targets = TargetSet( @@ -134,8 +123,9 @@ def test_mean_target_one_below_current_mean_is_rejected(feasible_frame) -> None: ), ) ) - with pytest.raises(ValueError, match="mean_one_below"): - build_constraint_matrix(frame, targets) + problem = build_constraint_matrix(frame, targets) + assert problem.names == ("mean_one_below@0",) + assert problem.target_vector.tolist() == pytest.approx([-1.0]) def test_normal_target_value_compiles_fine(feasible_frame) -> None: diff --git a/packages/populace-calibrate/tests/test_registry.py b/packages/populace-calibrate/tests/test_registry.py index db97c3e..ac1fd21 100644 --- a/packages/populace-calibrate/tests/test_registry.py +++ b/packages/populace-calibrate/tests/test_registry.py @@ -169,7 +169,7 @@ def test_registry_drives_a_real_calibration(self, feasible_frame) -> None: country="us", ) result = calibrate(frame, registry.to_target_set(), epochs=300, seed=0) - assert result.final_loss < result.initial_loss * 1e-2 + assert result.final_loss < result.initial_loss * 0.05 for diag in result.diagnostics: assert abs(diag.relative_error) < 0.02 diff --git a/packages/populace-calibrate/tests/test_solve.py b/packages/populace-calibrate/tests/test_solve.py index 632ea6c..c879bed 100644 --- a/packages/populace-calibrate/tests/test_solve.py +++ b/packages/populace-calibrate/tests/test_solve.py @@ -12,7 +12,13 @@ import numpy as np import pytest -from populace.calibrate import Target, TargetSet, calibrate +from populace.calibrate import ( + Target, + TargetSet, + calibrate, + default_target_loss_scales, + relative_error_loss, +) from populace.frame import WeightKind @@ -45,7 +51,7 @@ def test_calibration_reduces_loss_and_hits_feasible_targets(feasible_frame) -> N ) ) result = calibrate(frame, targets, epochs=400, seed=0) - assert result.final_loss < result.initial_loss * 1e-3 + assert result.final_loss < result.initial_loss * 0.01 for diag in result.diagnostics: assert abs(diag.relative_error) < 0.01 # within 1% @@ -160,6 +166,61 @@ def test_target_loss_weights_must_survive_skipped_targets(feasible_frame) -> Non ) +def test_target_loss_scales_follow_skipped_targets(feasible_frame) -> None: + frame, truths = feasible_frame() + targets = TargetSet( + ( + Target( + name="missing_measure", + entity="household", + aggregation="sum", + value=1.0, + measure="missing_measure", + ), + _population_target(truths["population"], 1.0), + ) + ) + + result = calibrate( + frame, + targets, + epochs=50, + seed=0, + target_loss_scales=np.asarray([10.0, 20.0]), + ) + + assert [skip.target.name for skip in result.skipped] == ["missing_measure"] + assert result.options["target_loss_scales"]["kind"] == "provided" + assert result.options["target_loss_scales"]["n"] == 1 + assert result.options["target_loss_scales"]["min"] == pytest.approx(20.0) + assert result.options["target_loss_scales"]["max"] == pytest.approx(20.0) + + +def test_target_loss_scales_must_survive_skipped_targets(feasible_frame) -> None: + frame, truths = feasible_frame() + targets = TargetSet( + ( + Target( + name="missing_measure", + entity="household", + aggregation="sum", + value=1.0, + measure="missing_measure", + ), + _population_target(truths["population"], 1.0), + ) + ) + + with pytest.raises(ValueError, match="target_loss_scales"): + calibrate( + frame, + targets, + epochs=50, + seed=0, + target_loss_scales=np.asarray([10.0, 0.0]), + ) + + def test_max_weight_ratio_is_respected(feasible_frame) -> None: frame, truths = feasible_frame() w0 = frame.resolve_weights("household").values @@ -348,7 +409,7 @@ def test_prune_with_conserve_and_cap_keeps_pruned_records_pruned( _income_target(truths["income"], 1.0), ) ) - free = calibrate(frame, targets, epochs=400, seed=0, l0_lambda=3e-3) + free = calibrate(frame, targets, epochs=400, seed=0, l0_lambda=3e-2) free_nonzero = int((free.frame.resolve_weights("household").values > 1e-6).sum()) assert free_nonzero < 100 # free mass prunes hard @@ -461,12 +522,17 @@ def test_final_loss_describes_the_returned_weights(feasible_frame) -> None: seed=0, mass="conserve", max_weight_ratio=2.0, + learning_rate=0.1, ) - # Recompute the bounded relative-error loss on the returned weights directly. + # Recompute the capped weighted-MAPE loss on the returned weights directly. # final_loss is a float64 closing eval, so it matches to machine epsilon. b = result.problem.target_vector + scales = default_target_loss_scales( + b, + result.problem.estimates(result.initial_weights), + ) est = result.problem.estimates(result.weights) - true_loss = float((((est - b) / (b + 1.0)) ** 2).mean()) + true_loss = relative_error_loss(est, b, target_loss_scales=scales) assert abs(result.final_loss - true_loss) < 1e-9 # final_loss is NOT merely the trajectory tail (which is pre-projection): on # this conserve+cap run they differ materially. @@ -474,7 +540,7 @@ def test_final_loss_describes_the_returned_weights(feasible_frame) -> None: # initial_loss still describes the input weights (the trajectory head), # to float32 precision since the trajectory is computed in float32 torch. est0 = result.problem.estimates(result.initial_weights) - true_initial = float((((est0 - b) / (b + 1.0)) ** 2).mean()) + true_initial = relative_error_loss(est0, b, target_loss_scales=scales) assert abs(result.initial_loss - true_initial) < 1e-5 diff --git a/tools/build_us_fiscal_refresh_release.py b/tools/build_us_fiscal_refresh_release.py index 440c0c1..6aadcbc 100644 --- a/tools/build_us_fiscal_refresh_release.py +++ b/tools/build_us_fiscal_refresh_release.py @@ -68,10 +68,25 @@ CALIBRATION_FILENAME = "populace_us_2024_calibration.npz" POST_EXPORT_ABSOLUTE_TOLERANCE = 1_000_000.0 POST_EXPORT_RELATIVE_TOLERANCE = 5e-4 -US_FISCAL_TARGET_LOSS_WEIGHTING = "absolute_target_value_with_ctc_priority" -US_CTC_TARGET_LOSS_WEIGHT_MULTIPLIER = 10.0 +US_FISCAL_TARGET_LOSS_WEIGHTING = ( + "semantic_weighted_mape_initial_or_target_scale_cap_1000pct" +) +US_FISCAL_TARGET_LOSS_CAP = 10.0 +US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER = 25.0 +US_STATE_TARGET_LOSS_MULTIPLIER = 0.25 US_FISCAL_TARGET_ROLE_LOSS_MULTIPLIERS = { - "ctc_total": US_CTC_TARGET_LOSS_WEIGHT_MULTIPLIER, + "aca_spending": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "ctc_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "eitc_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "federal_income_tax_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "income_tax_before_credits_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "medicare_part_b_premium_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "refundable_ctc_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "snap_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "social_security_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "ssi_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "tanf_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, + "unemployment_compensation_total": US_NATIONAL_TOTAL_TARGET_LOSS_MULTIPLIER, } US_CRITICAL_TARGET_FIT_REQUIREMENTS = ( { @@ -257,7 +272,7 @@ def _parse_args() -> argparse.Namespace: parser.add_argument("--out", type=Path, required=True) parser.add_argument("--release-id") parser.add_argument("--epochs", type=int, default=512) - parser.add_argument("--learning-rate", type=float, default=0.12) + parser.add_argument("--learning-rate", type=float, default=0.02) parser.add_argument("--max-weight-ratio", type=float, default=5.0) parser.add_argument("--seed", type=int, default=0) parser.add_argument( @@ -1285,19 +1300,30 @@ def _write_npz(path: Path, *, result, registry: TargetRegistry) -> None: ) -def _target_value_loss_weights(registry: TargetRegistry) -> np.ndarray: - values = np.asarray([abs(spec.value) for spec in registry.specs], dtype=np.float64) - weights = np.maximum(values, 1.0) +def _fiscal_target_loss_weights(registry: TargetRegistry) -> np.ndarray: + weights = np.ones(len(registry.specs), dtype=np.float64) + state_multipliers = np.asarray( + [ + US_STATE_TARGET_LOSS_MULTIPLIER if spec.metadata.get("state_fips") else 1.0 + for spec in registry.specs + ], + dtype=np.float64, + ) multipliers = np.asarray( [ - US_FISCAL_TARGET_ROLE_LOSS_MULTIPLIERS.get( - spec.metadata.get("target_role", ""), - 1.0, + ( + 1.0 + if spec.metadata.get("state_fips") + else US_FISCAL_TARGET_ROLE_LOSS_MULTIPLIERS.get( + spec.metadata.get("target_role", ""), + 1.0, + ) ) for spec in registry.specs ], dtype=np.float64, ) + weights *= state_multipliers weights *= multipliers return weights / weights.mean() @@ -1909,7 +1935,8 @@ def main() -> None: max_weight_ratio=args.max_weight_ratio, seed=args.seed, mass="conserve", - target_loss_weights=_target_value_loss_weights(registry), + target_loss_weights=_fiscal_target_loss_weights(registry), + target_loss_cap=US_FISCAL_TARGET_LOSS_CAP, ) _assert_release_gates( result, @@ -1934,6 +1961,7 @@ def main() -> None: "base_dataset_sha256": _sha256(base_h5), "target_compilation": compilation, "target_loss_weighting": US_FISCAL_TARGET_LOSS_WEIGHTING, + "target_loss_cap": US_FISCAL_TARGET_LOSS_CAP, "target_profile_coverage": { "passed": target_profile_gate.passed, "failures": list(target_profile_gate.failures),