From 013f95456d24956204cbd53f38204a8e0ab78d03 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mar=C3=ADa=20Juaristi?= <127882282+juaristi22@users.noreply.github.com> Date: Tue, 13 Oct 2026 15:47:20 +0100 Subject: [PATCH 1/3] Add disclosure control for UK geography outputs Add policyengine.disclosure, a self-contained array-first layer: a floor on distinct survey households, a dominance rule or keyed noise, bootstrap intervals, precision-adaptive rounding, a release check for numbers a reader could recover by subtraction, and a format validator for published files. Add compute_uk_geography_disclosure as the Simulation adapter, keep source_household_id and region on UK household outputs, and deprecate the two unprotected constituency and local-authority functions. Fixes #589 Co-Authored-By: Claude Opus 5.5 --- changelog.d/589.added.md | 1 + changelog.d/589.changed.md | 1 + docs/_quarto.yml | 1 + docs/disclosure.md | 185 +++++ docs/examples.md | 7 + docs/outputs.md | 23 +- docs/regions.md | 5 + examples/uk_geography_disclosure.py | 84 +++ src/policyengine/disclosure/__init__.py | 65 ++ src/policyengine/disclosure/__main__.py | 37 + src/policyengine/disclosure/areas.py | 205 ++++++ src/policyengine/disclosure/estimate.py | 446 ++++++++++++ src/policyengine/disclosure/noise.py | 102 +++ src/policyengine/disclosure/release.py | 668 ++++++++++++++++++ src/policyengine/disclosure/schema.py | 255 +++++++ src/policyengine/disclosure/support.py | 129 ++++ src/policyengine/disclosure/validate.py | 196 +++++ src/policyengine/outputs/__init__.py | 6 + .../outputs/constituency_impact.py | 13 +- .../outputs/local_authority_impact.py | 13 +- .../outputs/uk_geography_disclosure.py | 267 +++++++ .../outputs/uk_geography_impact.py | 12 + .../tax_benefit_models/uk/model.py | 6 + tests/test_constituency_impact.py | 25 + tests/test_disclosure_areas.py | 123 ++++ tests/test_disclosure_estimate.py | 494 +++++++++++++ tests/test_disclosure_noise.py | 108 +++ tests/test_disclosure_release.py | 512 ++++++++++++++ tests/test_disclosure_schema.py | 154 ++++ tests/test_disclosure_support.py | 157 ++++ tests/test_disclosure_validate.py | 204 ++++++ tests/test_extra_variables.py | 3 + tests/test_local_authority_impact.py | 27 + tests/test_uk_geography_disclosure.py | 273 +++++++ 34 files changed, 4804 insertions(+), 3 deletions(-) create mode 100644 changelog.d/589.added.md create mode 100644 changelog.d/589.changed.md create mode 100644 docs/disclosure.md create mode 100644 examples/uk_geography_disclosure.py create mode 100644 src/policyengine/disclosure/__init__.py create mode 100644 src/policyengine/disclosure/__main__.py create mode 100644 src/policyengine/disclosure/areas.py create mode 100644 src/policyengine/disclosure/estimate.py create mode 100644 src/policyengine/disclosure/noise.py create mode 100644 src/policyengine/disclosure/release.py create mode 100644 src/policyengine/disclosure/schema.py create mode 100644 src/policyengine/disclosure/support.py create mode 100644 src/policyengine/disclosure/validate.py create mode 100644 src/policyengine/outputs/uk_geography_disclosure.py create mode 100644 tests/test_disclosure_areas.py create mode 100644 tests/test_disclosure_estimate.py create mode 100644 tests/test_disclosure_noise.py create mode 100644 tests/test_disclosure_release.py create mode 100644 tests/test_disclosure_schema.py create mode 100644 tests/test_disclosure_support.py create mode 100644 tests/test_disclosure_validate.py create mode 100644 tests/test_uk_geography_disclosure.py diff --git a/changelog.d/589.added.md b/changelog.d/589.added.md new file mode 100644 index 00000000..8921726f --- /dev/null +++ b/changelog.d/589.added.md @@ -0,0 +1 @@ +Added `policyengine.disclosure`, a self-contained layer for statistical disclosure control and uncertainty on area-level outputs, and `compute_uk_geography_disclosure` for UK regions, constituencies and local authorities. Each published number rests on at least 10 distinct survey households (clones counted once through `source_household_id`), carries a 90% bootstrap interval and is rounded to the order of its standard error. Rows are `reported`, `approximate` or `withheld`. A release check reports numbers a reader could recover by subtraction across geography levels, a package and its components, and earlier publications. Keyed noise per survey household is available behind `DisclosureSettings(noise=True)`. `validate_public` and `python -m policyengine.disclosure validate` check the format of published files. UK runs now keep `source_household_id` and `region` on the household output when the dataset provides them. diff --git a/changelog.d/589.changed.md b/changelog.d/589.changed.md new file mode 100644 index 00000000..c06e4f1e --- /dev/null +++ b/changelog.d/589.changed.md @@ -0,0 +1 @@ +Deprecated `compute_uk_constituency_impacts`, `compute_uk_local_authority_impacts` and their `ConstituencyImpact` and `LocalAuthorityImpact` classes. Their output is not disclosure-controlled; use `compute_uk_geography_disclosure` instead. They now emit a `DeprecationWarning` and otherwise behave as before. diff --git a/docs/_quarto.yml b/docs/_quarto.yml index 179f8834..09e631bd 100644 --- a/docs/_quarto.yml +++ b/docs/_quarto.yml @@ -39,6 +39,7 @@ website: - impact-analysis.md - outputs.md - regions.md + - disclosure.md - visualisation.md - section: "Methodology" contents: diff --git a/docs/disclosure.md b/docs/disclosure.md new file mode 100644 index 00000000..2b240c51 --- /dev/null +++ b/docs/disclosure.md @@ -0,0 +1,185 @@ +--- +title: "Disclosure control" +--- + +Local results from survey microdata rest on few households, and a licence +binds every published number. `policyengine.disclosure` turns per-household +results into area tables that follow the licence's rules and carry their +uncertainty. Use it for every UK result below the national level that leaves +your machine. + +```python +from policyengine.outputs import compute_uk_geography_disclosure + +impact = compute_uk_geography_disclosure( + baseline_simulation=baseline, + reform_simulation=reform, + year=2026, + levels=("region", "constituency"), + reform="budget", +) +impact.findings # what blocks publication, if anything +impact.public # the rows that may be published, or None while blocked +``` + +The UK datasets derive from surveys under the UK Data Service End User +Licence. The defaults below follow its disclosure rules. You can raise them +for a study with stricter conditions; you cannot lower them. + +## What a published row holds + +Each row is one outcome for one area, reform and year. + +| Column | Content | +|---|---| +| `level`, `area_code`, `reform`, `year`, `outcome` | The key of the row | +| `value`, `lo`, `hi` | The rounded estimate and its 90% interval | +| `direction` | `gain`, `loss`, `too_close_to_call` or `no_change` | +| `tier` | `reported`, `approximate` or `withheld` | +| `reason` | `thin_support` for an approximate row, `disclosure_control` for a withheld row | +| `support_band` | `<10`, `10–29` or `30+` survey households, never an exact count | +| `noise` | Whether the row carries keyed noise | +| `estimator` | `direct`: the estimate comes from the area's own records | + +The outcomes are `average_change`, `relative_change`, `total_change`, +`share_gaining`, `share_losing` and `share_no_change`. A withheld row carries +no value, interval, direction or support band. + +## The rules + +**Count floor.** A number must rest on at least 10 distinct survey households +that contribute to it. For a change, a household contributes when its income +changes by more than `change_tolerance` (£1 a year). Copies and clones of one +survey household count once, through the dataset's `source_household_id`. For +a share, both sides must pass: the households that gain and those that do +not. + +**Dominance.** While noise is off, a number is withheld when one survey +household carries more than half of it. + +**Effective sample.** Weights that differ a lot make an estimate less precise +than its household count suggests. Where each record sits in one area, a low +effective sample makes a row `approximate`: it keeps its value and interval, +but the interval understates the uncertainty. Where a record carries weight in +several areas (an area × household weight matrix), the household count passes +by construction, so the effective sample must also be at least 10. + +**Zero cells.** A cell with no contributing household is published as zero. + +**Linked outcomes.** Average, relative and total change share a numerator, so +they are published or withheld together. Gaining, losing and no change add to +one, so the third share is published only when it is safe to derive. + +## Uncertainty and rounding + +Intervals come from a Poisson bootstrap over survey households with 1,000 +replicates. A ratio resamples its numerator and denominator together. The +weights stay fixed, so the intervals leave out calibration uncertainty and are +too narrow, most of all where support is thin. + +Values are rounded to the order of magnitude of their standard error. The +interval is rounded outward on the same base. + +## The release check + +Published numbers constrain each other. If a region is public and all but one +of its constituencies are public, subtraction gives the last one. The check +looks for such remainders: + +- child areas against their parent area (constituencies in regions, regions in + the nation); +- the components of a package against the package; +- earlier publications, passed as `prior_public`. + +`compute_uk_geography_disclosure` always checks constituencies and local +authorities against their region, and regions against the nation, also when +you publish one level only. A parent row counts as public when it passes the +rules. + +The check reports; it does not withhold cells by itself. A blocking finding +names the group and suggests a cell. Withhold it, then check again. When a +finding suggests no cell, withhold the parent (`finding.group` at +`finding.group_level`) instead: + +```python +for finding in impact.findings: + print(finding.code, finding.group, finding.suggested) + +impact.release.withhold( + level="constituency", + area_code="E14001234", + reform="budget", + year=2026, + outcome="average_change", +) +impact.refresh() # checks again and fills impact.public once nothing blocks +``` + +Two findings can be accepted after review, with `acknowledge=[...]`: +`unchecked_overlap` (levels that overlap without nesting, such as +constituencies and local authorities, are published together) and +`noise_reference_public` (an earlier table shows the same cells under another +noise setting). + +## Noise + +Noise is off by default. With `DisclosureSettings(noise=True)`, each survey +household's change is scaled by a fixed factor between 8% and 12% away from 1, +derived from its id and a secret key. The same household has the same factor +in every table and every rebuild. Noise replaces the dominance rule; the count +floor still applies. + +- **Size.** The noise in a total is about a tenth of that total's sampling + error. +- **Key.** Pass `noise_key=` or set `POLICYENGINE_DISCLOSURE_NOISE_KEY`. There + is no default key. Keep the key out of repositories and CI logs, and do not + change it between releases of one dataset: two noised versions of a number + average towards the truth. +- **Limits.** Denominators carry no noise. The factor cancels in a ratio of + two numbers that the same household dominates, and it cancels in every + subtraction, so noise never protects a remainder. + +## Validating published files + +`validate_public` checks the format of a table, a CSV file or a directory. It +needs no key and no microdata, so it can run on every pull request: + +```bash +python -m policyengine.disclosure validate public/data --reviewed notes.csv +``` + +It rejects private audit columns, unknown tiers, a value in a withheld row, +unrounded values and files it does not recognise. It cannot tell whether a +number rested on enough households; only the release that made the file can. + +## The private audit + +`audit_dir=` writes the exact values, household counts and shares behind every +row. The directory must be outside any git repository. + +## Arrays without a simulation + +The package works on arrays and imports nothing from the rest of the library: + +```python +from policyengine.disclosure import AreaWeights, Release, estimate_change + +seats = estimate_change( + level="constituency", + reform="budget", + year=2026, + areas=AreaWeights.from_labels(seat_codes, household_weights), + source_ids=source_household_ids, + baseline=baseline_net_income, + reformed=reformed_net_income, +) +release = Release() +release.add(seats) +release.check() +public = release.to_public() +``` + +`AreaWeights.from_matrix(matrix, area_codes)` takes an area × household weight +matrix instead. A dataset without `source_household_id` needs `source_ids=` +from the caller; the layer does not guess which records share a survey +household. diff --git a/docs/examples.md b/docs/examples.md index c1417a59..fba2cb0c 100644 --- a/docs/examples.md +++ b/docs/examples.md @@ -16,6 +16,13 @@ Applying parametric reforms, comparing baseline and reform with `ChangeAggregate ```{.python include="../examples/policy_change_uk.py"} ``` +## UK geography results under disclosure control + +Estimating results by area from per-household arrays, checking the release, withholding a cell and validating the public table. The data is synthetic. + +```{.python include="../examples/uk_geography_disclosure.py"} +``` + ## UK income bands Calculating net income and tax by income decile using representative microdata and `Aggregate` with quantile filters. diff --git a/docs/outputs.md b/docs/outputs.md index ff1708e1..33611631 100644 --- a/docs/outputs.md +++ b/docs/outputs.md @@ -260,7 +260,28 @@ for row in impacts.district_results: print(row["district_geoid"], row["avg_change"], row["winner_percentage"]) ``` -### UK constituencies / local authorities +### UK regions, constituencies and local authorities + +`compute_uk_geography_disclosure` returns results by area under disclosure +control: each number passes a floor on the survey households behind it and +carries a 90% interval. Use it for anything you publish. See +[Disclosure control](disclosure.md). + +```python +from policyengine.outputs import compute_uk_geography_disclosure + +impact = compute_uk_geography_disclosure( + baseline_simulation=baseline, + reform_simulation=reform, + year=2026, + levels=("region", "constituency"), +) +impact.public.to_frame() +``` + +The two older functions below are deprecated. Their output is not +disclosure-controlled, so do not publish it when the data is licensed survey +microdata. Constituency and local-authority breakdowns group household output rows by dataset-provided longwise geography columns. Constituencies use diff --git a/docs/regions.md b/docs/regions.md index 3371b76f..a6f3e237 100644 --- a/docs/regions.md +++ b/docs/regions.md @@ -88,6 +88,11 @@ national UK dataset. ## UK parliamentary constituencies +To publish results by constituency, local authority or region, use +`compute_uk_geography_disclosure`; see [Disclosure control](disclosure.md). +The functions in this section and the next are deprecated: their output is +not disclosure-controlled. + Constituency-level impacts group household output rows by the longwise `constituency_code_oa` column carried by the dataset. If the constituency CSV is available locally or from the PolicyEngine UK GCS bucket, PolicyEngine uses it diff --git a/examples/uk_geography_disclosure.py b/examples/uk_geography_disclosure.py new file mode 100644 index 00000000..d84dbc24 --- /dev/null +++ b/examples/uk_geography_disclosure.py @@ -0,0 +1,84 @@ +"""Example: Publish area results under disclosure control. + +This script demonstrates: +1. Building estimates for two geography levels from per-household arrays +2. Checking a release for numbers a reader could recover by subtraction +3. Resolving a finding by withholding the suggested cell +4. Writing the public table and validating its format + +The data is synthetic, so the script needs no dataset. + +Run: python examples/uk_geography_disclosure.py +""" + +import tempfile +from pathlib import Path + +import numpy as np + +from policyengine.disclosure import ( + AreaWeights, + Release, + estimate_change, + validate_public, +) + + +def main() -> Path: + rng = np.random.default_rng(0) + # Two regions of three seats, 400 households a seat. In seat S1 the + # reform reaches only a handful of households. + seats = np.repeat(["S1", "S2", "S3", "S4", "S5", "S6"], 400) + regions = np.repeat(["North", "South"], 1_200) + reach = np.where(seats == "S1", 0.01, np.where(seats == "S2", 0.08, 0.35)) + affected = rng.random(len(seats)) < reach + weights = rng.uniform(150.0, 450.0, len(seats)) + baseline = rng.lognormal(10.4, 0.5, len(seats)) + reformed = baseline + np.where(affected, rng.uniform(100.0, 900.0, len(seats)), 0) + source_ids = np.arange(len(seats)) + + release = Release() + for level, codes in (("region", regions), ("constituency", seats)): + release.add( + estimate_change( + level=level, + reform="example", + year=2026, + areas=AreaWeights.from_labels(codes, weights), + source_ids=source_ids, + baseline=baseline, + reformed=reformed, + ) + ) + release.add_nesting("constituency", "region") + + # S1 is withheld. Its region is public, so subtraction would recover it + # unless a second seat is withheld too. + for finding in release.check(): + print(f"{finding.code}: {finding.message} Suggested: {finding.suggested}") + release.withhold( + level="constituency", + area_code=finding.suggested, + reform="example", + year=2026, + outcome="average_change" if finding.family == "change" else finding.family, + ) + assert release.check() == [] + + public = release.to_public() + frame = public.to_frame() + print( + frame[frame["outcome"] == "average_change"][ + ["level", "area_code", "value", "lo", "hi", "tier", "support_band"] + ].to_string(index=False) + ) + + output = Path(tempfile.mkdtemp(prefix="disclosure-example-")) + public.to_csv(output / "example.csv") + assert validate_public(output) == [] + print(f"Wrote and validated {output / 'example.csv'}") + return output + + +if __name__ == "__main__": + main() diff --git a/src/policyengine/disclosure/__init__.py b/src/policyengine/disclosure/__init__.py new file mode 100644 index 00000000..b4b6885b --- /dev/null +++ b/src/policyengine/disclosure/__init__.py @@ -0,0 +1,65 @@ +"""Disclosure control and uncertainty for area-level outputs. + +Survey microdata under a licence binds every published number. This package +turns per-record results into area tables that carry the licence's rules: a +floor on the number of distinct survey households behind a number, a dominance +rule or keyed noise, bootstrap intervals, precision-adaptive rounding, and a +check on what a reader can recover by subtraction. + +The package is self-contained. It works on arrays and imports nothing from +the rest of the library, so it can be lifted into its own distribution. + +.. code-block:: python + + from policyengine.disclosure import AreaWeights, Release, estimate_change + + seats = estimate_change( + level="constituency", + reform="budget", + year=2026, + areas=AreaWeights.from_labels(seat_codes, household_weights), + source_ids=source_household_ids, + baseline=baseline_net_income, + reformed=reformed_net_income, + ) + release = Release() + release.add(seats) + findings = release.check() + public = release.to_public() # raises while a blocking finding remains +""" + +from .areas import UNASSIGNED as UNASSIGNED +from .areas import AreaWeights as AreaWeights +from .estimate import AreaEstimates as AreaEstimates +from .estimate import estimate_change as estimate_change +from .noise import NOISE_KEY_ENV as NOISE_KEY_ENV +from .noise import DisclosureKeyError as DisclosureKeyError +from .noise import noise_factors as noise_factors +from .noise import resolve_noise_key as resolve_noise_key +from .release import DisclosureError as DisclosureError +from .release import Release as Release +from .schema import DisclosureSettings as DisclosureSettings +from .schema import Finding as Finding +from .schema import PriorPublic as PriorPublic +from .schema import PublicRow as PublicRow +from .schema import PublicTable as PublicTable +from .validate import validate_public as validate_public + +__all__ = [ + "NOISE_KEY_ENV", + "UNASSIGNED", + "AreaEstimates", + "AreaWeights", + "DisclosureError", + "DisclosureKeyError", + "DisclosureSettings", + "Finding", + "PriorPublic", + "PublicRow", + "PublicTable", + "Release", + "estimate_change", + "noise_factors", + "resolve_noise_key", + "validate_public", +] diff --git a/src/policyengine/disclosure/__main__.py b/src/policyengine/disclosure/__main__.py new file mode 100644 index 00000000..ec8d3983 --- /dev/null +++ b/src/policyengine/disclosure/__main__.py @@ -0,0 +1,37 @@ +"""Command line: ``python -m policyengine.disclosure validate ...``.""" + +import argparse +import sys +from collections.abc import Sequence +from typing import Optional + +from .validate import validate_public + + +def main(argv: Optional[Sequence[str]] = None) -> int: + parser = argparse.ArgumentParser(prog="python -m policyengine.disclosure") + commands = parser.add_subparsers(dest="command", required=True) + validate = commands.add_parser( + "validate", help="Check the format of published tables." + ) + validate.add_argument("paths", nargs="+", help="CSV files or directories.") + validate.add_argument( + "--reviewed", + action="append", + default=[], + metavar="NAME", + help="A file in a directory that a person checked. Repeat as needed.", + ) + arguments = parser.parse_args(argv) + + problems: list[str] = [] + for path in arguments.paths: + problems.extend(validate_public(path, reviewed=arguments.reviewed)) + for problem in problems: + print(problem) + print(f"{len(problems)} problem(s) found.") + return 1 if problems else 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/policyengine/disclosure/areas.py b/src/policyengine/disclosure/areas.py new file mode 100644 index 00000000..9ebe1152 --- /dev/null +++ b/src/policyengine/disclosure/areas.py @@ -0,0 +1,205 @@ +"""Area x record weights. + +Two layouts hold the same thing. In the label layout each record sits in one +area with one weight. In the matrix layout a record can carry weight in many +areas. The rules read both through the same operations, so they are written +once. A matrix with at most one non-zero weight per record is stored in the +label layout. +""" + +from collections.abc import Iterable, Sequence +from typing import Optional + +import numpy as np +import pandas as pd + +UNASSIGNED = "" + + +def normalise_area_codes(area_codes: object) -> np.ndarray: + """Return codes as stripped strings, with missing codes as ``""``.""" + + def _one(value: object) -> str: + if value is None: + return UNASSIGNED + if isinstance(value, bytes): + return value.decode("utf-8").strip() + if isinstance(value, float) and np.isnan(value): + return UNASSIGNED + if value is pd.NA: + return UNASSIGNED + return str(value).strip() + + return np.array([_one(value) for value in pd.Series(area_codes)], dtype=object) + + +def _validated_weights(weights: object) -> np.ndarray: + array = np.asarray(weights, dtype=float) + if array.ndim != 1: + raise ValueError("Weights must be one-dimensional.") + if not np.isfinite(array).all(): + raise ValueError("Weights must be finite.") + if (array < 0).any(): + raise ValueError("Weights must be non-negative.") + return array + + +class AreaWeights: + """The weight of every record in every area of one geography level.""" + + def __init__( + self, + *, + area_codes: Sequence[str], + n_records: int, + label_index: Optional[np.ndarray] = None, + weights: Optional[np.ndarray] = None, + matrix: Optional[np.ndarray] = None, + ) -> None: + self.area_codes = tuple(area_codes) + self.n_records = int(n_records) + self._position = {code: i for i, code in enumerate(self.area_codes)} + # The label layout fills the first two; the matrix layout the third. + self._label_index = ( + np.empty(0, dtype=int) if label_index is None else label_index + ) + self._weights = np.empty(0) if weights is None else weights + self._matrix = matrix + if label_index is not None: + # Records sorted by area, so each area is one slice. + assigned = np.flatnonzero(label_index >= 0) + order = assigned[np.argsort(label_index[assigned], kind="stable")] + self._order = order + counts = np.bincount(label_index[assigned], minlength=len(self.area_codes)) + self._starts = np.concatenate(([0], np.cumsum(counts)[:-1])).astype(int) + self._counts = counts + + @classmethod + def from_labels(cls, area_codes: object, weights: object) -> "AreaWeights": + """One area and one weight per record. Blank codes stay unassigned.""" + codes = normalise_area_codes(area_codes) + weight_array = _validated_weights(weights) + if len(codes) != len(weight_array): + raise ValueError("Area codes and weights differ in length.") + return cls._from_codes(codes, weight_array) + + @classmethod + def _from_codes( + cls, + codes: np.ndarray, + weights: np.ndarray, + order: Optional[Sequence[str]] = None, + ) -> "AreaWeights": + weighted = {code for code in codes[weights > 0] if code != UNASSIGNED} + ordered = order if order is not None else sorted(weighted) + area_codes = [code for code in ordered if code in weighted] + position = {code: i for i, code in enumerate(area_codes)} + # A record with no code, or in an area with no weight, is unassigned. + label_index = np.array([position.get(code, -1) for code in codes], dtype=int) + return cls( + area_codes=area_codes, + n_records=len(codes), + label_index=label_index, + weights=weights, + ) + + @classmethod + def from_matrix(cls, matrix: object, area_codes: Sequence[object]) -> "AreaWeights": + """An areas x records weight matrix. The matrix is not copied.""" + array = np.asarray(matrix) + if array.ndim != 2: + raise ValueError("A weight matrix must be two-dimensional.") + codes = [str(code).strip() for code in normalise_area_codes(area_codes)] + if array.shape[0] != len(codes): + raise ValueError("A weight matrix needs one row per area code.") + if len(set(codes)) != len(codes) or UNASSIGNED in codes: + raise ValueError("Matrix area codes must be unique and non-blank.") + if not np.isfinite(array).all(): + raise ValueError("Weights must be finite.") + if (array < 0).any(): + raise ValueError("Weights must be non-negative.") + if array.shape[0] and np.count_nonzero(array, axis=0).max(initial=0) <= 1: + # One non-zero per record: the same data as one label per record. + weights = array.max(axis=0).astype(float) + row = array.argmax(axis=0) + labels = np.array( + [codes[r] if w > 0 else UNASSIGNED for r, w in zip(row, weights)], + dtype=object, + ) + return cls._from_codes(labels, weights, order=codes) + return cls(area_codes=codes, n_records=array.shape[1], matrix=array) + + @property + def n_areas(self) -> int: + return len(self.area_codes) + + @property + def one_area_per_record(self) -> bool: + return self._matrix is None + + @property + def labels(self) -> np.ndarray: + """The area code of each record (``""`` if unassigned).""" + if self._matrix is not None: + raise ValueError("A weight matrix has no single area per record.") + codes = np.array(self.area_codes + (UNASSIGNED,), dtype=object) + return np.asarray(codes[self._label_index]) + + def area_records(self, area: str) -> tuple[np.ndarray, np.ndarray]: + """The records with weight in one area, and their weights.""" + position = self._position[area] + if self._matrix is not None: + row = self._matrix[position] + records = np.flatnonzero(row) + return records, row[records].astype(float) + start = self._starts[position] + records = self._order[start : start + self._counts[position]] + weights = self._weights[records] + keep = weights > 0 + return records[keep], weights[keep] + + def record_weights( + self, areas: Iterable[str], *, include_unassigned: bool = False + ) -> np.ndarray: + """Each record's weight, summed over the given areas.""" + positions = [self._position[area] for area in areas] + if self._matrix is not None: + if not positions: + return np.zeros(self.n_records) + return self._matrix[positions].sum(axis=0, dtype=float) + mask = np.isin(self._label_index, positions) + if include_unassigned: + mask |= self._label_index < 0 + return np.where(mask, self._weights, 0.0) + + def totals( + self, values: object, multipliers: Optional[np.ndarray] = None + ) -> np.ndarray: + """Weighted sums of a per-record value, by area. + + With ``multipliers`` (replicates x records) the result has one column + per replicate. + """ + value_array = np.asarray(values, dtype=float) + if value_array.shape != (self.n_records,): + raise ValueError("Values must have one entry per record.") + if self._matrix is not None: + dtype = self._matrix.dtype + if multipliers is None: + return np.asarray(self._matrix @ value_array.astype(dtype), dtype=float) + # One product per chunk of replicates; no matrix-sized temporary. + scaled = np.ascontiguousarray( + (multipliers * value_array[None, :]).T, dtype=dtype + ) + return np.asarray(self._matrix @ scaled, dtype=float) + weighted = self._weights * value_array + if multipliers is None: + return np.bincount( + self._label_index[self._order], + weights=weighted[self._order], + minlength=self.n_areas, + ) + if self.n_areas == 0: + return np.zeros((0, multipliers.shape[0])) + by_area = multipliers[:, self._order] * weighted[self._order] + return np.add.reduceat(by_area, self._starts, axis=1).T diff --git a/src/policyengine/disclosure/estimate.py b/src/policyengine/disclosure/estimate.py new file mode 100644 index 00000000..1244ff54 --- /dev/null +++ b/src/policyengine/disclosure/estimate.py @@ -0,0 +1,446 @@ +"""Private estimates for one geography level, reform and year. + +Everything here is private: exact values, counts and shares. Only a +:class:`~.release.Release` turns estimates into rows that may be published. +""" + +import math +import warnings +from typing import Optional, Union + +import numpy as np +import pandas as pd + +from .areas import AreaWeights +from .noise import canonical_source_ids, noise_factors, resolve_noise_key +from .schema import ( + APPROXIMATE, + AVERAGE_CHANGE, + CHANGE_OUTCOMES, + FAMILIES, + FAMILY_CHANGE, + FAMILY_GAINING, + FAMILY_LOSING, + GAIN, + LOSS, + NO_CHANGE, + PRIVATE_SHARE_SPLIT, + PRIVATE_UNDEFINED, + REASON_THIN_SUPPORT, + RELATIVE_CHANGE, + REPORTED, + SHARE_GAINING, + SHARE_LOSING, + SHARE_NO_CHANGE, + TOO_CLOSE_TO_CALL, + TOTAL_CHANGE, + WITHHELD, + DisclosureSettings, +) +from .support import ( + Support, + binding_support, + cell_tier, + family_supports, + primary_tier, + support, + support_band, +) + +_REPLICATE_CHUNK = 100 +_TIER_ORDER = {REPORTED: 0, APPROXIMATE: 1, WITHHELD: 2} + + +class AreaEstimates: + """Estimates and support for the areas of one level. Private.""" + + def __init__( + self, + *, + level: str, + reform: str, + year: int, + areas: AreaWeights, + source_index: np.ndarray, + settings: DisclosureSettings, + noise: bool, + contributions: dict[str, np.ndarray], + no_change: np.ndarray, + table: pd.DataFrame, + ) -> None: + self.level = level + self.reform = reform + self.year = year + self.areas = areas + self.source_index = source_index + self.settings = settings + self.noise = noise + self.contributions = contributions + self.no_change = no_change + self.table = table + + def __repr__(self) -> str: + # The table holds exact values and counts, so it stays out of reprs. + return ( + f"AreaEstimates(level={self.level!r}, reform={self.reform!r}, " + f"year={self.year!r}, areas={self.areas.n_areas}, private)" + ) + + def family_tiers(self, family: str) -> dict[str, str]: + """The primary tier of each area's cell in one family.""" + outcome = AVERAGE_CHANGE if family == FAMILY_CHANGE else family + rows = self.table[self.table["outcome"] == outcome] + return dict(zip(rows["area_code"], rows["family_tier"])) + + def family_sources(self, family: str) -> dict[str, int]: + outcome = AVERAGE_CHANGE if family == FAMILY_CHANGE else family + rows = self.table[self.table["outcome"] == outcome] + return dict(zip(rows["area_code"], rows["n_sources"])) + + +def _per_record(values: object, n_records: int, name: str) -> np.ndarray: + array = np.asarray(values, dtype=float) + if array.shape != (n_records,): + raise ValueError(f"{name} must have one entry per record.") + if not np.isfinite(array).all(): + raise ValueError(f"{name} must be finite.") + return array + + +def _rounding_base(standard_error: float, minimum: float) -> float: + """The order of magnitude of the standard error, never below the minimum.""" + if not math.isfinite(standard_error) or standard_error <= 0: + return minimum + return max(10.0 ** math.floor(math.log10(standard_error)), minimum) + + +def _round(value: float, base: float, mode: str = "nearest") -> float: + quotient = value / base + if mode == "down": + steps = math.floor(quotient + 1e-9) + elif mode == "up": + steps = math.ceil(quotient - 1e-9) + else: + steps = round(quotient) + digits = max(0, -int(round(math.log10(base)))) + return round(steps * base, digits) + 0.0 + + +def _replicates( + areas: AreaWeights, + contributions: dict[str, np.ndarray], + baseline: np.ndarray, + source_index: np.ndarray, + n_sources: int, + settings: DisclosureSettings, +) -> dict[str, np.ndarray]: + """Poisson bootstrap over source households: areas x replicates per outcome. + + A replicate reweights every record of a source household by one draw. A + ratio resamples its numerator and denominator with the same draws. + """ + rng = np.random.default_rng(settings.bootstrap_seed) + ones = np.ones(areas.n_records) + out: dict[str, list[np.ndarray]] = { + outcome: [] + for outcome in (*CHANGE_OUTCOMES, SHARE_GAINING, SHARE_LOSING, SHARE_NO_CHANGE) + } + done = 0 + while done < settings.n_replicates: + size = min(_REPLICATE_CHUNK, settings.n_replicates - done) + done += size + # Replicate-major draws, so the result does not depend on the chunk. + draws = rng.poisson(1.0, size=(size, n_sources)).astype(np.uint8) + multipliers = draws[:, source_index] + change = areas.totals(contributions[FAMILY_CHANGE], multipliers) + households = areas.totals(ones, multipliers) + income = areas.totals(baseline, multipliers) + gaining = areas.totals(contributions[FAMILY_GAINING], multipliers) + losing = areas.totals(contributions[FAMILY_LOSING], multipliers) + with np.errstate(divide="ignore", invalid="ignore"): + share_gaining = np.clip(gaining / households, 0.0, 1.0) + share_losing = np.clip(losing / households, 0.0, 1.0 - share_gaining) + out[AVERAGE_CHANGE].append(change / households) + out[RELATIVE_CHANGE].append(np.where(income > 0, change / income, np.nan)) + out[TOTAL_CHANGE].append(change) + out[SHARE_GAINING].append(share_gaining) + out[SHARE_LOSING].append(share_losing) + out[SHARE_NO_CHANGE].append(1.0 - share_gaining - share_losing) + return {outcome: np.concatenate(chunks, axis=1) for outcome, chunks in out.items()} + + +def _direction(value: float, lo: float, hi: float) -> str: + if lo > 0: + return GAIN + if hi < 0: + return LOSS + if value == 0 and lo == 0 and hi == 0: + return NO_CHANGE + return TOO_CLOSE_TO_CALL + + +def estimate_change( + *, + level: str, + reform: str, + year: int, + areas: AreaWeights, + source_ids: object, + baseline: object, + reformed: object, + settings: Optional[DisclosureSettings] = None, + noise_key: Optional[Union[str, bytes]] = None, +) -> AreaEstimates: + """Estimate the change between two runs for every area of one level. + + ``baseline`` and ``reformed`` hold one value per record, for example + household net income. ``source_ids`` names the survey household each + record was copied or cloned from. + """ + settings = settings or DisclosureSettings() + # Resolve the key first, so a missing key fails before any computation. + key = resolve_noise_key(noise_key) if settings.noise else None + n_records = areas.n_records + baseline_array = _per_record(baseline, n_records, "baseline") + reformed_array = _per_record(reformed, n_records, "reformed") + canonical = canonical_source_ids(source_ids) + if len(canonical) != n_records: + raise ValueError("source_ids must have one entry per record.") + source_index, sources = pd.factorize(canonical) + n_sources = len(sources) + + change = reformed_array - baseline_array + change = np.where(np.abs(change) > settings.change_tolerance, change, 0.0) + # Noise goes on the change only, never on a denominator. + factor = ( + noise_factors(canonical, key, distance=settings.noise_distance) + if key is not None + else np.ones(n_records) + ) + contributions = { + FAMILY_CHANGE: factor * change, + FAMILY_GAINING: factor * (change > 0), + FAMILY_LOSING: factor * (change < 0), + } + no_change = (change == 0).astype(float) + + households = areas.totals(np.ones(n_records)) + income = areas.totals(baseline_array) + total = areas.totals(contributions[FAMILY_CHANGE]) + with np.errstate(divide="ignore", invalid="ignore"): + share_gaining = np.clip( + areas.totals(contributions[FAMILY_GAINING]) / households, 0.0, 1.0 + ) + share_losing = np.clip( + areas.totals(contributions[FAMILY_LOSING]) / households, + 0.0, + 1.0 - share_gaining, + ) + exact = { + AVERAGE_CHANGE: total / households, + RELATIVE_CHANGE: np.where(income > 0, total / income, np.nan), + TOTAL_CHANGE: total, + SHARE_GAINING: share_gaining, + SHARE_LOSING: share_losing, + SHARE_NO_CHANGE: 1.0 - share_gaining - share_losing, + } + + replicates = _replicates( + areas, contributions, baseline_array, source_index, n_sources, settings + ) + alpha = (1.0 - settings.interval) / 2.0 + bounds: dict[str, tuple[np.ndarray, np.ndarray, np.ndarray]] = {} + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) + for outcome, draws in replicates.items(): + lower = np.nanpercentile(draws, 100.0 * alpha, axis=1) + upper = np.nanpercentile(draws, 100.0 * (1.0 - alpha), axis=1) + bounds[outcome] = ( + np.fmin(lower, exact[outcome]), + np.fmax(upper, exact[outcome]), + np.nanstd(draws, axis=1, ddof=1), + ) + + one_area = areas.one_area_per_record + rules = {"one_area_per_record": one_area, "noise": settings.noise} + rows: list[dict] = [] + for position, code in enumerate(areas.area_codes): + records, weights = areas.area_records(code) + sources_here = source_index[records] + supports: dict[str, tuple[Support, ...]] = {} + tiers: dict[str, str] = {} + reasons: dict[str, Optional[str]] = {} + for family in FAMILIES: + supports[family] = family_supports( + weights, contributions[family][records], sources_here, family + ) + tiers[family], reasons[family] = cell_tier( + supports[family], settings, **rules + ) + + # Gaining, losing and no change add to one. Two published shares give + # the third, so the third cell must also pass. + unchanged = support(weights, no_change[records], sources_here) + third_tier, _ = primary_tier(unchanged, settings, **rules) + both = tiers[FAMILY_GAINING] != WITHHELD and tiers[FAMILY_LOSING] != WITHHELD + if both and third_tier == WITHHELD: + weaker = min( + (FAMILY_LOSING, FAMILY_GAINING), + key=lambda family: binding_support(supports[family]).n_sources, + ) + tiers[weaker], reasons[weaker] = WITHHELD, PRIVATE_SHARE_SPLIT + both = False + if both: + third = max( + (tiers[FAMILY_GAINING], tiers[FAMILY_LOSING], third_tier), + key=_TIER_ORDER.__getitem__, + ) + third_reason = None if third == REPORTED else REASON_THIN_SUPPORT + else: + third, third_reason = WITHHELD, PRIVATE_SHARE_SPLIT + + rounded: dict[str, tuple[float, float, float, float]] = {} + for outcome in (*CHANGE_OUTCOMES, SHARE_GAINING, SHARE_LOSING): + value = float(exact[outcome][position]) + if not math.isfinite(value): + continue + lo, hi, error = (float(part[position]) for part in bounds[outcome]) + base = _rounding_base(error, settings.min_rounding_base[outcome]) + rounded[outcome] = ( + _round(value, base), + _round(lo, base, "down"), + _round(hi, base, "up"), + base, + ) + for outcome in (SHARE_GAINING, SHARE_LOSING): + value, lo, hi, base = rounded[outcome] + rounded[outcome] = (value, max(lo, 0.0), min(hi, 1.0), base) + # No change is the remainder of the two rounded shares, so the three + # published shares add to one. Its value keeps their finer base: a + # reader can derive it at that precision anyway. Its interval is + # rounded on its own standard error. + base = min(rounded[SHARE_GAINING][3], rounded[SHARE_LOSING][3]) + value = _round(1.0 - rounded[SHARE_GAINING][0] - rounded[SHARE_LOSING][0], base) + value = min(max(value, 0.0), 1.0) + lo, hi, error = (float(part[position]) for part in bounds[SHARE_NO_CHANGE]) + own = max( + _rounding_base(error, settings.min_rounding_base[SHARE_NO_CHANGE]), base + ) + rounded[SHARE_NO_CHANGE] = ( + value, + max(min(_round(lo, own, "down"), value), 0.0), + min(max(_round(hi, own, "up"), value), 1.0), + base, + ) + + for outcome in (*CHANGE_OUTCOMES, *(SHARE_GAINING, SHARE_LOSING)): + family = FAMILY_CHANGE if outcome in CHANGE_OUTCOMES else outcome + rows.append( + _row( + level, + code, + reform, + year, + outcome, + family, + position, + exact, + bounds, + rounded.get(outcome), + tiers[family], + reasons[family], + binding_support(supports[family]), + one_area, + ) + ) + rows.append( + _row( + level, + code, + reform, + year, + SHARE_NO_CHANGE, + SHARE_NO_CHANGE, + position, + exact, + bounds, + rounded[SHARE_NO_CHANGE], + third, + third_reason, + binding_support( + ( + *supports[FAMILY_GAINING], + *supports[FAMILY_LOSING], + unchanged, + ) + ), + one_area, + ) + ) + + return AreaEstimates( + level=level, + reform=reform, + year=int(year), + areas=areas, + source_index=source_index, + settings=settings, + noise=settings.noise, + contributions=contributions, + no_change=no_change, + table=pd.DataFrame(rows), + ) + + +def _row( + level: str, + code: str, + reform: str, + year: int, + outcome: str, + family: str, + position: int, + exact: dict[str, np.ndarray], + bounds: dict[str, tuple[np.ndarray, np.ndarray, np.ndarray]], + rounded: Optional[tuple[float, float, float, float]], + family_tier: str, + family_reason: Optional[str], + cell: Support, + one_area_per_record: bool, +) -> dict: + """One private row: exact and rounded values with the cell's support.""" + tier, reason = family_tier, family_reason + value = lo = hi = base = None + direction = None + if rounded is None: + # No defined value, for example a relative change on zero income. + tier, reason = WITHHELD, PRIVATE_UNDEFINED + else: + value, lo, hi, base = rounded + if outcome in CHANGE_OUTCOMES: + direction = _direction(value, lo, hi) + exact_lo, exact_hi, error = (float(part[position]) for part in bounds[outcome]) + return { + "level": level, + "area_code": code, + "reform": reform, + "year": int(year), + "outcome": outcome, + "family": family, + "exact_value": float(exact[outcome][position]), + "exact_lo": exact_lo, + "exact_hi": exact_hi, + "standard_error": error, + "rounding_base": base, + "value": value, + "lo": lo, + "hi": hi, + "direction": direction, + "tier": tier, + "family_tier": family_tier, + "private_reason": reason, + "n_sources": cell.n_sources, + "n_eff": cell.n_eff, + "top_share": cell.top_share, + "support_band": support_band(cell, one_area_per_record=one_area_per_record), + } diff --git a/src/policyengine/disclosure/noise.py b/src/policyengine/disclosure/noise.py new file mode 100644 index 00000000..81c20a61 --- /dev/null +++ b/src/policyengine/disclosure/noise.py @@ -0,0 +1,102 @@ +"""Keyed noise: one fixed factor for each source household. + +Every source household gets a multiplicative factor a set distance from 1, +derived from its id and a secret key with HMAC-SHA256. The same household +has the same factor in every table and in every rebuild, so repeated +publication cannot average the noise away. +""" + +import hashlib +import hmac +import os +from numbers import Integral, Real +from typing import Optional, Union + +import numpy as np +import pandas as pd + +NOISE_KEY_ENV = "POLICYENGINE_DISCLOSURE_NOISE_KEY" +MIN_KEY_BYTES = 16 +_FACTOR_LABEL = b"policyengine.disclosure.noise-factor.v1\x00" + + +class DisclosureKeyError(ValueError): + """The noise key is missing or unusable. The message never holds the key.""" + + +def resolve_noise_key(key: Optional[Union[str, bytes]] = None) -> bytes: + """Return the noise key from the argument, else from the environment. + + There is no default key: a key in code is a public key. + """ + candidate = key if key is not None else os.environ.get(NOISE_KEY_ENV) + if candidate is None: + raise DisclosureKeyError( + f"Noise is on but no key is set: pass noise_key= or set {NOISE_KEY_ENV}." + ) + raw = candidate.encode("utf-8") if isinstance(candidate, str) else bytes(candidate) + raw = raw.strip() + if not raw: + raise DisclosureKeyError(f"The noise key is empty: check {NOISE_KEY_ENV}.") + if len(raw) < MIN_KEY_BYTES: + raise DisclosureKeyError( + f"The noise key is shorter than {MIN_KEY_BYTES} bytes." + ) + return raw + + +def _canonical_id(value: object) -> str: + if isinstance(value, bytes): + value = value.decode("utf-8") + if isinstance(value, str): + text = value.strip() + if not text: + raise ValueError("A source household id is blank.") + return text + if isinstance(value, bool | np.bool_): + raise ValueError("A source household id is a boolean.") + if isinstance(value, Integral): + return str(int(value)) + if isinstance(value, Real): + number = float(value) + if not np.isfinite(number) or not number.is_integer(): + raise ValueError("A source household id is not a whole number.") + return str(int(number)) + raise ValueError(f"Unsupported source household id type: {type(value).__name__}.") + + +def canonical_source_ids(source_ids: object) -> np.ndarray: + """Return one canonical string per row. + + 123, 123.0, "123" and b"123" name the same household. "0123" does not. + A missing or blank id raises: a row with no source cannot be counted. + """ + series = pd.Series(source_ids) + if series.isna().any(): + raise ValueError("A source household id is missing.") + codes, uniques = pd.factorize(series.astype(object), use_na_sentinel=True) + canonical = np.array([_canonical_id(value) for value in uniques], dtype=object) + return np.asarray(canonical[codes]) + + +def _factor(canonical_id: str, key: bytes, distance: tuple[float, float]) -> float: + digest = hmac.new( + key, _FACTOR_LABEL + canonical_id.encode("utf-8"), hashlib.sha256 + ).digest() + fraction = int.from_bytes(digest[:8], "big") / 2**64 + low, high = distance + magnitude = low + (high - low) * fraction + return 1.0 + magnitude if digest[8] & 1 else 1.0 - magnitude + + +def noise_factors( + source_ids: object, + key: bytes, + *, + distance: tuple[float, float] = (0.08, 0.12), +) -> np.ndarray: + """Return the noise factor of each row's source household.""" + canonical = canonical_source_ids(source_ids) + codes, uniques = pd.factorize(canonical) + factors = np.array([_factor(value, key, distance) for value in uniques]) + return np.asarray(factors[codes]) diff --git a/src/policyengine/disclosure/release.py b/src/policyengine/disclosure/release.py new file mode 100644 index 00000000..a9136b0f --- /dev/null +++ b/src/policyengine/disclosure/release.py @@ -0,0 +1,668 @@ +"""A release: the only object that turns private estimates into public rows. + +A release collects the estimates that will be published together, checks what +a reader could recover by subtraction, and refuses to publish while a blocking +finding remains. The check reports; it does not withhold cells by itself. A +person resolves a finding with :meth:`Release.withhold` and checks again. +""" + +import os +from collections.abc import Iterable, Mapping, Sequence +from itertools import combinations +from pathlib import Path +from typing import Any, Optional, Union + +import numpy as np +import pandas as pd + +from .estimate import AreaEstimates +from .schema import ( + ACKNOWLEDGEABLE_FINDINGS, + ALREADY_PUBLIC, + APPROXIMATE, + CHANGE_OUTCOMES, + FAMILIES, + FAMILY_CHANGE, + FAMILY_GAINING, + FAMILY_LOSING, + NOISE_REFERENCE_PUBLIC, + OUTCOME_FAMILY, + PRIVATE_SECONDARY, + REASON_DISCLOSURE_CONTROL, + REASON_THIN_SUPPORT, + RECOVERABLE_PARENT, + RECOVERABLE_REMAINDER, + SHARE_NO_CHANGE, + UNCHECKED_OVERLAP, + WITHHELD, + DisclosureSettings, + Finding, + PriorPublic, + PublicRow, + PublicTable, +) +from .support import cell_tier, family_supports + +_Key = tuple[str, str, int] +_FAMILY_OUTCOMES = { + FAMILY_CHANGE: CHANGE_OUTCOMES, + FAMILY_GAINING: (FAMILY_GAINING,), + FAMILY_LOSING: (FAMILY_LOSING,), +} + + +class DisclosureError(RuntimeError): + """The release cannot be checked or published as it stands.""" + + +class Release: + """Estimates that are published together, and the checks across them.""" + + def __init__( + self, + settings: Optional[DisclosureSettings] = None, + *, + prior_public: Sequence[PriorPublic] = (), + acknowledge: Iterable[str] = (), + ) -> None: + self.settings = settings or DisclosureSettings() + self._acknowledged = frozenset(acknowledge) + unknown = sorted(self._acknowledged - set(ACKNOWLEDGEABLE_FINDINGS)) + if unknown: + raise ValueError(f"These findings cannot be acknowledged: {unknown}") + self._prior = tuple(prior_public) + self._estimates: dict[_Key, AreaEstimates] = {} + self._published: dict[_Key, bool] = {} + self._nestings: list[tuple[str, str, Optional[dict[str, str]]]] = [] + self._packages: list[tuple[str, tuple[str, ...]]] = [] + self._manual: set[tuple[str, str, int, str, str]] = set() + self._findings: Optional[list[Finding]] = None + + def add(self, estimates: AreaEstimates, *, publish: bool = True) -> None: + """Add the estimates of one level, reform and year. + + With ``publish=False`` the estimates are a reference only: the release + does not emit their rows, and treats the rows that pass the rules as + public elsewhere. + """ + if estimates.settings != self.settings: + raise DisclosureError( + "The estimates were made under other settings than the release's." + ) + key = (estimates.level, estimates.reform, estimates.year) + if key in self._estimates: + raise DisclosureError(f"The release already holds estimates for {key}.") + self._estimates[key] = estimates + self._published[key] = publish + self._findings = None + + def add_nesting( + self, + child_level: str, + parent_level: str, + parent_of: Optional[Mapping[str, str]] = None, + ) -> None: + """Declare that the areas of one level tile the areas of another. + + ``parent_of`` (child area code to parent area code) is needed only when + the child level is a weight matrix and the parent has several areas. + """ + mapping = dict(parent_of) if parent_of is not None else None + self._nestings.append((child_level, parent_level, mapping)) + self._findings = None + + def add_package(self, package: str, components: Sequence[str]) -> None: + """Declare that one reform is the sum of other reforms in the release.""" + self._packages.append((package, tuple(components))) + self._findings = None + + def withhold( + self, *, level: str, area_code: str, reform: str, year: int, outcome: str + ) -> None: + """Withhold a cell by hand, for example to resolve a finding. + + Outcomes that share a numerator are withheld together. + """ + family = OUTCOME_FAMILY.get(outcome) + if family is None: + raise ValueError( + f"Withhold {FAMILY_GAINING!r} or {FAMILY_LOSING!r} instead of " + f"{outcome!r}." + ) + key = (level, reform, int(year)) + estimates = self._estimates.get(key) + if estimates is None or area_code not in estimates.areas.area_codes: + raise DisclosureError(f"The release has no cell for {key}, {area_code!r}.") + self._manual.add((*key, family, area_code)) + self._findings = None + + @property + def findings(self) -> Optional[list[Finding]]: + """The findings of the last check, or None if the release changed since.""" + return self._findings + + def check(self) -> list[Finding]: + """Report what a reader could recover from the rows as they stand.""" + tiers = { + (key, family): self._family_tiers(key, family) + for key in self._estimates + for family in FAMILIES + } + known = self._known_cells() + findings: list[Finding] = [] + self._check_nesting(tiers, known, findings) + self._check_packages(tiers, known, findings) + self._check_overlap(findings) + self._check_prior(tiers, findings) + self._findings = findings + return findings + + def to_public(self) -> PublicTable: + """The rows that may be published. Raises while anything blocks.""" + if self._findings is None: + raise DisclosureError("Run check() before to_public().") + blocking = sorted( + {finding.code for finding in self._findings if finding.blocking} + ) + if blocking: + raise DisclosureError( + f"Blocking findings remain: {blocking}. Resolve them, then run " + "check() again." + ) + rows: list[PublicRow] = [] + for key in sorted(self._estimates): + if not self._published[key]: + continue + estimates = self._estimates[key] + for record, tier, _ in self._effective_rows(key): + common = { + "level": record["level"], + "area_code": record["area_code"], + "reform": record["reform"], + "year": int(record["year"]), + "outcome": record["outcome"], + "noise": estimates.noise, + } + if tier == WITHHELD: + rows.append( + PublicRow( + **common, tier=WITHHELD, reason=REASON_DISCLOSURE_CONTROL + ) + ) + continue + direction = record["direction"] + rows.append( + PublicRow( + **common, + value=float(record["value"]), + lo=float(record["lo"]), + hi=float(record["hi"]), + direction=None if pd.isna(direction) else str(direction), + tier=tier, + reason=REASON_THIN_SUPPORT if tier == APPROXIMATE else None, + support_band=record["support_band"], + ) + ) + return PublicTable(rows=rows) + + def write_audit(self, private_dir: Union[str, Path]) -> Path: + """Write exact values, counts and shares to a private directory. + + The directory must be outside any git repository, so the audit cannot + be committed or published with a bundle. + """ + directory = Path(private_dir).expanduser().resolve() + for folder in (directory, *directory.parents): + if (folder / ".git").exists(): + raise DisclosureError( + "The audit path is inside a git repository. Choose a private " + "path outside any repository." + ) + directory.mkdir(parents=True, exist_ok=True) + frames = [] + for key in sorted(self._estimates): + frame = self._estimates[key].table.drop(columns=["family_tier"]).copy() + effective = list(self._effective_rows(key)) + frame["tier"] = [tier for _, tier, _ in effective] + frame["private_reason"] = [reason for _, _, reason in effective] + frame["published"] = self._published[key] + frames.append(frame) + audit = pd.concat(frames, ignore_index=True) if frames else pd.DataFrame() + path = directory / "disclosure_audit.csv" + descriptor = os.open(path, os.O_WRONLY | os.O_CREAT | os.O_TRUNC, 0o600) + with os.fdopen(descriptor, "w", newline="") as handle: + audit.to_csv(handle, index=False) + os.chmod(path, 0o600) + return path + + # -- tiers --------------------------------------------------------------- + + def _family_tiers(self, key: _Key, family: str) -> dict[str, str]: + tiers = self._estimates[key].family_tiers(family) + for level, reform, year, manual_family, area in self._manual: + if (level, reform, year) == key and manual_family == family: + tiers[area] = WITHHELD + return tiers + + def _effective_rows(self, key: _Key) -> Iterable[tuple[dict, str, Optional[str]]]: + """Each private row with its tier and reason after manual withholding.""" + manual = { + (family, area) + for level, reform, year, family, area in self._manual + if (level, reform, year) == key + } + for record in self._estimates[key].table.to_dict("records"): + outcome, area = record["outcome"], record["area_code"] + families = ( + (FAMILY_GAINING, FAMILY_LOSING) + if outcome == SHARE_NO_CHANGE + else (OUTCOME_FAMILY[outcome],) + ) + reason = record["private_reason"] + reason = None if pd.isna(reason) else reason + if any((family, area) in manual for family in families): + yield record, WITHHELD, PRIVATE_SECONDARY + else: + yield record, record["tier"], reason + + def _known_cells(self) -> set[tuple[str, str, str, int, str]]: + """Family cells an earlier publication already made public.""" + known = set() + for prior in self._prior: + for level, area, reform, year, outcome in prior.cells: + family = OUTCOME_FAMILY.get(outcome) + if family is not None: + known.add((level, area, reform, int(year), family)) + return known + + # -- relations ----------------------------------------------------------- + + def _passes( + self, estimates: AreaEstimates, weights: np.ndarray, family: str + ) -> tuple[bool, bool]: + """Whether a remainder may be recovered, and whether it is empty.""" + records = np.flatnonzero(weights) + supports = family_supports( + weights[records], + estimates.contributions[family][records], + estimates.source_index[records], + family, + ) + tier, _ = cell_tier( + supports, + self.settings, + one_area_per_record=estimates.areas.one_area_per_record, + noise=estimates.noise, + ) + empty = supports[0].n_sources == 0 + return tier != WITHHELD, empty + + @staticmethod + def _require_same_records(first: AreaEstimates, second: AreaEstimates) -> None: + same = first.areas.n_records == second.areas.n_records and np.array_equal( + first.source_index, second.source_index + ) + if not same: + raise DisclosureError( + "A relation needs both sets of estimates on the same records " + "and source households." + ) + + def _groups( + self, + child: AreaEstimates, + parent: AreaEstimates, + parent_of: Optional[dict[str, str]], + ) -> dict[str, tuple[list[str], Optional[np.ndarray]]]: + """For each parent area: its child areas and its records with no child.""" + parents = parent.areas.area_codes + if not child.areas.one_area_per_record: + if len(parents) == 1: + return {parents[0]: (list(child.areas.area_codes), None)} + if parent_of is None: + raise DisclosureError( + f"Nesting a weight matrix ({child.level}) in {parent.level} " + "needs parent_of." + ) + missing = sorted(set(child.areas.area_codes) - set(parent_of)) + if missing: + raise DisclosureError( + f"parent_of lacks {len(missing)} {child.level} areas." + ) + return { + code: ( + [c for c in child.areas.area_codes if parent_of[c] == code], + None, + ) + for code in parents + } + try: + parent_labels = parent.areas.labels + except ValueError as error: + raise DisclosureError( + f"The parent level {parent.level} must have one area per record." + ) from error + child_labels = child.areas.labels + pairs = pd.DataFrame({"parent": parent_labels, "child": child_labels}) + pairs = pairs[pairs["child"] != ""].drop_duplicates() + if pairs["child"].duplicated().any(): + raise DisclosureError( + f"A {child.level} area lies in more than one {parent.level} area." + ) + if (pairs["parent"] == "").any(): + raise DisclosureError( + f"A {child.level} area has records outside every {parent.level} area." + ) + children = pairs.groupby("parent")["child"].apply(list).to_dict() + loose_weights = child.areas.record_weights([], include_unassigned=True) + return { + code: ( + sorted(children.get(code, [])), + np.where(parent_labels == code, loose_weights, 0.0), + ) + for code in parents + } + + def _check_nesting( + self, + tiers: dict[tuple[_Key, str], dict[str, str]], + known: set[tuple[str, str, str, int, str]], + findings: list[Finding], + ) -> None: + for child_level, parent_level, parent_of in self._nestings: + for key in sorted(self._estimates): + level, reform, year = key + if level != child_level or not self._published[key]: + continue + child = self._estimates[key] + parent_key = (parent_level, reform, year) + parent = self._estimates.get(parent_key) + if parent is None: + raise DisclosureError( + f"Nesting {child_level} in {parent_level} needs " + f"{parent_level} estimates for reform {reform!r}, {year}." + ) + self._require_same_records(child, parent) + groups = self._groups(child, parent, parent_of) + for family in FAMILIES: + child_tiers = tiers[(key, family)] + sources = child.family_sources(family) + for code, (children, loose) in groups.items(): + hidden = [ + c + for c in children + if child_tiers[c] == WITHHELD + and (level, c, reform, year, family) not in known + ] + weights = child.areas.record_weights(hidden) + if loose is not None: + weights = weights + loose + passes, empty = self._passes(child, weights, family) + location: dict[str, Any] = { + "level": child_level, + "reform": reform, + "year": year, + "family": family, + "group": code, + "group_level": parent_level, + } + parent_withheld = ( + tiers[(parent_key, family)].get(code) == WITHHELD + and (parent_level, code, reform, year, family) not in known + ) + if parent_withheld: + # There is no public parent to subtract from. The + # risk is the reverse: published children that add + # up to a parent this release withholds. + if empty and self._published[parent_key]: + findings.append( + Finding( + code=RECOVERABLE_PARENT, + blocking=True, + relation="nesting", + message=( + f"{parent_level} {code!r} is withheld " + f"but all its {child_level} cells are " + "published, so addition recovers it." + ), + **location, + ) + ) + continue + if not passes: + # Withholding a cell with no contributor adds + # nothing to the remainder. + shown = [ + c + for c in children + if child_tiers[c] != WITHHELD and sources[c] > 0 + ] + findings.append( + Finding( + code=RECOVERABLE_REMAINDER, + blocking=True, + relation="nesting", + message=( + f"The withheld {child_level} cells in " + f"{parent_level} {code!r} fail the rules " + "together, so subtraction recovers them. " + "Withhold another cell of the group, or " + "the parent." + ), + members=sorted(hidden), + suggested=min( + shown, + key=lambda c: (sources[c], c), + default=None, + ), + **location, + ) + ) + + def _check_packages( + self, + tiers: dict[tuple[_Key, str], dict[str, str]], + known: set[tuple[str, str, str, int, str]], + findings: list[Finding], + ) -> None: + family = FAMILY_CHANGE + for package, components in self._packages: + for key in sorted(self._estimates): + level, reform, year = key + if reform != package or not self._published[key]: + continue + whole = self._estimates[key] + parts = [] + for component in components: + part_key = (level, component, year) + part = self._estimates.get(part_key) + if part is None: + raise DisclosureError( + f"Package {package!r} needs estimates for component " + f"{component!r} at {level}, {year}." + ) + self._require_same_records(whole, part) + if part.areas.area_codes != whole.areas.area_codes: + raise DisclosureError( + "A package and its components need the same areas." + ) + parts.append((part_key, part)) + for area in whole.areas.area_codes: + hidden = [ + (part_key, part) + for part_key, part in parts + if self._published[part_key] + and tiers[(part_key, family)][area] == WITHHELD + and (level, area, part.reform, year, family) not in known + ] + whole_hidden = ( + tiers[(key, family)][area] == WITHHELD + and (level, area, package, year, family) not in known + ) + location: dict[str, Any] = { + "level": level, + "reform": package, + "year": year, + "family": family, + "group": area, + "group_level": level, + } + if whole_hidden and not hidden: + findings.append( + Finding( + code=RECOVERABLE_PARENT, + blocking=True, + relation="package", + message=( + f"Package {package!r} is withheld in {area!r} " + "but all its components are published, so " + "addition recovers it." + ), + **location, + ) + ) + elif hidden and not whole_hidden: + records, weights = whole.areas.area_records(area) + combined = sum( + part.contributions[family][records] for _, part in hidden + ) + supports = family_supports( + weights, combined, whole.source_index[records], family + ) + tier, _ = cell_tier( + supports, + self.settings, + one_area_per_record=whole.areas.one_area_per_record, + noise=whole.noise, + ) + if tier != WITHHELD or supports[0].n_sources == 0: + continue + shown = [ + part + for part_key, part in parts + if tiers[(part_key, family)][area] != WITHHELD + and part.family_sources(family)[area] > 0 + ] + findings.append( + Finding( + code=RECOVERABLE_REMAINDER, + blocking=True, + relation="package", + message=( + f"The withheld components of {package!r} in " + f"{area!r} fail the rules together, so " + "subtraction recovers them. Withhold another " + "component, or the package." + ), + members=sorted(part.reform for _, part in hidden), + suggested=min( + (part.reform for part in shown), + key=lambda name: ( + self._estimates[ + (level, name, year) + ].family_sources(family)[area], + name, + ), + default=None, + ), + **location, + ) + ) + + def _check_overlap(self, findings: list[Finding]) -> None: + """Levels published together need a declared relation between them.""" + related: dict[str, set[str]] = {} + for child_level, parent_level, _ in self._nestings: + related.setdefault(child_level, set()).add(parent_level) + + def nests(child: str, parent: str) -> bool: + seen, queue = set(), [child] + while queue: + current = queue.pop() + for above in related.get(current, ()): + if above == parent: + return True + if above not in seen: + seen.add(above) + queue.append(above) + return False + + by_run: dict[tuple[str, int], set[str]] = {} + for (level, reform, year), publish in self._published.items(): + if publish: + by_run.setdefault((reform, year), set()).add(level) + pairs = { + pair + for levels in by_run.values() + for pair in combinations(sorted(levels), 2) + if not nests(*pair) and not nests(*reversed(pair)) + } + for first, second in sorted(pairs): + findings.append( + Finding( + code=UNCHECKED_OVERLAP, + blocking=UNCHECKED_OVERLAP not in self._acknowledged, + relation="release", + message=( + f"{first} and {second} are published together with no " + "declared nesting. Differences between overlapping areas " + "are not checked." + ), + level=f"{first}|{second}", + ) + ) + + def _check_prior( + self, + tiers: dict[tuple[_Key, str], dict[str, str]], + findings: list[Finding], + ) -> None: + for prior in self._prior: + cells: dict[tuple[_Key, str], set[str]] = {} + for level, area, reform, year, outcome in prior.cells: + family = OUTCOME_FAMILY.get(outcome) + key = (level, reform, int(year)) + if family is not None and self._published.get(key): + cells.setdefault((key, family), set()).add(area) + for (key, family), areas in sorted(cells.items()): + level, reform, year = key + family_tiers = tiers[(key, family)] + present = sorted(areas & set(family_tiers)) + withheld = [a for a in present if family_tiers[a] == WITHHELD] + shown = [a for a in present if family_tiers[a] != WITHHELD] + location: dict[str, Any] = { + "level": level, + "reform": reform, + "year": year, + "family": family, + "group": prior.name, + } + if withheld: + findings.append( + Finding( + code=ALREADY_PUBLIC, + blocking=False, + relation="prior_public", + message=( + f"{prior.name!r} already published cells that " + "this release withholds." + ), + members=withheld, + **location, + ) + ) + if shown and prior.noise != self._estimates[key].noise: + findings.append( + Finding( + code=NOISE_REFERENCE_PUBLIC, + blocking=NOISE_REFERENCE_PUBLIC not in self._acknowledged, + relation="prior_public", + message=( + f"{prior.name!r} published the same cells with a " + "different noise setting. A reader can compare " + "the two and learn the noise factors." + ), + members=shown, + **location, + ) + ) diff --git a/src/policyengine/disclosure/schema.py b/src/policyengine/disclosure/schema.py new file mode 100644 index 00000000..565558b7 --- /dev/null +++ b/src/policyengine/disclosure/schema.py @@ -0,0 +1,255 @@ +"""Settings, public models and constants of the disclosure layer.""" + +import math +from collections.abc import Sequence +from pathlib import Path +from typing import Optional, Union + +import pandas as pd +from pydantic import BaseModel, ConfigDict, Field, field_validator + +REPORTED = "reported" +APPROXIMATE = "approximate" +WITHHELD = "withheld" +TIERS = (REPORTED, APPROXIMATE, WITHHELD) + +# The middle band uses an en dash, as the published tables do. +BAND_UNDER_10 = "<10" +BAND_10_TO_29 = "10–29" +BAND_30_PLUS = "30+" +SUPPORT_BANDS = (BAND_UNDER_10, BAND_10_TO_29, BAND_30_PLUS) + +# Reasons a published table may show. A withheld row always shows the same +# reason: a finer one would tell a reader which rows are secondary. +REASON_THIN_SUPPORT = "thin_support" +REASON_DISCLOSURE_CONTROL = "disclosure_control" +PUBLIC_REASONS = (REASON_THIN_SUPPORT, REASON_DISCLOSURE_CONTROL) + +# Reasons kept in the private audit only. +PRIVATE_BELOW_FLOOR = "below_floor" +PRIVATE_LOW_N_EFF = "low_effective_sample" +PRIVATE_DOMINANT = "dominant_household" +PRIVATE_SHARE_SPLIT = "share_split" +PRIVATE_SECONDARY = "secondary" +PRIVATE_UNDEFINED = "undefined_denominator" + +GAIN = "gain" +LOSS = "loss" +TOO_CLOSE_TO_CALL = "too_close_to_call" +NO_CHANGE = "no_change" +DIRECTIONS = (GAIN, LOSS, TOO_CLOSE_TO_CALL, NO_CHANGE) + +AVERAGE_CHANGE = "average_change" +RELATIVE_CHANGE = "relative_change" +TOTAL_CHANGE = "total_change" +SHARE_GAINING = "share_gaining" +SHARE_LOSING = "share_losing" +SHARE_NO_CHANGE = "share_no_change" +CHANGE_OUTCOMES = (AVERAGE_CHANGE, RELATIVE_CHANGE, TOTAL_CHANGE) +SHARE_OUTCOMES = (SHARE_GAINING, SHARE_LOSING, SHARE_NO_CHANGE) +OUTCOMES = CHANGE_OUTCOMES + SHARE_OUTCOMES + +# Outcomes that share a numerator share a tier. A reader can turn one into the +# other with public denominators, so they are published or withheld together. +FAMILY_CHANGE = "change" +FAMILY_GAINING = SHARE_GAINING +FAMILY_LOSING = SHARE_LOSING +FAMILIES = (FAMILY_CHANGE, FAMILY_GAINING, FAMILY_LOSING) +OUTCOME_FAMILY = { + AVERAGE_CHANGE: FAMILY_CHANGE, + RELATIVE_CHANGE: FAMILY_CHANGE, + TOTAL_CHANGE: FAMILY_CHANGE, + SHARE_GAINING: FAMILY_GAINING, + SHARE_LOSING: FAMILY_LOSING, +} + +ESTIMATOR_DIRECT = "direct" + +PUBLIC_COLUMNS = ( + "level", + "area_code", + "reform", + "year", + "outcome", + "value", + "lo", + "hi", + "direction", + "tier", + "reason", + "support_band", + "noise", + "estimator", +) +# Label columns a publisher may add beside the public columns. +OPTIONAL_LABEL_COLUMNS = ("area_name",) +# Columns of the private audit. None may appear in a published file. +AUDIT_COLUMNS = ( + "family", + "exact_value", + "exact_lo", + "exact_hi", + "standard_error", + "rounding_base", + "n_sources", + "n_eff", + "top_share", + "private_reason", +) + +DEFAULT_MIN_ROUNDING_BASE = { + AVERAGE_CHANGE: 1.0, + RELATIVE_CHANGE: 0.0001, + TOTAL_CHANGE: 1000.0, + SHARE_GAINING: 0.001, + SHARE_LOSING: 0.001, + SHARE_NO_CHANGE: 0.001, +} + +RECOVERABLE_REMAINDER = "recoverable_remainder" +RECOVERABLE_PARENT = "recoverable_parent" +UNCHECKED_OVERLAP = "unchecked_overlap" +NOISE_REFERENCE_PUBLIC = "noise_reference_public" +ALREADY_PUBLIC = "already_public" +# Findings a caller may accept after review. The others always block. +ACKNOWLEDGEABLE_FINDINGS = (UNCHECKED_OVERLAP, NOISE_REFERENCE_PUBLIC) + + +def is_power_of_ten(value: float) -> bool: + if not math.isfinite(value) or value <= 0: + return False + exponent = math.log10(value) + return abs(exponent - round(exponent)) < 1e-9 + + +class DisclosureSettings(BaseModel): + """Publication rules. The floors can be raised, never lowered.""" + + model_config = ConfigDict(frozen=True, extra="forbid") + + min_sources: int = Field(default=10, ge=10) + max_top_share: float = Field(default=0.5, gt=0.0, le=0.5) + min_n_eff: float = Field(default=10.0, ge=10.0) + approximate_below_n_eff: float = Field(default=30.0, ge=0.0) + # A household contributes, gains or loses when its change is larger than + # this. Smaller changes count as no change and are left out of every sum. + change_tolerance: float = Field(default=1.0, gt=0.0) + noise: bool = False + # Distance of a noise factor from 1: uniform on this range, either side. + noise_distance: tuple[float, float] = (0.08, 0.12) + n_replicates: int = Field(default=1000, ge=200) + interval: float = Field(default=0.90, gt=0.5, lt=1.0) + bootstrap_seed: int = 0 + min_rounding_base: dict[str, float] = Field( + default_factory=lambda: dict(DEFAULT_MIN_ROUNDING_BASE) + ) + + @field_validator("noise_distance") + @classmethod + def _noise_distance_keeps_clear_of_one( + cls, value: tuple[float, float] + ) -> tuple[float, float]: + low, high = value + if not 0.0 < low < high < 1.0: + raise ValueError("noise_distance must satisfy 0 < low < high < 1") + return value + + @field_validator("min_rounding_base") + @classmethod + def _rounding_bases_are_complete_powers_of_ten( + cls, value: dict[str, float] + ) -> dict[str, float]: + unknown = sorted(set(value) - set(OUTCOMES)) + if unknown: + raise ValueError(f"unknown outcomes in min_rounding_base: {unknown}") + merged = {**DEFAULT_MIN_ROUNDING_BASE, **value} + for outcome, base in merged.items(): + if not is_power_of_ten(base): + raise ValueError( + f"min_rounding_base[{outcome!r}] must be a power of ten" + ) + return merged + + +class PublicRow(BaseModel): + """One published row: an outcome for an area, a reform and a year.""" + + model_config = ConfigDict(frozen=True, extra="forbid") + + level: str + area_code: str + reform: str + year: int + outcome: str + value: Optional[float] = None + lo: Optional[float] = None + hi: Optional[float] = None + direction: Optional[str] = None + tier: str + reason: Optional[str] = None + support_band: Optional[str] = None + noise: bool = False + estimator: str = ESTIMATOR_DIRECT + + +class PublicTable(BaseModel): + """The rows a release may publish.""" + + model_config = ConfigDict(frozen=True, extra="forbid") + + rows: list[PublicRow] + + def to_frame(self) -> pd.DataFrame: + records = [row.model_dump() for row in self.rows] + return pd.DataFrame(records, columns=list(PUBLIC_COLUMNS)) + + def to_csv(self, path: Union[str, Path]) -> Path: + path = Path(path) + self.to_frame().to_csv(path, index=False) + return path + + +class Finding(BaseModel): + """A result of the release check. It carries labels, never values.""" + + model_config = ConfigDict(frozen=True, extra="forbid") + + code: str + blocking: bool + relation: str + message: str + level: Optional[str] = None + reform: Optional[str] = None + year: Optional[int] = None + family: Optional[str] = None + group: Optional[str] = None + group_level: Optional[str] = None + members: list[str] = Field(default_factory=list) + suggested: Optional[str] = None + + +class PriorPublic(BaseModel): + """Cells an earlier publication already made public.""" + + model_config = ConfigDict(frozen=True, extra="forbid") + + name: str + noise: bool = False + cells: list[tuple[str, str, str, int, str]] + + @classmethod + def from_frame( + cls, name: str, frame: pd.DataFrame, *, noise: bool = False + ) -> "PriorPublic": + """Read cells from columns level, area_code, reform, year, outcome.""" + columns: Sequence[str] = ("level", "area_code", "reform", "year", "outcome") + missing = [column for column in columns if column not in frame.columns] + if missing: + raise ValueError(f"prior publication {name!r} lacks columns {missing}") + cells = [ + (str(level), str(area), str(reform), int(year), str(outcome)) + for level, area, reform, year, outcome in frame[list(columns)].itertuples( + index=False, name=None + ) + ] + return cls(name=name, noise=noise, cells=cells) diff --git a/src/policyengine/disclosure/support.py b/src/policyengine/disclosure/support.py new file mode 100644 index 00000000..7d6b50c5 --- /dev/null +++ b/src/policyengine/disclosure/support.py @@ -0,0 +1,129 @@ +"""Support behind a number: how many source households it rests on. + +One function serves a cell, the two sides of a share, and a remainder that a +reader could recover by subtraction. Clones and copies of one survey household +count once: weights and contributions are summed by source first. +""" + +from typing import NamedTuple, Optional + +import numpy as np + +from .schema import ( + APPROXIMATE, + BAND_10_TO_29, + BAND_30_PLUS, + BAND_UNDER_10, + FAMILY_CHANGE, + PRIVATE_BELOW_FLOOR, + PRIVATE_DOMINANT, + PRIVATE_LOW_N_EFF, + REASON_THIN_SUPPORT, + REPORTED, + WITHHELD, + DisclosureSettings, +) + + +class Support(NamedTuple): + n_sources: int + n_eff: float + top_share: float + + +NO_SUPPORT = Support(0, 0.0, 0.0) + + +def support( + weights: np.ndarray, contributions: np.ndarray, source_index: np.ndarray +) -> Support: + """Support of one cell, from aligned per-record arrays. + + A record contributes when it has weight and a non-zero contribution. + """ + contributing = (weights > 0) & (contributions != 0) + if not contributing.any(): + return NO_SUPPORT + _, by_source = np.unique(source_index[contributing], return_inverse=True) + source_weight = np.bincount(by_source, weights=weights[contributing]) + magnitude = np.bincount( + by_source, + weights=np.abs(weights[contributing] * contributions[contributing]), + ) + n_eff = float(source_weight.sum() ** 2 / (source_weight**2).sum()) + top_share = float(magnitude.max() / magnitude.sum()) + return Support(len(source_weight), n_eff, top_share) + + +def family_supports( + weights: np.ndarray, + contributions: np.ndarray, + source_index: np.ndarray, + family: str, +) -> tuple[Support, ...]: + """The supports a family's cell must pass: one for a change, two for a share.""" + cell = support(weights, contributions, source_index) + if family == FAMILY_CHANGE: + return (cell,) + complement = (contributions == 0).astype(float) + return (cell, support(weights, complement, source_index)) + + +def primary_tier( + cell: Support, + settings: DisclosureSettings, + *, + one_area_per_record: bool, + noise: bool, +) -> tuple[str, Optional[str]]: + """Tier and private reason of one side. The count rule comes first.""" + if cell.n_sources == 0: + return REPORTED, None + if cell.n_sources < settings.min_sources: + return WITHHELD, PRIVATE_BELOW_FLOOR + if not one_area_per_record and cell.n_eff < settings.min_n_eff: + return WITHHELD, PRIVATE_LOW_N_EFF + if not noise and cell.top_share > settings.max_top_share: + return WITHHELD, PRIVATE_DOMINANT + if cell.n_eff < settings.approximate_below_n_eff: + return APPROXIMATE, REASON_THIN_SUPPORT + return REPORTED, None + + +def binding_support(supports: tuple[Support, ...]) -> Support: + """The weakest non-empty side, which sets the band and the audit columns.""" + present = [cell for cell in supports if cell.n_sources > 0] + if not present: + return NO_SUPPORT + return min(present, key=lambda cell: (cell.n_sources, cell.n_eff)) + + +def cell_tier( + supports: tuple[Support, ...], + settings: DisclosureSettings, + *, + one_area_per_record: bool, + noise: bool, +) -> tuple[str, Optional[str]]: + """Tier of a cell from all its sides: the worst side decides.""" + tier: str = REPORTED + reason: Optional[str] = None + for cell in supports: + side_tier, side_reason = primary_tier( + cell, settings, one_area_per_record=one_area_per_record, noise=noise + ) + if side_tier == WITHHELD: + return side_tier, side_reason + if side_tier == APPROXIMATE: + tier, reason = side_tier, side_reason + return tier, reason + + +def support_band(cell: Support, *, one_area_per_record: bool) -> str: + """The published band. It never shows an exact count.""" + basis = cell.n_sources if one_area_per_record else cell.n_eff + if basis < 10: + return BAND_UNDER_10 + if basis < 30: + return BAND_10_TO_29 + return BAND_30_PLUS diff --git a/src/policyengine/disclosure/validate.py b/src/policyengine/disclosure/validate.py new file mode 100644 index 00000000..e6294930 --- /dev/null +++ b/src/policyengine/disclosure/validate.py @@ -0,0 +1,196 @@ +"""Format checks on a published table. + +``validate_public`` needs no key and no microdata, so it can run on every +pull request. It checks necessary conditions only: the shape of the file, not +the support behind it. It cannot tell whether a number rested on enough +households. Only the release that made the file can. +""" + +import math +from collections.abc import Iterable, Mapping +from pathlib import Path +from typing import Optional, Union + +import pandas as pd + +from .schema import ( + APPROXIMATE, + AUDIT_COLUMNS, + BAND_UNDER_10, + CHANGE_OUTCOMES, + DEFAULT_MIN_ROUNDING_BASE, + DIRECTIONS, + ESTIMATOR_DIRECT, + OPTIONAL_LABEL_COLUMNS, + OUTCOMES, + PUBLIC_COLUMNS, + REASON_DISCLOSURE_CONTROL, + REASON_THIN_SUPPORT, + SHARE_GAINING, + SHARE_LOSING, + SHARE_NO_CHANGE, + SHARE_OUTCOMES, + SUPPORT_BANDS, + TIERS, + WITHHELD, + PublicTable, +) + +_KEY = ["level", "area_code", "reform", "year", "outcome"] +_VALUES = ["value", "lo", "hi"] +# An interval this many rounding steps wide means the row was not rounded. +_MAX_INTERVAL_STEPS = 100.0 + + +def _shared_base(numbers: Iterable[float]) -> Optional[float]: + """The largest power of ten that divides every number, if any.""" + values = [number for number in numbers if number != 0] + if not values: + return None + for exponent in range(15, -13, -1): + base = 10.0**exponent + steps = [v / base for v in values] + # A base larger than a number does not divide it. + if all(round(s) != 0 and abs(s - round(s)) < 1e-6 for s in steps): + return base + return 0.0 + + +def _check_frame( + frame: pd.DataFrame, minimum_base: Mapping[str, float], label: str +) -> list[str]: + problems: list[str] = [] + columns = set(frame.columns) + leaked = sorted(columns & set(AUDIT_COLUMNS)) + if leaked: + problems.append(f"{label}: private audit columns present: {leaked}") + missing = [column for column in PUBLIC_COLUMNS if column not in columns] + if missing: + problems.append(f"{label}: missing columns: {missing}") + extra = sorted(columns - set(PUBLIC_COLUMNS) - set(OPTIONAL_LABEL_COLUMNS)) + if extra: + problems.append(f"{label}: columns that are not allowed: {extra}") + if missing: + return problems + + if frame.duplicated(subset=_KEY).any(): + problems.append(f"{label}: duplicate rows for one area, reform and outcome") + + shares: dict[tuple, dict[str, float]] = {} + for number, row in enumerate(frame.to_dict("records"), start=2): + where = f"{label}: row {number}" + tier, outcome = row["tier"], row["outcome"] + if tier not in TIERS: + problems.append(f"{where}: unknown tier") + continue + if outcome not in OUTCOMES: + problems.append(f"{where}: unknown outcome") + continue + if row["noise"] not in ("True", "False"): + problems.append(f"{where}: noise must be True or False") + if row["estimator"] != ESTIMATOR_DIRECT: + problems.append(f"{where}: unknown estimator") + if not row["year"].lstrip("-").isdigit(): + problems.append(f"{where}: year is not a whole number") + + if tier == WITHHELD: + filled = [c for c in (*_VALUES, "direction", "support_band") if row[c]] + if filled: + problems.append(f"{where}: a withheld row has entries in {filled}") + if row["reason"] != REASON_DISCLOSURE_CONTROL: + problems.append(f"{where}: a withheld row has the wrong reason") + continue + + expected = REASON_THIN_SUPPORT if tier == APPROXIMATE else "" + if row["reason"] != expected: + problems.append(f"{where}: the reason does not match the tier") + if row["support_band"] not in SUPPORT_BANDS: + problems.append(f"{where}: unknown support band") + direction = row["direction"] + if outcome in CHANGE_OUTCOMES and direction not in DIRECTIONS: + problems.append(f"{where}: unknown direction") + if outcome in SHARE_OUTCOMES and direction: + problems.append(f"{where}: a share has no direction") + try: + value, lo, hi = (float(row[column]) for column in _VALUES) + except ValueError: + problems.append(f"{where}: value, lo and hi must be numbers") + continue + if not all(math.isfinite(number) for number in (value, lo, hi)): + problems.append(f"{where}: value, lo and hi must be finite") + continue + if not lo <= value <= hi: + problems.append(f"{where}: the value is outside its interval") + if row["support_band"] == BAND_UNDER_10 and value != 0: + problems.append(f"{where}: a non-zero value has a support band under 10") + if outcome in SHARE_OUTCOMES: + if lo < 0 or hi > 1: + problems.append(f"{where}: a share is outside 0 to 1") + key = tuple(row[column] for column in _KEY[:-1]) + shares.setdefault(key, {})[outcome] = value + + base = _shared_base((value, lo, hi)) + if base is not None: + if base < minimum_base[outcome] * (1 - 1e-9): + problems.append(f"{where}: finer than the minimum rounding") + elif outcome != SHARE_NO_CHANGE and hi - lo > _MAX_INTERVAL_STEPS * base: + # No change is derived from the two rounded shares, so its + # value can be finer than its own interval. + problems.append(f"{where}: not rounded to its interval") + + for key, parts in shares.items(): + if len(parts) == 3: + total = parts[SHARE_GAINING] + parts[SHARE_LOSING] + parts[SHARE_NO_CHANGE] + if abs(total - 1.0) > 1e-9: + problems.append(f"{label}: the shares of {key} do not add to one") + elif SHARE_NO_CHANGE in parts: + problems.append(f"{label}: {key} shows no change without both shares") + return problems + + +def _read(path: Path) -> pd.DataFrame: + return pd.read_csv(path, dtype=str, keep_default_na=False) + + +def validate_public( + table: Union[str, Path, pd.DataFrame, PublicTable], + *, + reviewed: Iterable[str] = (), + min_rounding_base: Optional[Mapping[str, float]] = None, +) -> list[str]: + """Return the problems found in a published table, file or directory. + + An empty list means the format checks pass. In a directory every file must + be a table from this package, or be named in ``reviewed`` after a person + checked it. + """ + minimum_base = {**DEFAULT_MIN_ROUNDING_BASE, **(min_rounding_base or {})} + if isinstance(table, PublicTable): + table = table.to_frame() + if isinstance(table, pd.DataFrame): + frame = table.astype(object).where(table.notna(), "").astype(str) + return _check_frame(frame, minimum_base, "table") + + path = Path(table) + if not path.exists(): + return [f"{path}: no such file or directory"] + if path.is_file(): + return _check_frame(_read(path), minimum_base, path.name) + + allowed = set(reviewed) + problems: list[str] = [] + for entry in sorted(path.rglob("*")): + if entry.is_dir(): + continue + name = entry.relative_to(path).as_posix() + if name in allowed or entry.name in allowed: + continue + if entry.suffix.lower() != ".csv": + problems.append(f"{name}: not a disclosure table and not reviewed") + continue + frame = _read(entry) + if not set(PUBLIC_COLUMNS) <= set(frame.columns): + problems.append(f"{name}: not a disclosure table and not reviewed") + continue + problems.extend(_check_frame(frame, minimum_base, name)) + return problems diff --git a/src/policyengine/outputs/__init__.py b/src/policyengine/outputs/__init__.py index ca85f2e9..740227d5 100644 --- a/src/policyengine/outputs/__init__.py +++ b/src/policyengine/outputs/__init__.py @@ -78,6 +78,10 @@ default_uk_geography_asset_strategies, resolve_uk_geography_asset_paths, ) +from policyengine.outputs.uk_geography_disclosure import ( + UKGeographyDisclosure, + compute_uk_geography_disclosure, +) __all__ = [ "Output", @@ -129,6 +133,8 @@ "compute_uk_constituency_impacts", "LocalAuthorityImpact", "compute_uk_local_authority_impacts", + "UKGeographyDisclosure", + "compute_uk_geography_disclosure", "CONSTITUENCY_ASSET_SPEC", "LOCAL_AUTHORITY_ASSET_SPEC", "GCSUKGeographyAssetStrategy", diff --git a/src/policyengine/outputs/constituency_impact.py b/src/policyengine/outputs/constituency_impact.py index 2dbd04ee..214664a6 100644 --- a/src/policyengine/outputs/constituency_impact.py +++ b/src/policyengine/outputs/constituency_impact.py @@ -17,6 +17,7 @@ from policyengine.outputs.uk_geography_impact import ( compute_longwise_uk_geography_impacts, resolve_uk_geography_lookup_csv_path, + warn_unprotected_uk_geography_output, ) if TYPE_CHECKING: @@ -43,6 +44,10 @@ class ConstituencyImpact(Output): def run(self) -> None: """Group household output rows and compute per-constituency metrics.""" + warn_unprotected_uk_geography_output("ConstituencyImpact", "constituency") + self._compute() + + def _compute(self) -> None: baseline_hh = self.baseline_simulation.output_dataset.data.household reform_hh = self.reform_simulation.output_dataset.data.household @@ -66,6 +71,9 @@ def compute_uk_constituency_impacts( ) -> ConstituencyImpact: """Compute per-constituency income changes for UK. + Deprecated: the output is not disclosure-controlled. Use + ``compute_uk_geography_disclosure(levels=("constituency",))`` instead. + Args: baseline_simulation: Completed baseline simulation. reform_simulation: Completed reform simulation. @@ -94,5 +102,8 @@ def compute_uk_constituency_impacts( constituency_csv_path=lookup_csv_path, year=year, ) - impact.run() + warn_unprotected_uk_geography_output( + "compute_uk_constituency_impacts", "constituency" + ) + impact._compute() return impact diff --git a/src/policyengine/outputs/local_authority_impact.py b/src/policyengine/outputs/local_authority_impact.py index 40d50e8f..e38b48d1 100644 --- a/src/policyengine/outputs/local_authority_impact.py +++ b/src/policyengine/outputs/local_authority_impact.py @@ -17,6 +17,7 @@ from policyengine.outputs.uk_geography_impact import ( compute_longwise_uk_geography_impacts, resolve_uk_geography_lookup_csv_path, + warn_unprotected_uk_geography_output, ) if TYPE_CHECKING: @@ -43,6 +44,10 @@ class LocalAuthorityImpact(Output): def run(self) -> None: """Group household output rows and compute per-local-authority metrics.""" + warn_unprotected_uk_geography_output("LocalAuthorityImpact", "local_authority") + self._compute() + + def _compute(self) -> None: baseline_hh = self.baseline_simulation.output_dataset.data.household reform_hh = self.reform_simulation.output_dataset.data.household @@ -66,6 +71,9 @@ def compute_uk_local_authority_impacts( ) -> LocalAuthorityImpact: """Compute per-local-authority income changes for UK. + Deprecated: the output is not disclosure-controlled. Use + ``compute_uk_geography_disclosure(levels=("local_authority",))`` instead. + Args: baseline_simulation: Completed baseline simulation. reform_simulation: Completed reform simulation. @@ -94,5 +102,8 @@ def compute_uk_local_authority_impacts( local_authority_csv_path=lookup_csv_path, year=year, ) - impact.run() + warn_unprotected_uk_geography_output( + "compute_uk_local_authority_impacts", "local_authority" + ) + impact._compute() return impact diff --git a/src/policyengine/outputs/uk_geography_disclosure.py b/src/policyengine/outputs/uk_geography_disclosure.py new file mode 100644 index 00000000..48be972c --- /dev/null +++ b/src/policyengine/outputs/uk_geography_disclosure.py @@ -0,0 +1,267 @@ +"""Disclosure-controlled UK geography outputs. + +A thin adapter from two completed simulations to +:mod:`policyengine.disclosure`. It reads the household output, builds one +checked release for the requested levels, and exposes the rows that may be +published. +""" + +from collections.abc import Iterable, Mapping, Sequence +from pathlib import Path +from typing import TYPE_CHECKING, Optional, Union + +import numpy as np +import pandas as pd +from pydantic import ConfigDict, PrivateAttr + +from policyengine.core import Output +from policyengine.disclosure import ( + AreaWeights, + DisclosureSettings, + Finding, + PriorPublic, + PublicTable, + Release, + estimate_change, +) + +if TYPE_CHECKING: + from policyengine.core.simulation import Simulation + +UK_NATIONAL_LEVEL = "national" +UK_NATIONAL_AREA_CODE = "UK" +UK_REGION_LEVEL = "region" +# Household output columns that hold each level's area code, in order of +# preference. +UK_DISCLOSURE_LEVEL_COLUMNS = { + UK_REGION_LEVEL: ("region_code_oa", "region"), + "constituency": ("constituency_code_oa",), + "local_authority": ("la_code_oa",), +} +UK_SOURCE_HOUSEHOLD_ID_COLUMN = "source_household_id" + + +class UKGeographyDisclosure(Output): + """Income changes by UK area that are safe to publish. + + Every number passes the disclosure rules of + :class:`~policyengine.disclosure.DisclosureSettings` and carries a 90% + interval. ``public`` stays empty while a blocking finding remains. + """ + + model_config = ConfigDict(arbitrary_types_allowed=True) + + baseline_simulation: "Simulation" + reform_simulation: "Simulation" + year: int + levels: tuple[str, ...] = (UK_REGION_LEVEL,) + reform: str = "reform" + settings: Optional[DisclosureSettings] = None + income_variable: str = "household_net_income" + + # Results populated by run() + public: Optional[PublicTable] = None + findings: Optional[list[Finding]] = None + + _release: Optional[Release] = PrivateAttr(default=None) + + @property + def release(self) -> Release: + """The checked release, for example to withhold a cell by hand.""" + if self._release is None: + raise ValueError("Call run() first.") + return self._release + + def run( + self, + *, + noise_key: Optional[Union[str, bytes]] = None, + source_ids: Optional[Sequence] = None, + area_weights: Optional[Mapping[str, AreaWeights]] = None, + parent_of: Optional[Mapping[str, Mapping[str, str]]] = None, + acknowledge: Iterable[str] = (), + prior_public: Sequence[PriorPublic] = (), + audit_dir: Optional[Union[str, Path]] = None, + ) -> None: + """Estimate each level, check the release and fill the results.""" + baseline_hh = pd.DataFrame( + self.baseline_simulation.output_dataset.data.household + ) + reform_hh = pd.DataFrame(self.reform_simulation.output_dataset.data.household) + for name, frame in (("Baseline", baseline_hh), ("Reform", reform_hh)): + if self.income_variable not in frame.columns: + raise ValueError( + f"{name} household output must include {self.income_variable}" + ) + if "household_weight" not in baseline_hh.columns: + raise ValueError("Baseline household output must include household_weight") + if len(baseline_hh) != len(reform_hh): + raise ValueError( + "Baseline and reform household outputs must have the same row count" + ) + if ( + "household_id" in baseline_hh.columns + and "household_id" in reform_hh.columns + ): + same = np.array_equal( + baseline_hh["household_id"].to_numpy(), + reform_hh["household_id"].to_numpy(), + ) + if not same: + raise ValueError( + "Baseline and reform household outputs must hold the same " + "households in the same order" + ) + if source_ids is None: + if UK_SOURCE_HOUSEHOLD_ID_COLUMN not in baseline_hh.columns: + raise ValueError( + f"The household output has no {UK_SOURCE_HOUSEHOLD_ID_COLUMN} " + "column, so support cannot be counted by survey household. " + "Pass source_ids= with one id per household row." + ) + source_ids = baseline_hh[UK_SOURCE_HOUSEHOLD_ID_COLUMN].to_numpy() + + weights = baseline_hh["household_weight"].to_numpy(dtype=float) + custom = dict(area_weights or {}) + parents = dict(parent_of or {}) + + def areas_for(level: str) -> Optional[AreaWeights]: + if level in custom: + return custom[level] + if level == UK_NATIONAL_LEVEL: + codes = np.full(len(baseline_hh), UK_NATIONAL_AREA_CODE, dtype=object) + return AreaWeights.from_labels(codes, weights) + columns = UK_DISCLOSURE_LEVEL_COLUMNS.get(level) + if columns is None: + raise ValueError(f"Unknown level {level!r}") + for column in columns: + if column in baseline_hh.columns: + return AreaWeights.from_labels(baseline_hh[column], weights) + return None + + settings = self.settings or DisclosureSettings() + release = Release(settings, prior_public=prior_public, acknowledge=acknowledge) + # The nation and the regions are always estimated: a level that is not + # requested still serves as the parent in the remainder check. + estimated: dict[str, AreaWeights] = {} + for level in dict.fromkeys((UK_NATIONAL_LEVEL, UK_REGION_LEVEL, *self.levels)): + areas = areas_for(level) + if areas is None: + if level in self.levels: + columns = UK_DISCLOSURE_LEVEL_COLUMNS[level] + raise ValueError( + f"The household output has none of {list(columns)} for " + f"level {level!r}" + ) + continue + estimated[level] = areas + release.add( + estimate_change( + level=level, + reform=self.reform, + year=self.year, + areas=areas, + source_ids=source_ids, + baseline=baseline_hh[self.income_variable].to_numpy(dtype=float), + reformed=reform_hh[self.income_variable].to_numpy(dtype=float), + settings=settings, + noise_key=noise_key, + ), + publish=level in self.levels, + ) + + for level, areas in estimated.items(): + if level == UK_NATIONAL_LEVEL: + continue + under_region = ( + level != UK_REGION_LEVEL + and UK_REGION_LEVEL in estimated + and estimated[UK_REGION_LEVEL].one_area_per_record + and (areas.one_area_per_record or level in parents) + ) + if under_region: + release.add_nesting(level, UK_REGION_LEVEL, parents.get(level)) + else: + release.add_nesting(level, UK_NATIONAL_LEVEL) + + self._release = release + self.refresh() + if audit_dir is not None: + release.write_audit(audit_dir) + + def refresh(self) -> Optional[PublicTable]: + """Check the release again, for example after withholding a cell.""" + release = self.release + self.findings = release.check() + blocked = any(finding.blocking for finding in self.findings) + self.public = None if blocked else release.to_public() + return self.public + + +def compute_uk_geography_disclosure( + baseline_simulation: "Simulation", + reform_simulation: "Simulation", + *, + year: int, + levels: Sequence[str] = (UK_REGION_LEVEL,), + reform: str = "reform", + settings: Optional[DisclosureSettings] = None, + income_variable: str = "household_net_income", + noise_key: Optional[Union[str, bytes]] = None, + source_ids: Optional[Sequence] = None, + area_weights: Optional[Mapping[str, AreaWeights]] = None, + parent_of: Optional[Mapping[str, Mapping[str, str]]] = None, + acknowledge: Iterable[str] = (), + prior_public: Sequence[PriorPublic] = (), + audit_dir: Optional[Union[str, Path]] = None, +) -> UKGeographyDisclosure: + """Compute income changes by UK area under disclosure control. + + Args: + baseline_simulation: Completed baseline simulation. + reform_simulation: Completed reform simulation. + year: Year of the results, used as a label. + levels: Levels to publish: ``national``, ``region``, ``constituency`` + or ``local_authority``. + reform: Label of the reform. A package and its components need + distinct labels. + settings: Publication rules. The defaults follow the UK Data Service + licence: at least 10 distinct survey households behind a number. + income_variable: Household output column whose change is measured. + noise_key: Secret key for the noise, needed only when + ``settings.noise`` is on. Else read from + ``POLICYENGINE_DISCLOSURE_NOISE_KEY``. + source_ids: The survey household of each household row. Needed when + the dataset has no ``source_household_id`` column. + area_weights: Weights for a level, in place of the level's area code + column, for example an area x household weight matrix. + parent_of: For a level given as a weight matrix, the region code of + each of its areas. + acknowledge: Findings the caller accepts after review. + prior_public: Earlier publications that count as already public. + audit_dir: Private directory, outside any repository, for the exact + values and counts. + + Returns: + UKGeographyDisclosure with ``public`` and ``findings`` populated. + ``public`` is None while a blocking finding remains. + """ + impact = UKGeographyDisclosure.model_construct( + baseline_simulation=baseline_simulation, + reform_simulation=reform_simulation, + year=int(year), + levels=tuple(levels), + reform=reform, + settings=settings, + income_variable=income_variable, + ) + impact.run( + noise_key=noise_key, + source_ids=source_ids, + area_weights=area_weights, + parent_of=parent_of, + acknowledge=acknowledge, + prior_public=prior_public, + audit_dir=audit_dir, + ) + return impact diff --git a/src/policyengine/outputs/uk_geography_impact.py b/src/policyengine/outputs/uk_geography_impact.py index ff2528ea..2f48f252 100644 --- a/src/policyengine/outputs/uk_geography_impact.py +++ b/src/policyengine/outputs/uk_geography_impact.py @@ -1,5 +1,6 @@ """Longwise UK geography impact helpers.""" +import warnings from pathlib import Path from typing import Optional @@ -12,6 +13,17 @@ ) +def warn_unprotected_uk_geography_output(name: str, level: str) -> None: + """Warn that an output without disclosure control is deprecated.""" + warnings.warn( + f"{name} is deprecated. Its output is not disclosure-controlled, so it " + "must not be published when the data is licensed survey microdata. Use " + f"compute_uk_geography_disclosure(levels=({level!r},)) instead.", + DeprecationWarning, + stacklevel=3, + ) + + def _normalise_code(value) -> str: if value is None or pd.isna(value): return "" diff --git a/src/policyengine/tax_benefit_models/uk/model.py b/src/policyengine/tax_benefit_models/uk/model.py index 4aafcc0a..d0d1ed65 100644 --- a/src/policyengine/tax_benefit_models/uk/model.py +++ b/src/policyengine/tax_benefit_models/uk/model.py @@ -35,6 +35,12 @@ "constituency_code_oa", "la_code_oa", "region_code_oa", + # The stored region, for region-level outputs on datasets that carry no + # ``region_code_oa``. + "region", + # The survey household a row was copied or cloned from. Disclosure + # control counts support by this id, never by ``household_id``. + "source_household_id", ] diff --git a/tests/test_constituency_impact.py b/tests/test_constituency_impact.py index 4b7e3d1d..c0dcdc92 100644 --- a/tests/test_constituency_impact.py +++ b/tests/test_constituency_impact.py @@ -7,10 +7,17 @@ from microdf import MicroDataFrame from policyengine.outputs.constituency_impact import ( + ConstituencyImpact, compute_uk_constituency_impacts, ) from policyengine.outputs.uk_geography_assets import CONSTITUENCY_ASSET_SPEC +# These outputs are deprecated: they are not disclosure-controlled. The tests +# below still pin their behaviour until they are removed. +pytestmark = pytest.mark.filterwarnings( + "ignore:.*not disclosure-controlled.*:DeprecationWarning" +) + def _make_sim(household_data: dict) -> MagicMock: """Create a mock Simulation with household-level data.""" @@ -248,3 +255,21 @@ def test_compute_constituency_impacts_requires_longwise_geography_column(): reform, download_missing_assets=False, ) + + +def test_constituency_impacts_warn_that_they_are_not_disclosure_controlled(): + household = { + "household_net_income": [100.0], + "household_weight": [1.0], + "constituency_code_oa": ["A001"], + } + baseline, reform = _make_sim(household), _make_sim(household) + message = "not disclosure-controlled.*compute_uk_geography_disclosure" + + with pytest.warns(DeprecationWarning, match=message): + compute_uk_constituency_impacts(baseline, reform, download_missing_assets=False) + impact = ConstituencyImpact.model_construct( + baseline_simulation=baseline, reform_simulation=reform + ) + with pytest.warns(DeprecationWarning, match=message): + impact.run() diff --git a/tests/test_disclosure_areas.py b/tests/test_disclosure_areas.py new file mode 100644 index 00000000..883b2d20 --- /dev/null +++ b/tests/test_disclosure_areas.py @@ -0,0 +1,123 @@ +"""Area x record weights in the label and the matrix layout. Synthetic data.""" + +import numpy as np +import pytest + +from policyengine.disclosure import UNASSIGNED, AreaWeights + + +def _labels_and_matrix(): + codes = np.array(["B", "A", "B", "C", "A", "C", "B"], dtype=object) + weights = np.array([10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0]) + order = ["A", "B", "C"] + matrix = np.zeros((3, len(codes))) + for record, (code, weight) in enumerate(zip(codes, weights)): + matrix[order.index(code), record] = weight + return codes, weights, order, matrix + + +def test_labels_and_a_one_nonzero_matrix_are_the_same_data(): + codes, weights, order, matrix = _labels_and_matrix() + values = np.arange(1.0, 8.0) + multipliers = np.random.default_rng(0).poisson(1.0, size=(5, 7)).astype(np.uint8) + + from_labels = AreaWeights.from_labels(codes, weights) + from_matrix = AreaWeights.from_matrix(matrix, order) + + assert from_matrix.one_area_per_record + assert from_labels.area_codes == from_matrix.area_codes == ("A", "B", "C") + assert np.array_equal(from_labels.labels, from_matrix.labels) + assert np.allclose(from_labels.totals(values), from_matrix.totals(values)) + assert np.allclose( + from_labels.totals(values, multipliers), + from_matrix.totals(values, multipliers), + ) + + +def test_totals_equal_a_direct_weighted_sum(): + codes, weights, _, _ = _labels_and_matrix() + values = np.arange(1.0, 8.0) + areas = AreaWeights.from_labels(codes, weights) + + totals = areas.totals(values) + + assert totals.tolist() == [ + float((weights * values)[codes == code].sum()) for code in ("A", "B", "C") + ] + + +def test_dense_matrix_totals_match_the_matrix_product(): + rng = np.random.default_rng(1) + matrix = rng.uniform(0.0, 2.0, size=(4, 30)).astype(np.float32) + values = rng.normal(size=30) + multipliers = rng.poisson(1.0, size=(6, 30)).astype(np.uint8) + + areas = AreaWeights.from_matrix(matrix, ["a", "b", "c", "d"]) + + assert not areas.one_area_per_record + assert np.shares_memory(areas._matrix, matrix) + assert np.allclose(areas.totals(values), matrix @ values, rtol=1e-5) + expected = matrix @ (multipliers * values).T + assert np.allclose(areas.totals(values, multipliers), expected, rtol=1e-5) + with pytest.raises(ValueError, match="single area"): + _ = areas.labels + + +def test_area_records_and_record_weights(): + codes, weights, _, matrix = _labels_and_matrix() + areas = AreaWeights.from_labels(codes, weights) + + records, record_weights = areas.area_records("B") + + assert sorted(records.tolist()) == [0, 2, 6] + assert sorted(record_weights.tolist()) == [10.0, 30.0, 70.0] + assert areas.record_weights(["A", "C"]).tolist() == [0, 20, 0, 40, 50, 60, 0] + dense = AreaWeights.from_matrix(matrix + 1.0, ["A", "B", "C"]) + assert ( + dense.record_weights(["A", "C"]).tolist() + == (matrix[0] + matrix[2] + 2).tolist() + ) + + +def test_blank_codes_stay_as_unassigned_records(): + areas = AreaWeights.from_labels( + ["A", "", None, b"A", " B ", float("nan")], [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] + ) + + assert areas.area_codes == ("A", "B") + assert areas.labels.tolist() == ["A", UNASSIGNED, UNASSIGNED, "A", "B", UNASSIGNED] + assert areas.record_weights([], include_unassigned=True).tolist() == [ + 0.0, + 2.0, + 3.0, + 0.0, + 0.0, + 6.0, + ] + assert areas.totals(np.ones(6)).tolist() == [5.0, 5.0] + + +def test_an_area_with_no_weight_is_not_an_area(): + areas = AreaWeights.from_labels(["A", "Z", "A"], [1.0, 0.0, 2.0]) + + assert areas.area_codes == ("A",) + assert areas.labels.tolist() == ["A", UNASSIGNED, "A"] + + +@pytest.mark.parametrize( + "weights", [[1.0, -1.0], [1.0, float("nan")], [[1.0, 2.0]], [1.0, float("inf")]] +) +def test_bad_weights_raise(weights): + with pytest.raises(ValueError): + AreaWeights.from_labels(["A", "A"], weights) + + +def test_bad_matrices_raise(): + with pytest.raises(ValueError, match="one row per area"): + AreaWeights.from_matrix(np.ones((2, 3)), ["A"]) + with pytest.raises(ValueError, match="unique"): + AreaWeights.from_matrix(np.ones((2, 3)), ["A", "A"]) + with pytest.raises(ValueError, match="non-negative"): + AreaWeights.from_matrix(-np.ones((2, 3)), ["A", "B"]) + with pytest.raises(ValueError, match="one entry per record"): + AreaWeights.from_labels(["A", "A"], [1.0, 1.0]).totals(np.ones(3)) diff --git a/tests/test_disclosure_estimate.py b/tests/test_disclosure_estimate.py new file mode 100644 index 00000000..6c6a21ce --- /dev/null +++ b/tests/test_disclosure_estimate.py @@ -0,0 +1,494 @@ +"""Private estimates: values, support, intervals and rounding. Synthetic data.""" + +import math + +import numpy as np +import pandas as pd +import pytest + +from policyengine.disclosure import ( + NOISE_KEY_ENV, + AreaWeights, + DisclosureKeyError, + DisclosureSettings, + Release, + estimate_change, + validate_public, +) +from policyengine.disclosure import estimate as estimate_module +from policyengine.disclosure.schema import ( + APPROXIMATE, + BAND_30_PLUS, + BAND_UNDER_10, + PRIVATE_BELOW_FLOOR, + PRIVATE_DOMINANT, + PRIVATE_LOW_N_EFF, + PRIVATE_SHARE_SPLIT, + PRIVATE_UNDEFINED, + REPORTED, + WITHHELD, +) +from policyengine.outputs.uk_geography_impact import ( + compute_longwise_uk_geography_impacts, +) + +TEST_KEY = "test-key-for-unit-tests-only" + + +def _population(spec, *, clones=1, weight=100.0, change=200.0, seed=0): + """Build records from ``{area: (households, affected)}``. + + Each household is its own survey household unless ``clones`` is above 1, + in which case each survey household appears that many times in the area. + """ + rng = np.random.default_rng(seed) + codes, deltas = [], [] + for code, (households, affected) in spec.items(): + codes += [code] * households + deltas += [change] * affected + [0.0] * (households - affected) + n = len(codes) + baseline = rng.uniform(20_000.0, 60_000.0, n) + return { + "codes": np.array(codes, dtype=object), + "weights": np.full(n, weight), + "source_ids": np.arange(n) // clones, + "baseline": baseline, + "reformed": baseline + np.array(deltas), + } + + +def _estimate(data, *, settings=None, noise_key=None, areas=None): + return estimate_change( + level="constituency", + reform="reform", + year=2026, + areas=areas or AreaWeights.from_labels(data["codes"], data["weights"]), + source_ids=data["source_ids"], + baseline=data["baseline"], + reformed=data["reformed"], + settings=settings, + noise_key=noise_key, + ) + + +def _cell(estimates, area, outcome): + table = estimates.table + rows = table[(table["area_code"] == area) & (table["outcome"] == outcome)] + assert len(rows) == 1 + return rows.iloc[0] + + +def test_exact_values_match_the_unprotected_computation(): + rng = np.random.default_rng(3) + n = 240 + codes = np.repeat(np.array(["A", "B", "C"], dtype=object), 80) + weights = rng.uniform(50.0, 400.0, n) + baseline = rng.uniform(15_000.0, 80_000.0, n) + change = np.where(rng.random(n) < 0.5, rng.choice([-300.0, 250.0, 900.0], n), 0.0) + data = { + "codes": codes, + "weights": weights, + "source_ids": np.arange(n), + "baseline": baseline, + "reformed": baseline + change, + } + + estimates = _estimate(data) + household = pd.DataFrame( + {"code": codes, "household_weight": weights, "household_net_income": baseline} + ) + expected = compute_longwise_uk_geography_impacts( + baseline_household=household, + reform_household=household.assign(household_net_income=baseline + change), + geography_column="code", + result_key_prefix="area", + ) + + for row in expected: + average = _cell(estimates, row["area_code"], "average_change") + relative = _cell(estimates, row["area_code"], "relative_change") + total = _cell(estimates, row["area_code"], "total_change") + assert average["exact_value"] == pytest.approx( + row["average_household_income_change"] + ) + assert relative["exact_value"] == pytest.approx( + row["relative_household_income_change"] + ) + assert total["exact_value"] == pytest.approx( + row["average_household_income_change"] * row["population"] + ) + + +def test_the_floor_withholds_under_ten_contributing_households(): + estimates = _estimate(_population({"A": (100, 9), "B": (100, 10)})) + + below = _cell(estimates, "A", "average_change") + at_floor = _cell(estimates, "B", "average_change") + + assert (below["tier"], below["private_reason"]) == (WITHHELD, PRIVATE_BELOW_FLOOR) + assert below["support_band"] == BAND_UNDER_10 + assert at_floor["tier"] == APPROXIMATE + assert at_floor["n_sources"] == 10 + + +def test_clones_of_one_survey_household_count_once(): + estimates = _estimate(_population({"A": (40, 40)}, clones=8)) + + cell = _cell(estimates, "A", "average_change") + + assert cell["n_sources"] == 5 + assert cell["tier"] == WITHHELD + + +def test_outcomes_with_one_numerator_share_a_tier(): + estimates = _estimate(_population({"A": (100, 9), "B": (100, 60)})) + + for area in ("A", "B"): + tiers = { + _cell(estimates, area, outcome)["tier"] + for outcome in ("average_change", "relative_change", "total_change") + } + assert len(tiers) == 1 + + +def test_a_cell_with_no_contributing_household_publishes_zero(): + estimates = _estimate(_population({"A": (50, 0), "B": (100, 60)})) + + cell = _cell(estimates, "A", "average_change") + losing = _cell(estimates, "B", "share_losing") + + assert (cell["value"], cell["lo"], cell["hi"]) == (0.0, 0.0, 0.0) + assert (cell["tier"], cell["direction"]) == (REPORTED, "no_change") + assert cell["support_band"] == BAND_UNDER_10 + assert (losing["value"], losing["tier"]) == (0.0, REPORTED) + assert losing["support_band"] == BAND_30_PLUS + + +def test_changes_inside_the_tolerance_count_as_no_change(): + data = _population({"A": (100, 60)}, change=0.5) + + cell = _cell(_estimate(data), "A", "average_change") + + assert cell["n_sources"] == 0 + assert cell["exact_value"] == 0.0 + + +def test_thin_support_is_approximate_and_strong_support_is_reported(): + estimates = _estimate(_population({"A": (200, 12), "B": (200, 90)})) + + assert _cell(estimates, "A", "average_change")["tier"] == APPROXIMATE + assert _cell(estimates, "B", "average_change")["tier"] == REPORTED + assert _cell(estimates, "B", "average_change")["support_band"] == BAND_30_PLUS + + +def test_a_heavy_contributor_is_withheld_only_while_noise_is_off(): + data = _population({"A": (60, 60)}) + data["reformed"][0] += 500_000.0 + + without_noise = _cell(_estimate(data), "A", "average_change") + with_noise = _cell( + _estimate(data, settings=DisclosureSettings(noise=True), noise_key=TEST_KEY), + "A", + "average_change", + ) + + assert without_noise["top_share"] > 0.5 + assert (without_noise["tier"], without_noise["private_reason"]) == ( + WITHHELD, + PRIVATE_DOMINANT, + ) + assert with_noise["tier"] == REPORTED + + +def test_one_heavy_record_in_many_areas_fails_the_effective_sample_gate(): + """A distinct-count floor passes this cell; the n_eff gate must not.""" + n = 300 + matrix = np.ones((3, n)) + matrix[:, 0] = 1_000_000.0 + data = _population({"X": (n, n)}) + areas = AreaWeights.from_matrix(matrix, ["a", "b", "c"]) + + estimates = _estimate(data, areas=areas) + + for area in ("a", "b", "c"): + cell = _cell(estimates, area, "average_change") + assert cell["n_sources"] == n + assert cell["n_eff"] < 2 + assert (cell["tier"], cell["private_reason"]) == (WITHHELD, PRIVATE_LOW_N_EFF) + assert cell["support_band"] == BAND_UNDER_10 + + +def test_an_even_weight_matrix_passes_and_bands_by_effective_sample(): + n = 300 + data = _population({"X": (n, n)}) + areas = AreaWeights.from_matrix(np.ones((2, n)), ["a", "b"]) + + cell = _cell(_estimate(data, areas=areas), "a", "average_change") + + assert cell["tier"] == REPORTED + assert cell["support_band"] == BAND_30_PLUS + + +def test_noise_changes_the_change_but_never_a_denominator(): + data = _population({"A": (400, 300), "B": (400, 200)}, seed=5) + data["reformed"] += np.where(data["reformed"] > data["baseline"], 150.0, 0.0) + plain = _estimate(data) + noised = _estimate( + data, settings=DisclosureSettings(noise=True), noise_key=TEST_KEY + ) + + for area in ("A", "B"): + plain_total = _cell(plain, area, "total_change")["exact_value"] + noised_total = _cell(noised, area, "total_change")["exact_value"] + assert noised_total != plain_total + assert 0.88 <= noised_total / plain_total <= 1.12 + for outcome in ("average_change", "relative_change"): + plain_denominator = plain_total / _cell(plain, area, outcome)["exact_value"] + noised_denominator = ( + noised_total / _cell(noised, area, outcome)["exact_value"] + ) + assert noised_denominator == pytest.approx(plain_denominator) + + +def test_noise_is_about_a_tenth_of_the_sampling_error_of_a_total(): + rng = np.random.default_rng(11) + data = _population({"A": (2_000, 2_000)}) + data["reformed"] = data["baseline"] + rng.uniform(100.0, 900.0, 2_000) + plain = _cell(_estimate(data), "A", "total_change") + noised = _cell( + _estimate(data, settings=DisclosureSettings(noise=True), noise_key=TEST_KEY), + "A", + "total_change", + ) + + shift = (noised["exact_value"] - plain["exact_value"]) / plain["standard_error"] + + assert abs(shift) < 0.4 + assert noised["standard_error"] / plain["standard_error"] == pytest.approx( + 1.0, abs=0.05 + ) + + +def test_noise_without_a_key_fails_before_anything_is_computed(monkeypatch): + monkeypatch.delenv(NOISE_KEY_ENV, raising=False) + data = _population({"A": (50, 30)}) + data["baseline"] = data["baseline"][:10] + + with pytest.raises(DisclosureKeyError): + _estimate(data, settings=DisclosureSettings(noise=True)) + + +def test_noise_reads_the_key_from_the_environment(monkeypatch): + monkeypatch.setenv(NOISE_KEY_ENV, TEST_KEY) + data = _population({"A": (80, 60)}) + settings = DisclosureSettings(noise=True) + + from_environment = _estimate(data, settings=settings) + from_argument = _estimate(data, settings=settings, noise_key=TEST_KEY) + + assert from_environment.noise + pd.testing.assert_frame_equal(from_environment.table, from_argument.table) + + +def test_shares_stay_between_zero_and_one_under_noise(): + data = _population({"A": (80, 80)}) + + estimates = _estimate( + data, settings=DisclosureSettings(noise=True), noise_key=TEST_KEY + ) + + gaining = _cell(estimates, "A", "share_gaining") + unchanged = _cell(estimates, "A", "share_no_change") + assert 0.0 <= gaining["exact_value"] <= 1.0 + assert 0.0 <= gaining["lo"] <= gaining["value"] <= gaining["hi"] <= 1.0 + assert unchanged["exact_value"] >= 0.0 + + +def test_published_shares_add_to_one(): + rng = np.random.default_rng(8) + n = 600 + baseline = rng.uniform(20_000.0, 50_000.0, n) + data = { + "codes": np.repeat(np.array(["A", "B"], dtype=object), n // 2), + "weights": rng.uniform(80.0, 300.0, n), + "source_ids": np.arange(n), + "baseline": baseline, + "reformed": baseline + rng.choice([-400.0, 0.0, 300.0], n), + } + + estimates = _estimate(data) + + for area in ("A", "B"): + parts = [ + _cell(estimates, area, outcome) + for outcome in ("share_gaining", "share_losing", "share_no_change") + ] + assert all(part["tier"] != WITHHELD for part in parts) + assert sum(part["value"] for part in parts) == pytest.approx(1.0, abs=1e-12) + assert sum(part["exact_value"] for part in parts) == pytest.approx(1.0) + + +def test_published_shares_add_to_one_when_their_rounding_differs(): + """A common share and a rare one round on different bases.""" + n = 1_000 + baseline = np.full(n, 30_000.0) + change = np.array([300.0] * 500 + [-300.0] * 20 + [0.0] * 480) + data = { + "codes": np.full(n, "A", dtype=object), + "weights": np.random.default_rng(9).uniform(80.0, 300.0, n), + "source_ids": np.arange(n), + "baseline": baseline, + "reformed": baseline + change, + } + + estimates = _estimate(data) + + gaining = _cell(estimates, "A", "share_gaining") + losing = _cell(estimates, "A", "share_losing") + unchanged = _cell(estimates, "A", "share_no_change") + assert gaining["rounding_base"] > losing["rounding_base"] + assert unchanged["tier"] != WITHHELD + assert gaining["value"] + losing["value"] + unchanged["value"] == pytest.approx( + 1.0, abs=1e-12 + ) + assert unchanged["lo"] <= unchanged["value"] <= unchanged["hi"] + release = Release() + release.add(estimates) + release.check() + assert validate_public(release.to_public()) == [] + + +def test_a_share_needs_both_sides_above_the_floor(): + estimates = _estimate(_population({"A": (100, 95), "B": (100, 50)})) + + few_others = _cell(estimates, "A", "share_gaining") + balanced = _cell(estimates, "B", "share_gaining") + + assert few_others["tier"] == WITHHELD + assert balanced["tier"] != WITHHELD + + +def test_two_shares_are_not_published_when_they_give_a_small_third(): + baseline = np.full(83, 30_000.0) + change = np.array([300.0] * 40 + [-300.0] * 40 + [0.0] * 3) + data = { + "codes": np.full(83, "A", dtype=object), + "weights": np.full(83, 100.0), + "source_ids": np.arange(83), + "baseline": baseline, + "reformed": baseline + change, + } + + estimates = _estimate(data) + + gaining = _cell(estimates, "A", "share_gaining") + losing = _cell(estimates, "A", "share_losing") + unchanged = _cell(estimates, "A", "share_no_change") + assert gaining["tier"] != WITHHELD + assert (losing["tier"], losing["private_reason"]) == (WITHHELD, PRIVATE_SHARE_SPLIT) + assert unchanged["tier"] == WITHHELD + + +def test_a_uniform_change_has_no_sampling_error_in_its_average(): + estimates = _estimate(_population({"A": (120, 120)}, change=100.0)) + + cell = _cell(estimates, "A", "average_change") + + assert cell["standard_error"] == pytest.approx(0.0, abs=1e-9) + assert (cell["value"], cell["lo"], cell["hi"]) == (100.0, 100.0, 100.0) + assert cell["direction"] == "gain" + + +def test_intervals_are_reproducible_and_do_not_depend_on_the_chunk(monkeypatch): + rng = np.random.default_rng(2) + data = _population({"A": (150, 90), "B": (150, 70)}) + data["reformed"] = data["reformed"] + rng.normal(0.0, 80.0, 300) * ( + data["reformed"] != data["baseline"] + ) + + first = _estimate(data).table + again = _estimate(data).table + monkeypatch.setattr(estimate_module, "_REPLICATE_CHUNK", 37) + chunked = _estimate(data).table + other_seed = _estimate(data, settings=DisclosureSettings(bootstrap_seed=1)).table + + pd.testing.assert_frame_equal(first, again) + pd.testing.assert_frame_equal(first, chunked) + assert not first["exact_lo"].equals(other_seed["exact_lo"]) + + +def test_values_are_rounded_to_the_order_of_their_standard_error(): + rng = np.random.default_rng(4) + data = _population({"A": (300, 200), "B": (300, 120)}) + data["reformed"] = data["baseline"] + rng.uniform(50.0, 2_000.0, 600) * ( + data["reformed"] != data["baseline"] + ) + settings = DisclosureSettings() + + table = _estimate(data).table + + for row in table[table["tier"] != WITHHELD].to_dict("records"): + base = row["rounding_base"] + assert base >= settings.min_rounding_base[row["outcome"]] + assert math.log10(base) == pytest.approx(round(math.log10(base))) + for column in ("value", "lo", "hi"): + assert row[column] / base == pytest.approx(round(row[column] / base)) + assert row["lo"] <= row["value"] <= row["hi"] + if row["outcome"] != "share_no_change" and row["standard_error"] > base: + assert base == 10.0 ** math.floor(math.log10(row["standard_error"])) + average = _cell(_estimate(data), "A", "average_change") + assert average["exact_lo"] < average["exact_value"] < average["exact_hi"] + assert average["direction"] == "gain" + + +def test_an_interval_across_zero_is_too_close_to_call(): + rng = np.random.default_rng(6) + data = _population({"A": (200, 0)}) + data["reformed"] = data["baseline"] + rng.choice([-500.0, 500.0], 200) + + cell = _cell(_estimate(data), "A", "average_change") + + assert cell["lo"] < 0 < cell["hi"] + assert cell["direction"] == "too_close_to_call" + + +def test_relative_change_is_withheld_where_it_is_undefined(): + data = _population({"A": (60, 40)}) + data["reformed"] = data["reformed"] - data["baseline"] + data["baseline"] = np.zeros(60) + + estimates = _estimate(data) + + relative = _cell(estimates, "A", "relative_change") + assert (relative["tier"], relative["private_reason"]) == ( + WITHHELD, + PRIVATE_UNDEFINED, + ) + assert pd.isna(relative["value"]) + assert _cell(estimates, "A", "average_change")["tier"] != WITHHELD + + +def test_repr_holds_no_values(): + estimates = _estimate(_population({"A": (50, 30)})) + + assert repr(estimates) == ( + "AreaEstimates(level='constituency', reform='reform', year=2026, " + "areas=1, private)" + ) + + +@pytest.mark.parametrize("field", ["baseline", "reformed", "source_ids"]) +def test_inputs_must_have_one_entry_per_record(field): + data = _population({"A": (50, 30)}) + data[field] = data[field][:10] + + with pytest.raises(ValueError, match="one entry per record"): + _estimate(data) + + +def test_non_finite_inputs_raise(): + data = _population({"A": (50, 30)}) + data["reformed"][3] = np.nan + + with pytest.raises(ValueError, match="finite"): + _estimate(data) diff --git a/tests/test_disclosure_noise.py b/tests/test_disclosure_noise.py new file mode 100644 index 00000000..8930316b --- /dev/null +++ b/tests/test_disclosure_noise.py @@ -0,0 +1,108 @@ +"""Keyed noise factors. All ids and keys here are synthetic.""" + +import numpy as np +import pytest + +from policyengine.disclosure import ( + NOISE_KEY_ENV, + DisclosureKeyError, + noise_factors, + resolve_noise_key, +) +from policyengine.disclosure.noise import canonical_source_ids + +TEST_KEY = b"test-key-for-unit-tests-only" + + +def test_factors_match_known_answers_for_the_test_key(): + """A change here re-draws the noise of every published table.""" + factors = noise_factors([123, "0123", 7, "2024:55"], TEST_KEY) + + assert factors.tolist() == pytest.approx( + [ + 0.88808850233229, + 1.1191625336822915, + 1.1148886386302803, + 1.0900808835347382, + ], + abs=1e-15, + ) + + +def test_factors_persist_across_calls_and_row_orders(): + ids = np.arange(500) + shuffled = np.random.default_rng(0).permutation(ids) + + first = noise_factors(ids, TEST_KEY) + second = noise_factors(shuffled, TEST_KEY) + + assert np.array_equal(first, noise_factors(ids, TEST_KEY)) + assert np.array_equal(first[shuffled], second) + + +def test_factors_stay_in_the_two_bands_with_nothing_near_one(): + factors = noise_factors(np.arange(20_000), TEST_KEY) + + low = (factors >= 0.88) & (factors <= 0.92) + high = (factors >= 1.08) & (factors <= 1.12) + assert (low | high).all() + assert 0.45 < low.mean() < 0.55 + + +def test_noise_distance_sets_the_bands(): + factors = noise_factors(np.arange(2_000), TEST_KEY, distance=(0.2, 0.3)) + + assert (np.abs(factors - 1.0) >= 0.2).all() + assert (np.abs(factors - 1.0) <= 0.3).all() + + +def test_id_types_name_the_same_household(): + factors = noise_factors([123, "123", 123.0, b"123", " 123 ", "0123"], TEST_KEY) + + assert len(set(factors[:5])) == 1 + assert factors[5] != factors[0] + + +def test_integer_and_string_arrays_agree(): + integers = noise_factors(np.array([5, 6, 7], dtype=np.int64), TEST_KEY) + strings = noise_factors(np.array(["5", "6", "7"]), TEST_KEY) + + assert np.array_equal(integers, strings) + + +@pytest.mark.parametrize("bad", [[1, None], [1.0, float("nan")], ["a", " "], [1.5]]) +def test_missing_blank_or_fractional_ids_raise(bad): + with pytest.raises(ValueError): + canonical_source_ids(bad) + + +def test_different_keys_give_different_factors(): + ids = np.arange(200) + + assert not np.array_equal( + noise_factors(ids, TEST_KEY), noise_factors(ids, b"another-key-0123456789") + ) + + +@pytest.mark.parametrize("key", [None, "", " ", b"", "short-key"]) +def test_key_resolution_fails_closed(key, monkeypatch): + monkeypatch.delenv(NOISE_KEY_ENV, raising=False) + + with pytest.raises(DisclosureKeyError): + resolve_noise_key(key) + + +def test_key_comes_from_the_environment_and_the_argument_wins(monkeypatch): + monkeypatch.setenv(NOISE_KEY_ENV, "environment-key-0123456789\n") + + assert resolve_noise_key() == b"environment-key-0123456789" + assert resolve_noise_key("argument-key-0123456789") == b"argument-key-0123456789" + + +def test_key_error_never_shows_the_key(monkeypatch): + monkeypatch.delenv(NOISE_KEY_ENV, raising=False) + + with pytest.raises(DisclosureKeyError) as error: + resolve_noise_key("tooshort-secret") + + assert "tooshort-secret" not in str(error.value) diff --git a/tests/test_disclosure_release.py b/tests/test_disclosure_release.py new file mode 100644 index 00000000..2b970dad --- /dev/null +++ b/tests/test_disclosure_release.py @@ -0,0 +1,512 @@ +"""The release check and the public rows. Synthetic data.""" + +import stat + +import numpy as np +import pandas as pd +import pytest + +from policyengine.disclosure import ( + AreaWeights, + DisclosureError, + DisclosureSettings, + PriorPublic, + Release, + estimate_change, +) +from policyengine.disclosure.schema import ( + ALREADY_PUBLIC, + APPROXIMATE, + AUDIT_COLUMNS, + NOISE_REFERENCE_PUBLIC, + PRIVATE_SECONDARY, + PUBLIC_COLUMNS, + REASON_DISCLOSURE_CONTROL, + REASON_THIN_SUPPORT, + RECOVERABLE_PARENT, + RECOVERABLE_REMAINDER, + REPORTED, + UNCHECKED_OVERLAP, + WITHHELD, +) + +TEST_KEY = "test-key-for-unit-tests-only" +REGION_OF = {"S1": "R1", "S2": "R1", "S3": "R1", "S4": "R2", "S5": "R2", "S6": "R2"} +STRONG = {"S1": 60, "S2": 60, "S3": 60, "S4": 60, "S5": 60, "S6": 60} + + +def _records(affected, *, households=120, change=200.0, loose=0): + """Households by seat, with the first ``affected[seat]`` of each changed. + + ``loose`` adds changed households in region R1 that have no seat. + """ + seats, regions, deltas = [], [], [] + for seat, region in REGION_OF.items(): + seats += [seat] * households + regions += [region] * households + deltas += [change] * affected[seat] + [0.0] * (households - affected[seat]) + seats += [""] * loose + regions += ["R1"] * loose + deltas += [change] * loose + n = len(seats) + baseline = np.random.default_rng(0).uniform(20_000.0, 60_000.0, n) + return { + "seats": np.array(seats, dtype=object), + "regions": np.array(regions, dtype=object), + "weights": np.full(n, 100.0), + "source_ids": np.arange(n), + "baseline": baseline, + "change": np.array(deltas), + } + + +def _estimates(data, level, *, reform="reform", change=None, settings=None, key=None): + codes = data["seats"] if level == "constituency" else data[f"{level}s"] + delta = data["change"] if change is None else change + return estimate_change( + level=level, + reform=reform, + year=2026, + areas=AreaWeights.from_labels(codes, data["weights"]), + source_ids=data["source_ids"], + baseline=data["baseline"], + reformed=data["baseline"] + delta, + settings=settings, + noise_key=key, + ) + + +def _release(affected, *, loose=0, settings=None, key=None, **options): + data = _records(affected, loose=loose) + release = Release(settings, **options) + release.add(_estimates(data, "region", settings=settings, key=key)) + release.add(_estimates(data, "constituency", settings=settings, key=key)) + release.add_nesting("constituency", "region") + return release + + +def _blocking(findings, code): + return [ + finding for finding in findings if finding.blocking and finding.code == code + ] + + +def _public_row(release, area, outcome="average_change", level="constituency"): + frame = release.to_public().to_frame() + rows = frame[ + (frame["level"] == level) + & (frame["area_code"] == area) + & (frame["outcome"] == outcome) + ] + assert len(rows) == 1 + return rows.iloc[0] + + +def test_to_public_needs_a_check_first(): + release = _release(STRONG) + + with pytest.raises(DisclosureError, match="check"): + release.to_public() + + assert release.check() == [] + release.to_public() + release.withhold( + level="constituency", + area_code="S1", + reform="reform", + year=2026, + outcome="average_change", + ) + assert release.findings is None + with pytest.raises(DisclosureError, match="check"): + release.to_public() + + +def test_public_rows_hold_only_public_columns(): + release = _release(STRONG) + release.check() + + frame = release.to_public().to_frame() + + assert list(frame.columns) == list(PUBLIC_COLUMNS) + assert not set(frame.columns) & set(AUDIT_COLUMNS) + assert set(frame["level"]) == {"region", "constituency"} + assert len(frame) == (2 + 6) * 6 + reported = _public_row(release, "S1") + assert reported["tier"] == REPORTED + assert pd.isna(reported["reason"]) + assert reported["lo"] <= reported["value"] <= reported["hi"] + + +def test_approximate_rows_keep_a_value_and_an_interval(): + release = _release({**STRONG, "S1": 14}) + release.check() + + row = _public_row(release, "S1") + + assert (row["tier"], row["reason"]) == (APPROXIMATE, REASON_THIN_SUPPORT) + assert row["lo"] <= row["value"] <= row["hi"] + assert row["direction"] == "gain" + + +def test_one_withheld_seat_in_a_region_is_recoverable_until_a_second_goes(): + release = _release({**STRONG, "S1": 4, "S2": 30}) + + findings = release.check() + + blocking = _blocking(findings, RECOVERABLE_REMAINDER) + assert {finding.family for finding in blocking} == {"change", "share_gaining"} + for finding in blocking: + assert (finding.relation, finding.group) == ("nesting", "R1") + assert finding.members == ["S1"] + assert finding.suggested == "S2" + with pytest.raises(DisclosureError, match="recoverable_remainder"): + release.to_public() + + for outcome in ("average_change", "share_gaining"): + release.withhold( + level="constituency", + area_code="S2", + reform="reform", + year=2026, + outcome=outcome, + ) + assert release.check() == [] + for seat in ("S1", "S2"): + row = _public_row(release, seat) + assert (row["tier"], row["reason"]) == (WITHHELD, REASON_DISCLOSURE_CONTROL) + for column in ("value", "lo", "hi", "direction", "support_band"): + assert pd.isna(row[column]) + assert _public_row(release, "S3")["tier"] == REPORTED + assert _public_row(release, "S2", "share_no_change")["tier"] == WITHHELD + assert _public_row(release, "S2", "share_losing")["tier"] == REPORTED + + +def test_findings_hold_labels_and_never_counts(): + release = _release({**STRONG, "S1": 4, "S2": 30}) + + for finding in release.check(): + assert not any( + character.isdigit() for character in finding.message.replace("R1", "") + ) + + +def test_withheld_seats_that_pass_together_need_no_more_withholding(): + release = _release({**STRONG, "S1": 6, "S2": 7}) + + findings = release.check() + + assert _blocking(findings, RECOVERABLE_REMAINDER) == [] + assert _public_row(release, "S1")["tier"] == WITHHELD + assert _public_row(release, "S2")["tier"] == WITHHELD + assert _public_row(release, "S3")["tier"] == REPORTED + + +def test_a_withheld_parent_gives_no_total_to_subtract_from(): + release = _release({**STRONG, "S1": 4, "S2": 0, "S3": 0}) + + assert release.check() == [] + assert _public_row(release, "R1", level="region")["tier"] == WITHHELD + assert _public_row(release, "S1")["tier"] == WITHHELD + assert _public_row(release, "S2")["value"] == 0.0 + + +def test_the_suggestion_skips_cells_that_add_nothing(): + release = _release({**STRONG, "S1": 4, "S2": 0, "S3": 60}) + + findings = _blocking(release.check(), RECOVERABLE_REMAINDER) + + assert {finding.suggested for finding in findings} == {"S3"} + assert {finding.group_level for finding in findings} == {"region"} + + +def test_noise_never_protects_a_remainder(): + settings = DisclosureSettings(noise=True) + release = _release({**STRONG, "S1": 4}, settings=settings, key=TEST_KEY) + + assert _blocking(release.check(), RECOVERABLE_REMAINDER) + + +def test_records_with_no_seat_are_part_of_the_remainder(): + release = _release(STRONG, loose=5) + + findings = _blocking(release.check(), RECOVERABLE_REMAINDER) + + assert {finding.group for finding in findings} == {"R1"} + assert all(finding.members == [] for finding in findings) + + +def test_a_withheld_parent_with_all_children_published_is_recoverable(): + release = _release(STRONG) + release.withhold( + level="region", + area_code="R1", + reform="reform", + year=2026, + outcome="average_change", + ) + + findings = _blocking(release.check(), RECOVERABLE_PARENT) + + assert [(f.relation, f.group, f.family) for f in findings] == [ + ("nesting", "R1", "change") + ] + + +def test_a_reference_level_is_checked_against_but_not_published(): + data = _records({**STRONG, "S1": 4}) + release = Release() + release.add(_estimates(data, "region"), publish=False) + release.add(_estimates(data, "constituency")) + release.add_nesting("constituency", "region") + + assert _blocking(release.check(), RECOVERABLE_REMAINDER) + release.withhold( + level="constituency", + area_code="S2", + reform="reform", + year=2026, + outcome="average_change", + ) + release.withhold( + level="constituency", + area_code="S2", + reform="reform", + year=2026, + outcome="share_gaining", + ) + release.check() + assert set(release.to_public().to_frame()["level"]) == {"constituency"} + + +def _package_release(): + data = _records(STRONG) + seat_one = data["seats"] == "S1" + broad = data["change"] + narrow = np.where(seat_one & (np.arange(len(seat_one)) % 120 < 4), 500.0, 0.0) + narrow = np.where(~seat_one & (np.arange(len(seat_one)) % 120 < 60), 500.0, narrow) + release = Release() + for reform, change in (("A", broad), ("B", narrow), ("P", broad + narrow)): + release.add(_estimates(data, "constituency", reform=reform, change=change)) + release.add_package("P", ["A", "B"]) + return release + + +def test_a_package_gives_away_a_withheld_component(): + release = _package_release() + + findings = _blocking(release.check(), RECOVERABLE_REMAINDER) + + assert [(f.relation, f.reform, f.group, f.members) for f in findings] == [ + ("package", "P", "S1", ["B"]) + ] + assert findings[0].suggested == "A" + + release.withhold( + level="constituency", + area_code="S1", + reform="P", + year=2026, + outcome="average_change", + ) + assert release.check() == [] + + +def test_a_withheld_package_with_all_components_published_is_recoverable(): + release = _package_release() + for reform in ("A", "P"): + release.withhold( + level="constituency", + area_code="S4", + reform=reform, + year=2026, + outcome="average_change", + ) + release.withhold( + level="constituency", + area_code="S1", + reform="P", + year=2026, + outcome="average_change", + ) + assert release.check() == [] + + release = _package_release() + release.withhold( + level="constituency", + area_code="S4", + reform="P", + year=2026, + outcome="average_change", + ) + findings = _blocking(release.check(), RECOVERABLE_PARENT) + assert [(f.relation, f.group) for f in findings] == [("package", "S4")] + + +def test_relations_need_their_estimates_on_the_same_records(): + data = _records(STRONG) + release = Release() + release.add(_estimates(data, "constituency")) + release.add_nesting("constituency", "region") + with pytest.raises(DisclosureError, match="needs region estimates"): + release.check() + + other = _records(STRONG, loose=3) + release.add(_estimates(other, "region")) + with pytest.raises(DisclosureError, match="same records"): + release.check() + + release = Release() + release.add(_estimates(data, "constituency", reform="P")) + release.add_package("P", ["A"]) + with pytest.raises(DisclosureError, match="component"): + release.check() + + +def test_a_seat_in_two_regions_is_not_a_nesting(): + data = _records(STRONG) + data["regions"][0] = "R2" + release = Release() + release.add(_estimates(data, "region")) + release.add(_estimates(data, "constituency")) + release.add_nesting("constituency", "region") + + with pytest.raises(DisclosureError, match="more than one"): + release.check() + + +def test_levels_published_together_need_a_declared_relation(): + data = _records(STRONG) + data["local_authoritys"] = np.where(data["seats"] < "S3", "L1", "L2").astype(object) + + def build(**options): + release = Release(**options) + release.add(_estimates(data, "constituency")) + release.add(_estimates(data, "local_authority")) + return release + + blocked = build() + findings = blocked.check() + assert [(f.code, f.blocking, f.level) for f in findings] == [ + (UNCHECKED_OVERLAP, True, "constituency|local_authority") + ] + with pytest.raises(DisclosureError, match=UNCHECKED_OVERLAP): + blocked.to_public() + + accepted = build(acknowledge=[UNCHECKED_OVERLAP]) + assert [f.blocking for f in accepted.check()] == [False] + accepted.to_public() + + with pytest.raises(ValueError, match="cannot be acknowledged"): + Release(acknowledge=[RECOVERABLE_REMAINDER]) + + +def test_a_cell_that_is_already_public_is_not_a_secret(): + cells = [ + ("constituency", "S1", "reform", 2026, "average_change"), + ("constituency", "S1", "reform", 2026, "share_gaining"), + ] + release = _release( + {**STRONG, "S1": 4}, prior_public=[PriorPublic(name="2025 table", cells=cells)] + ) + + findings = release.check() + + assert not any(finding.blocking for finding in findings) + assert {(f.code, f.group, tuple(f.members)) for f in findings} == { + (ALREADY_PUBLIC, "2025 table", ("S1",)) + } + assert _public_row(release, "S1")["tier"] == WITHHELD + + +def test_a_prior_with_another_noise_setting_blocks_until_acknowledged(): + cells = [("constituency", "S3", "reform", 2026, "average_change")] + prior = PriorPublic(name="noised table", noise=True, cells=cells) + + blocked = _release(STRONG, prior_public=[prior]) + findings = blocked.check() + assert [(f.code, f.blocking, f.members) for f in findings] == [ + (NOISE_REFERENCE_PUBLIC, True, ["S3"]) + ] + + accepted = _release( + STRONG, prior_public=[prior], acknowledge=[NOISE_REFERENCE_PUBLIC] + ) + assert [f.blocking for f in accepted.check()] == [False] + + +def test_a_weight_matrix_child_needs_parent_of_under_several_parents(): + data = _records(STRONG) + n = len(data["seats"]) + matrix = np.ones((2, n)) + seats = estimate_change( + level="constituency", + reform="reform", + year=2026, + areas=AreaWeights.from_matrix(matrix, ["M1", "M2"]), + source_ids=data["source_ids"], + baseline=data["baseline"], + reformed=data["baseline"] + data["change"], + ) + + release = Release() + release.add(_estimates(data, "region")) + release.add(seats) + release.add_nesting("constituency", "region") + with pytest.raises(DisclosureError, match="parent_of"): + release.check() + + release = Release() + release.add(_estimates(data, "region")) + release.add(seats) + release.add_nesting("constituency", "region", {"M1": "R1", "M2": "R2"}) + assert release.check() == [] + + +def test_the_release_refuses_mixed_settings_and_duplicates(): + data = _records(STRONG) + release = Release() + release.add(_estimates(data, "constituency")) + + with pytest.raises(DisclosureError, match="already holds"): + release.add(_estimates(data, "constituency")) + with pytest.raises(DisclosureError, match="other settings"): + release.add( + _estimates(data, "region", settings=DisclosureSettings(min_sources=30)) + ) + + +def test_withhold_takes_a_cell_of_the_release_and_a_family_outcome(): + release = _release(STRONG) + cell = {"level": "constituency", "reform": "reform", "year": 2026} + + with pytest.raises(ValueError, match="share_gaining"): + release.withhold(area_code="S1", outcome="share_no_change", **cell) + with pytest.raises(DisclosureError, match="no cell"): + release.withhold(area_code="S9", outcome="average_change", **cell) + + +def test_the_audit_is_private_and_stays_out_of_repositories(tmp_path): + release = _release({**STRONG, "S1": 4}) + release.withhold( + level="constituency", + area_code="S2", + reform="reform", + year=2026, + outcome="average_change", + ) + repository = tmp_path / "repo" + (repository / ".git").mkdir(parents=True) + + with pytest.raises(DisclosureError, match="git repository"): + release.write_audit(repository / "private" / "audit") + + path = release.write_audit(tmp_path / "private") + audit = pd.read_csv(path) + assert stat.S_IMODE(path.stat().st_mode) == 0o600 + assert set(AUDIT_COLUMNS) <= set(audit.columns) + manual = audit[ + (audit["area_code"] == "S2") & (audit["outcome"] == "average_change") + ].iloc[0] + assert (manual["tier"], manual["private_reason"]) == (WITHHELD, PRIVATE_SECONDARY) + assert not list(repository.rglob("*.csv")) diff --git a/tests/test_disclosure_schema.py b/tests/test_disclosure_schema.py new file mode 100644 index 00000000..782dc3ff --- /dev/null +++ b/tests/test_disclosure_schema.py @@ -0,0 +1,154 @@ +"""Settings and public models of the disclosure layer.""" + +import ast +from pathlib import Path + +import pandas as pd +import pytest +from pydantic import ValidationError + +from policyengine.disclosure import ( + DisclosureSettings, + PriorPublic, + PublicRow, + PublicTable, +) +from policyengine.disclosure.schema import ( + BAND_10_TO_29, + PUBLIC_COLUMNS, + REASON_DISCLOSURE_CONTROL, + REPORTED, + WITHHELD, +) + +PACKAGE = Path(__file__).parents[1] / "src" / "policyengine" / "disclosure" + + +@pytest.mark.parametrize( + "loosened", + [ + {"min_sources": 9}, + {"max_top_share": 0.6}, + {"min_n_eff": 9.9}, + {"change_tolerance": 0.0}, + {"n_replicates": 50}, + ], +) +def test_settings_cannot_be_loosened(loosened): + with pytest.raises(ValidationError): + DisclosureSettings(**loosened) + + +def test_settings_can_be_raised_for_special_conditions(): + settings = DisclosureSettings(min_sources=30, min_n_eff=30.0) + + assert settings.min_sources == 30 + + +def test_settings_are_frozen_and_reject_unknown_fields(): + settings = DisclosureSettings() + + with pytest.raises(ValidationError): + settings.min_sources = 30 + with pytest.raises(ValidationError): + DisclosureSettings(noise_key="secret") + + +def test_rounding_bases_must_be_powers_of_ten_for_known_outcomes(): + merged = DisclosureSettings(min_rounding_base={"average_change": 10.0}) + + assert merged.min_rounding_base["average_change"] == 10.0 + assert merged.min_rounding_base["share_gaining"] == 0.001 + with pytest.raises(ValidationError): + DisclosureSettings(min_rounding_base={"average_change": 5.0}) + with pytest.raises(ValidationError): + DisclosureSettings(min_rounding_base={"median_change": 1.0}) + + +def test_noise_distance_must_keep_clear_of_one(): + with pytest.raises(ValidationError): + DisclosureSettings(noise_distance=(0.0, 0.1)) + with pytest.raises(ValidationError): + DisclosureSettings(noise_distance=(0.2, 0.1)) + + +def test_public_table_round_trips_through_json_with_a_withheld_row(): + table = PublicTable( + rows=[ + PublicRow( + level="region", + area_code="A", + reform="r", + year=2026, + outcome="average_change", + value=120.0, + lo=100.0, + hi=140.0, + direction="gain", + tier=REPORTED, + support_band=BAND_10_TO_29, + ), + PublicRow( + level="region", + area_code="B", + reform="r", + year=2026, + outcome="average_change", + tier=WITHHELD, + reason=REASON_DISCLOSURE_CONTROL, + ), + ] + ) + + restored = PublicTable.model_validate_json(table.model_dump_json()) + + assert restored == table + assert list(table.to_frame().columns) == list(PUBLIC_COLUMNS) + assert pd.isna(table.to_frame().loc[1, "value"]) + + +def test_middle_support_band_uses_an_en_dash(): + assert BAND_10_TO_29 == "10–29" + + +def test_prior_public_reads_cells_from_a_frame(): + frame = pd.DataFrame( + { + "level": ["constituency"], + "area_code": ["C1"], + "reform": ["r"], + "year": [2026], + "outcome": ["average_change"], + "value": [1.0], + } + ) + + prior = PriorPublic.from_frame("2025 dashboard", frame) + + assert prior.cells == [("constituency", "C1", "r", 2026, "average_change")] + with pytest.raises(ValueError, match="lacks columns"): + PriorPublic.from_frame("bad", frame.drop(columns=["outcome"])) + + +def test_package_imports_nothing_from_the_rest_of_the_library(): + """The package must stay liftable into its own distribution.""" + for path in sorted(PACKAGE.glob("*.py")): + for node in ast.walk(ast.parse(path.read_text())): + if isinstance(node, ast.Import): + names = [alias.name for alias in node.names] + elif isinstance(node, ast.ImportFrom): + # Relative imports stay inside the package. + names = [node.module or ""] if node.level == 0 else [] + else: + continue + for name in names: + assert not name.startswith("policyengine"), (path.name, name) + assert not name.startswith("scipy"), (path.name, name) + + +def test_package_never_calls_the_builtin_hash(): + """``hash()`` is salted per process, so it cannot give stable factors.""" + for path in sorted(PACKAGE.glob("*.py")): + for node in ast.walk(ast.parse(path.read_text())): + if isinstance(node, ast.Call) and isinstance(node.func, ast.Name): + assert node.func.id != "hash", path.name diff --git a/tests/test_disclosure_support.py b/tests/test_disclosure_support.py new file mode 100644 index 00000000..8bbb4e19 --- /dev/null +++ b/tests/test_disclosure_support.py @@ -0,0 +1,157 @@ +"""Support counts and the primary rules. Synthetic data.""" + +import numpy as np +import pytest + +from policyengine.disclosure import DisclosureSettings +from policyengine.disclosure.schema import ( + APPROXIMATE, + BAND_10_TO_29, + BAND_30_PLUS, + BAND_UNDER_10, + PRIVATE_BELOW_FLOOR, + PRIVATE_DOMINANT, + PRIVATE_LOW_N_EFF, + REPORTED, + WITHHELD, +) +from policyengine.disclosure.support import ( + Support, + cell_tier, + family_supports, + primary_tier, + support, + support_band, +) + +SETTINGS = DisclosureSettings() +LABELS = {"one_area_per_record": True, "noise": False} + + +def test_support_counts_clones_of_one_source_once(): + weights = np.full(20, 100.0) + contributions = np.full(20, 50.0) + source_index = np.repeat(np.arange(5), 4) + + cell = support(weights, contributions, source_index) + + assert cell.n_sources == 5 + assert cell.n_eff == pytest.approx(5.0) + assert cell.top_share == pytest.approx(0.2) + + +def test_support_ignores_records_with_no_change_or_no_weight(): + weights = np.array([100.0, 100.0, 0.0, 100.0]) + contributions = np.array([50.0, 0.0, 50.0, -20.0]) + source_index = np.arange(4) + + cell = support(weights, contributions, source_index) + + assert cell.n_sources == 2 + assert cell.top_share == pytest.approx(50.0 / 70.0) + + +def test_support_of_nothing_is_empty(): + cell = support(np.ones(3), np.zeros(3), np.arange(3)) + + assert cell == Support(0, 0.0, 0.0) + + +def test_effective_sample_falls_with_unequal_weights(): + weights = np.array([1000.0] + [1.0] * 99) + cell = support(weights, np.ones(100), np.arange(100)) + + assert cell.n_sources == 100 + assert cell.n_eff == pytest.approx((1099.0**2) / (1000.0**2 + 99.0)) + assert cell.n_eff < 2 + + +def test_a_share_has_two_sides(): + weights = np.full(50, 10.0) + gaining = np.concatenate([np.ones(12), np.zeros(38)]) + + gainers, others = family_supports(weights, gaining, np.arange(50), "share_gaining") + (change,) = family_supports(weights, gaining, np.arange(50), "change") + + assert (gainers.n_sources, others.n_sources) == (12, 38) + assert change.n_sources == 12 + + +@pytest.mark.parametrize( + ("n_sources", "tier", "reason"), + [ + (0, REPORTED, None), + (9, WITHHELD, PRIVATE_BELOW_FLOOR), + (10, APPROXIMATE, "thin_support"), + (29, APPROXIMATE, "thin_support"), + (30, REPORTED, None), + ], +) +def test_the_count_floor_and_the_thin_tier_on_label_data(n_sources, tier, reason): + cell = Support(n_sources, float(n_sources), 0.1) + + assert primary_tier(cell, SETTINGS, **LABELS) == (tier, reason) + + +def test_effective_sample_withholds_only_on_weight_matrix_data(): + cell = Support(500, 4.0, 0.3) + + on_labels = primary_tier(cell, SETTINGS, one_area_per_record=True, noise=False) + on_matrix = primary_tier(cell, SETTINGS, one_area_per_record=False, noise=False) + + assert on_labels == (APPROXIMATE, "thin_support") + assert on_matrix == (WITHHELD, PRIVATE_LOW_N_EFF) + + +def test_dominance_withholds_only_while_noise_is_off(): + cell = Support(40, 35.0, 0.6) + + without_noise = primary_tier(cell, SETTINGS, one_area_per_record=True, noise=False) + with_noise = primary_tier(cell, SETTINGS, one_area_per_record=True, noise=True) + + assert without_noise == (WITHHELD, PRIVATE_DOMINANT) + assert with_noise == (REPORTED, None) + + +def test_the_count_rule_comes_before_dominance(): + cell = Support(3, 3.0, 0.9) + + assert primary_tier(cell, SETTINGS, **LABELS) == (WITHHELD, PRIVATE_BELOW_FLOOR) + + +def test_the_floor_can_be_raised_to_thirty(): + settings = DisclosureSettings(min_sources=30) + + tier, _ = primary_tier(Support(20, 20.0, 0.1), settings, **LABELS) + + assert tier == WITHHELD + + +def test_the_worst_side_of_a_cell_decides(): + passing, thin, small, empty = ( + Support(40, 40.0, 0.1), + Support(15, 15.0, 0.1), + Support(4, 4.0, 0.4), + Support(0, 0.0, 0.0), + ) + + assert cell_tier((passing, empty), SETTINGS, **LABELS) == (REPORTED, None) + assert cell_tier((passing, thin), SETTINGS, **LABELS)[0] == APPROXIMATE + assert cell_tier((thin, small), SETTINGS, **LABELS)[0] == WITHHELD + + +@pytest.mark.parametrize( + ("n_sources", "band"), + [(0, BAND_UNDER_10), (9, BAND_UNDER_10), (10, BAND_10_TO_29), (29, BAND_10_TO_29)] + + [(30, BAND_30_PLUS), (5_000, BAND_30_PLUS)], +) +def test_support_band_never_shows_a_count(n_sources, band): + cell = Support(n_sources, 1_000.0, 0.1) + + assert support_band(cell, one_area_per_record=True) == band + + +def test_support_band_uses_the_effective_sample_on_weight_matrix_data(): + cell = Support(5_000, 12.0, 0.1) + + assert support_band(cell, one_area_per_record=False) == BAND_10_TO_29 diff --git a/tests/test_disclosure_validate.py b/tests/test_disclosure_validate.py new file mode 100644 index 00000000..a3cac99c --- /dev/null +++ b/tests/test_disclosure_validate.py @@ -0,0 +1,204 @@ +"""Format checks on published tables. Synthetic data.""" + +import numpy as np +import pandas as pd +import pytest + +from policyengine.disclosure import ( + AreaWeights, + Release, + estimate_change, + validate_public, +) +from policyengine.disclosure.__main__ import main + + +def _public_table(): + """A release with reported, approximate, withheld and zero rows.""" + rng = np.random.default_rng(0) + affected = {"A": 4, "B": 14, "C": 70, "D": 0, "E": 6} + codes, deltas = [], [] + for code, count in affected.items(): + codes += [code] * 150 + deltas += list(rng.uniform(100.0, 900.0, count)) + [0.0] * (150 - count) + n = len(codes) + baseline = rng.uniform(20_000.0, 60_000.0, n) + release = Release() + release.add( + estimate_change( + level="constituency", + reform="reform", + year=2026, + areas=AreaWeights.from_labels(codes, np.full(n, 100.0)), + source_ids=np.arange(n), + baseline=baseline, + reformed=baseline + np.array(deltas), + ) + ) + release.check() + return release.to_public() + + +@pytest.fixture(scope="module") +def table(): + return _public_table() + + +@pytest.fixture +def frame(table): + return pd.read_csv( + pd.io.common.StringIO(table.to_frame().to_csv(index=False)), + dtype=str, + keep_default_na=False, + ) + + +def _first(frame, **where): + mask = np.ones(len(frame), dtype=bool) + for column, value in where.items(): + mask &= (frame[column] == value).to_numpy() + return frame.index[mask][0] + + +def test_release_output_passes_as_a_table_a_frame_and_a_file(table, tmp_path): + path = table.to_csv(tmp_path / "constituency.csv") + + assert set(table.to_frame()["tier"]) == {"reported", "approximate", "withheld"} + assert validate_public(table) == [] + assert validate_public(table.to_frame()) == [] + assert validate_public(path) == [] + assert validate_public(tmp_path) == [] + + +def test_a_private_audit_column_is_rejected(frame): + frame["n_sources"] = "12" + + problems = validate_public(frame) + + assert any("private audit columns" in problem for problem in problems) + + +def test_an_unknown_column_is_rejected_but_an_area_name_is_allowed(frame): + frame["area_name"] = "Somewhere" + assert validate_public(frame) == [] + + frame["household_count"] = "100" + assert any("not allowed" in problem for problem in validate_public(frame)) + + +def test_a_missing_column_is_rejected(frame): + problems = validate_public(frame.drop(columns=["tier"])) + + assert any("missing columns" in problem for problem in problems) + + +def test_an_unknown_tier_is_rejected(frame): + frame.loc[0, "tier"] = "fallback" + + assert any("unknown tier" in problem for problem in validate_public(frame)) + + +@pytest.mark.parametrize("column", ["value", "lo", "hi", "direction", "support_band"]) +def test_a_withheld_row_must_be_empty(frame, column): + row = _first(frame, tier="withheld", outcome="average_change") + frame.loc[row, column] = "10–29" if column == "support_band" else "5" + + assert any("withheld row has entries" in p for p in validate_public(frame)) + + +def test_a_withheld_row_must_not_say_why_in_detail(frame): + row = _first(frame, tier="withheld") + frame.loc[row, "reason"] = "secondary" + + assert any("wrong reason" in problem for problem in validate_public(frame)) + + +def test_a_band_with_a_hyphen_is_rejected(frame): + row = _first(frame, tier="approximate") + frame.loc[row, "support_band"] = "10-29" + + assert any("unknown support band" in p for p in validate_public(frame)) + + +def test_an_unrounded_value_is_rejected(frame): + row = _first(frame, tier="reported", outcome="average_change", area_code="C") + frame.loc[row, ["value", "lo", "hi"]] = ["437.6213", "420.1182", "455.9921"] + + assert any("rounding" in p or "rounded" in p for p in validate_public(frame)) + + +def test_a_value_outside_its_interval_is_rejected(frame): + row = _first(frame, tier="reported", outcome="average_change", area_code="C") + frame.loc[row, "lo"] = str(float(frame.loc[row, "hi"]) + 100) + + assert any("outside its interval" in p for p in validate_public(frame)) + + +def test_a_duplicate_row_is_rejected(frame): + doubled = pd.concat([frame, frame.iloc[[0]]], ignore_index=True) + + assert any("duplicate" in problem for problem in validate_public(doubled)) + + +def test_a_non_zero_value_cannot_have_a_band_under_ten(frame): + row = _first(frame, tier="reported", outcome="average_change", area_code="C") + frame.loc[row, "support_band"] = "<10" + + assert any("under 10" in problem for problem in validate_public(frame)) + zero = _first(frame, outcome="average_change", area_code="D") + assert frame.loc[zero, "support_band"] == "<10" + assert frame.loc[zero, "value"] == "0.0" + + +def test_shares_must_add_to_one_and_stay_in_range(frame): + row = _first(frame, outcome="share_no_change", area_code="C") + frame.loc[row, "value"] = "0.5" + frame.loc[row, ["lo", "hi"]] = ["0.4", "0.6"] + + assert any("do not add to one" in p for p in validate_public(frame)) + + frame.loc[row, "hi"] = "1.4" + assert any("outside 0 to 1" in problem for problem in validate_public(frame)) + + +def test_no_change_needs_both_shares(frame): + row = _first(frame, outcome="share_gaining", area_code="C") + frame.loc[row, ["value", "lo", "hi", "support_band"]] = "" + frame.loc[row, ["tier", "reason"]] = ["withheld", "disclosure_control"] + + assert any("without both shares" in p for p in validate_public(frame)) + + +def test_a_directory_must_hold_only_checked_files(table, tmp_path): + table.to_csv(tmp_path / "constituency.csv") + pd.DataFrame({"constituency_code": ["A"], "average_gain": [12.34]}).to_csv( + tmp_path / "old_seat_table.csv", index=False + ) + (tmp_path / "notes.json").write_text("{}") + + problems = validate_public(tmp_path) + + assert problems == [ + "notes.json: not a disclosure table and not reviewed", + "old_seat_table.csv: not a disclosure table and not reviewed", + ] + assert ( + validate_public(tmp_path, reviewed=["notes.json", "old_seat_table.csv"]) == [] + ) + assert validate_public(tmp_path / "missing") != [] + + +def test_the_command_line_exits_with_the_number_of_problems(table, tmp_path, capsys): + good = tmp_path / "good" + good.mkdir() + table.to_csv(good / "constituency.csv") + bad = tmp_path / "bad" + bad.mkdir() + (bad / "seat_values.csv").write_text("constituency_code,average_gain\nA,12.34\n") + + assert main(["validate", str(good)]) == 0 + assert main(["validate", str(good), str(bad)]) == 1 + assert main(["validate", str(bad), "--reviewed", "seat_values.csv"]) == 0 + output = capsys.readouterr().out + assert "1 problem(s) found." in output + assert "12.34" not in output diff --git a/tests/test_extra_variables.py b/tests/test_extra_variables.py index 219c0602..4912e1f4 100644 --- a/tests/test_extra_variables.py +++ b/tests/test_extra_variables.py @@ -159,6 +159,7 @@ def _uk_fixture_dataset(tmp_path): "constituency_code_oa": ["E14001234"], "la_code_oa": ["E09000001"], "region_code_oa": ["E12000007"], + "source_household_id": [4321], "tenure_type": ["RENT_PRIVATELY"], "rent": [12_000.0], "council_tax": [1_500.0], @@ -209,6 +210,8 @@ def test__uk_extra_variables_appear_on_output_dataset(tmp_path) -> None: assert household["constituency_code_oa"].tolist() == ["E14001234"] assert household["la_code_oa"].tolist() == ["E09000001"] assert household["region_code_oa"].tolist() == ["E12000007"] + assert household["source_household_id"].tolist() == [4321] + assert household["region"].tolist() == ["LONDON"] def test__uk_resolve_entity_variables_raises_on_unknown_variable() -> None: diff --git a/tests/test_local_authority_impact.py b/tests/test_local_authority_impact.py index 8e54a50a..aa32d154 100644 --- a/tests/test_local_authority_impact.py +++ b/tests/test_local_authority_impact.py @@ -7,10 +7,17 @@ from microdf import MicroDataFrame from policyengine.outputs.local_authority_impact import ( + LocalAuthorityImpact, compute_uk_local_authority_impacts, ) from policyengine.outputs.uk_geography_assets import LOCAL_AUTHORITY_ASSET_SPEC +# These outputs are deprecated: they are not disclosure-controlled. The tests +# below still pin their behaviour until they are removed. +pytestmark = pytest.mark.filterwarnings( + "ignore:.*not disclosure-controlled.*:DeprecationWarning" +) + def _make_sim(household_data: dict) -> MagicMock: """Create a mock Simulation with household-level data.""" @@ -250,3 +257,23 @@ def test_compute_local_authority_impacts_requires_longwise_geography_column(): reform, download_missing_assets=False, ) + + +def test_local_authority_impacts_warn_that_they_are_not_disclosure_controlled(): + household = { + "household_net_income": [100.0], + "household_weight": [1.0], + "la_code_oa": ["A001"], + } + baseline, reform = _make_sim(household), _make_sim(household) + message = "not disclosure-controlled.*compute_uk_geography_disclosure" + + with pytest.warns(DeprecationWarning, match=message): + compute_uk_local_authority_impacts( + baseline, reform, download_missing_assets=False + ) + impact = LocalAuthorityImpact.model_construct( + baseline_simulation=baseline, reform_simulation=reform + ) + with pytest.warns(DeprecationWarning, match=message): + impact.run() diff --git a/tests/test_uk_geography_disclosure.py b/tests/test_uk_geography_disclosure.py new file mode 100644 index 00000000..ff89d906 --- /dev/null +++ b/tests/test_uk_geography_disclosure.py @@ -0,0 +1,273 @@ +"""The Simulation adapter for disclosure-controlled UK geography outputs. + +All data here is synthetic. +""" + +from unittest.mock import MagicMock + +import numpy as np +import pandas as pd +import pytest +from microdf import MicroDataFrame + +from policyengine import outputs +from policyengine.disclosure import ( + NOISE_KEY_ENV, + AreaWeights, + DisclosureError, + DisclosureKeyError, + DisclosureSettings, + PublicTable, + validate_public, +) +from policyengine.outputs.uk_geography_disclosure import ( + UKGeographyDisclosure, + compute_uk_geography_disclosure, +) + +TEST_KEY = "test-key-for-unit-tests-only" +SEATS = {"S1": "NORTH_EAST", "S2": "NORTH_EAST", "S3": "LONDON", "S4": "LONDON"} + + +def _make_sim(household_data: dict) -> MagicMock: + """Create a mock Simulation with household-level data.""" + hh_df = MicroDataFrame(pd.DataFrame(household_data), weights="household_weight") + sim = MagicMock() + sim.output_dataset.data.household = hh_df + sim.id = "test-sim" + return sim + + +def _households(affected=None, *, per_seat=150): + affected = affected or dict.fromkeys(SEATS, 80) + rng = np.random.default_rng(0) + seats, regions, deltas = [], [], [] + for seat, region in SEATS.items(): + seats += [seat] * per_seat + regions += [region] * per_seat + deltas += list(rng.uniform(100.0, 900.0, affected[seat])) + deltas += [0.0] * (per_seat - affected[seat]) + n = len(seats) + baseline = { + "household_id": np.arange(n), + "household_weight": np.full(n, 100.0), + "household_net_income": rng.uniform(20_000.0, 60_000.0, n), + "region": regions, + "constituency_code_oa": seats, + "source_household_id": np.arange(n) + 5_000, + } + reform = { + **baseline, + "household_net_income": baseline["household_net_income"] + np.array(deltas), + } + return baseline, reform + + +def _compute(baseline, reform, **options): + options.setdefault("year", 2026) + return compute_uk_geography_disclosure( + _make_sim(baseline), _make_sim(reform), **options + ) + + +def test_regions_are_published_with_intervals_and_pass_the_validator(): + impact = _compute(*_households()) + + frame = impact.public.to_frame() + + assert isinstance(impact, UKGeographyDisclosure) + assert impact.findings == [] + assert set(frame["level"]) == {"region"} + assert set(frame["area_code"]) == {"NORTH_EAST", "LONDON"} + averages = frame[frame["outcome"] == "average_change"] + assert (averages["tier"] == "reported").all() + assert (averages["lo"] <= averages["value"]).all() + assert (averages["value"] <= averages["hi"]).all() + assert validate_public(impact.public) == [] + + +def test_public_rows_round_trip_through_json(): + impact = _compute(*_households(), levels=("national", "region", "constituency")) + + restored = PublicTable.model_validate_json(impact.public.model_dump_json()) + + assert restored == impact.public + assert set(impact.public.to_frame()["level"]) == { + "national", + "region", + "constituency", + } + assert set(impact.public.to_frame().query("level == 'national'")["area_code"]) == { + "UK" + } + + +def test_region_codes_are_preferred_over_the_stored_region(): + baseline, reform = _households() + codes = [ + "E12000001" if r == "NORTH_EAST" else "E12000007" for r in baseline["region"] + ] + baseline["region_code_oa"] = codes + reform["region_code_oa"] = codes + + impact = _compute(baseline, reform) + + assert set(impact.public.to_frame()["area_code"]) == {"E12000001", "E12000007"} + + +def test_a_dataset_without_a_source_id_is_refused_unless_ids_are_passed(): + baseline, reform = _households() + ids = baseline.pop("source_household_id") + reform.pop("source_household_id") + + with pytest.raises(ValueError, match="source_ids="): + _compute(baseline, reform) + + impact = _compute(baseline, reform, source_ids=ids) + assert impact.public is not None + + +def test_passed_ids_count_clones_once(): + baseline, reform = _households() + clones = np.arange(len(baseline["household_id"])) // 30 + + impact = _compute(baseline, reform, source_ids=clones, levels=("constituency",)) + + assert ( + impact.public is None + or ( + impact.public.to_frame().query("outcome == 'average_change'")["tier"] + == "withheld" + ).all() + ) + + +def test_a_thin_seat_blocks_publication_until_it_is_resolved(): + baseline, reform = _households({"S1": 4, "S2": 40, "S3": 80, "S4": 80}) + + impact = _compute(baseline, reform, levels=("constituency",)) + + assert impact.public is None + blocking = [finding for finding in impact.findings if finding.blocking] + assert {finding.code for finding in blocking} == {"recoverable_remainder"} + assert {finding.suggested for finding in blocking} == {"S2"} + for outcome in ("average_change", "share_gaining"): + impact.release.withhold( + level="constituency", + area_code="S2", + reform="reform", + year=2026, + outcome=outcome, + ) + public = impact.refresh() + assert impact.public is public + assert impact.findings == [] + tiers = public.to_frame().query("outcome == 'average_change'") + assert dict(zip(tiers["area_code"], tiers["tier"])) == { + "S1": "withheld", + "S2": "withheld", + "S3": "reported", + "S4": "reported", + } + + +def test_frames_must_hold_the_same_households_in_the_same_order(): + baseline, reform = _households() + reform["household_id"] = reform["household_id"][::-1] + + with pytest.raises(ValueError, match="same order"): + _compute(baseline, reform) + + baseline, reform = _households() + with pytest.raises(ValueError, match="household_net_income"): + _compute(baseline, {k: v for k, v in reform.items() if "income" not in k}) + + +def test_an_unknown_or_missing_level_raises(): + baseline, reform = _households() + + with pytest.raises(ValueError, match="Unknown level"): + _compute(baseline, reform, levels=("ward",)) + with pytest.raises(ValueError, match="la_code_oa"): + _compute(baseline, reform, levels=("local_authority",)) + + +def test_seats_and_local_authorities_together_need_an_acknowledgement(): + baseline, reform = _households() + authorities = [ + "L1" if s in ("S1", "S3") else "L2" for s in baseline["constituency_code_oa"] + ] + regions = ["R"] * len(authorities) + for frame in (baseline, reform): + frame["la_code_oa"] = authorities + frame["region"] = regions + levels = ("constituency", "local_authority") + + blocked = _compute(baseline, reform, levels=levels) + accepted = _compute( + baseline, reform, levels=levels, acknowledge=["unchecked_overlap"] + ) + + assert blocked.public is None + assert [finding.code for finding in blocked.findings] == ["unchecked_overlap"] + assert accepted.public is not None + + +def test_noise_needs_a_key_and_the_key_never_reaches_the_results(monkeypatch): + monkeypatch.delenv(NOISE_KEY_ENV, raising=False) + baseline, reform = _households() + settings = DisclosureSettings(noise=True) + + with pytest.raises(DisclosureKeyError): + _compute(baseline, reform, settings=settings) + + impact = _compute(baseline, reform, settings=settings, noise_key=TEST_KEY) + plain = _compute(baseline, reform) + + assert impact.public.to_frame()["noise"].all() + assert not plain.public.to_frame()["noise"].any() + assert TEST_KEY not in impact.public.model_dump_json() + assert TEST_KEY not in repr(impact.__dict__) + assert TEST_KEY not in repr(impact.__pydantic_private__) + assert TEST_KEY not in repr(vars(impact.release)) + + +def test_a_weight_matrix_can_replace_the_area_codes(): + baseline, reform = _households() + n = len(baseline["household_id"]) + matrix = AreaWeights.from_matrix(np.full((2, n), 50.0), ["M1", "M2"]) + + impact = _compute( + baseline, + reform, + levels=("constituency",), + area_weights={"constituency": matrix}, + ) + + assert set(impact.public.to_frame()["area_code"]) == {"M1", "M2"} + + +def test_the_audit_goes_to_a_private_directory_outside_repositories(tmp_path): + baseline, reform = _households() + repository = tmp_path / "repo" + (repository / ".git").mkdir(parents=True) + + with pytest.raises(DisclosureError, match="git repository"): + _compute(baseline, reform, audit_dir=repository / "audit") + + _compute(baseline, reform, audit_dir=tmp_path / "audit") + assert (tmp_path / "audit" / "disclosure_audit.csv").exists() + + +def test_the_release_is_available_only_after_a_run(): + impact = UKGeographyDisclosure.model_construct( + baseline_simulation=MagicMock(), reform_simulation=MagicMock(), year=2026 + ) + + with pytest.raises(ValueError, match=r"run\(\) first"): + _ = impact.release + + +def test_the_adapter_is_exported_from_outputs(): + assert outputs.UKGeographyDisclosure is UKGeographyDisclosure + assert outputs.compute_uk_geography_disclosure is compute_uk_geography_disclosure From 051bd3b61c23b2a849f20a2cdf3fefbeccfd28a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mar=C3=ADa=20Juaristi?= <127882282+juaristi22@users.noreply.github.com> Date: Tue, 13 Oct 2026 17:29:02 +0100 Subject: [PATCH 2/3] Address review of the disclosure layer Add the p% rule on signed totals to the dominance check. Keep weight-matrix rules for weights given as a matrix. Treat withheld cells of reference estimates as hidden in the package and nesting checks. Check each area against the areas of another level wholly inside it. Validate settings again before use and bound the change tolerance. Require region parents for matrix levels. Send the deprecation message through the log. Make AreaEstimates private and unpicklable. Give zero cells the band of the households measured. Check directions in validate_public. Co-Authored-By: Claude Opus 5.5 --- docs/disclosure.md | 93 +++++--- src/policyengine/disclosure/areas.py | 25 ++- src/policyengine/disclosure/estimate.py | 76 ++++--- src/policyengine/disclosure/release.py | 201 +++++++++++++++--- src/policyengine/disclosure/schema.py | 18 +- src/policyengine/disclosure/support.py | 37 +++- src/policyengine/disclosure/validate.py | 15 ++ .../outputs/uk_geography_disclosure.py | 28 ++- .../outputs/uk_geography_impact.py | 17 +- tests/test_constituency_impact.py | 6 +- tests/test_disclosure_areas.py | 11 +- tests/test_disclosure_estimate.py | 80 ++++++- tests/test_disclosure_release.py | 86 +++++++- tests/test_disclosure_schema.py | 30 +++ tests/test_disclosure_support.py | 37 ++++ tests/test_disclosure_validate.py | 15 +- tests/test_local_authority_impact.py | 6 +- tests/test_uk_geography_disclosure.py | 10 +- 18 files changed, 659 insertions(+), 132 deletions(-) diff --git a/docs/disclosure.md b/docs/disclosure.md index 2b240c51..cb61fcd2 100644 --- a/docs/disclosure.md +++ b/docs/disclosure.md @@ -49,22 +49,38 @@ no value, interval, direction or support band. **Count floor.** A number must rest on at least 10 distinct survey households that contribute to it. For a change, a household contributes when its income -changes by more than `change_tolerance` (£1 a year). Copies and clones of one -survey household count once, through the dataset's `source_household_id`. For -a share, both sides must pass: the households that gain and those that do -not. +changes by more than `change_tolerance` (£1 a year, which is also the lowest +value the setting takes). Copies and clones of one survey household count +once, through the dataset's `source_household_id`. For a share, both sides +must pass: the households that gain and those that do not. -**Dominance.** While noise is off, a number is withheld when one survey -household carries more than half of it. +**Dominance.** While noise is off, a number is withheld in two cases: -**Effective sample.** Weights that differ a lot make an estimate less precise -than its household count suggests. Where each record sits in one area, a low -effective sample makes a row `approximate`: it keeps its value and interval, -but the interval understates the uncertainty. Where a record carries weight in -several areas (an area × household weight matrix), the household count passes -by construction, so the effective sample must also be at least 10. +- one survey household carries more than half of it, by size; +- what is left of it without its largest household, or without its two + largest, is less than 10% of the largest household's part (the p% rule). + This catches a total whose other parts cancel, and a total that two + households carry between them. + +These rules protect one or two households. They do not stop a number that +three or more households carry almost alone, as long as 10 contribute. -**Zero cells.** A cell with no contributing household is published as zero. +**Effective sample.** Weights that differ a lot make an estimate less precise +than its household count suggests. Where the data gives one area and one +weight per record, a low effective sample makes a row `approximate`: it keeps +its value and interval, but the interval understates the uncertainty. Where +the data is an area × household weight matrix, the household count passes by +construction, so the effective sample must also be at least 10. Weights given +as a matrix keep this rule whatever they hold: one row of a dense matrix has +one weight per record, yet every household of the nation sits in that area. + +**Zero cells.** A cell with no contributing household is published as zero, +with the interval [0, 0] and the direction `no_change`. It means that no +surveyed household in the area changed, not that no household there does. Its +support band counts the households it was measured on. + +**Settings.** The layer validates its settings again before it uses them, so +a copy made with `model_copy(update=...)` cannot lower a floor. **Linked outcomes.** Average, relative and total change share a numerator, so they are published or withheld together. Gaining, losing and no change add to @@ -89,12 +105,20 @@ looks for such remainders: - child areas against their parent area (constituencies in regions, regions in the nation); - the components of a package against the package; +- an area against the areas of another level that lie wholly inside it, when + two levels that do not nest are released together; - earlier publications, passed as `prior_public`. +A remainder gets the rules of a cell. Estimates added with `publish=False` +are a reference: the release does not emit their rows, but it treats the rows +that pass as public elsewhere and the withheld ones as secret, and checks +them like its own. + `compute_uk_geography_disclosure` always checks constituencies and local authorities against their region, and regions against the nation, also when -you publish one level only. A parent row counts as public when it passes the -rules. +you publish one level only. A level given as a weight matrix needs +`parent_of=` with the region of each area; without it the function raises, +because a check against the nation alone nearly always passes. The check reports; it does not withhold cells by itself. A blocking finding names the group and suggests a cell. Withhold it, then check again. When a @@ -116,10 +140,14 @@ impact.refresh() # checks again and fills impact.public once nothing blocks ``` Two findings can be accepted after review, with `acknowledge=[...]`: -`unchecked_overlap` (levels that overlap without nesting, such as -constituencies and local authorities, are published together) and -`noise_reference_public` (an earlier table shows the same cells under another -noise setting). + +- `unchecked_overlap`: levels that overlap without nesting, such as + constituencies and local authorities, are released together. The release + checks each area against the areas wholly inside it. It does not check + other differences between overlapping areas, and on weight-matrix data it + checks none. +- `noise_reference_public`: an earlier table shows the same cells under + another noise setting. ## Noise @@ -136,8 +164,11 @@ floor still applies. change it between releases of one dataset: two noised versions of a number average towards the truth. - **Limits.** Denominators carry no noise. The factor cancels in a ratio of - two numbers that the same household dominates, and it cancels in every - subtraction, so noise never protects a remainder. + two numbers that the same household dominates. A remainder that subtraction + recovers is the noised sum of the withheld cells, so it has the protection + of a noised cell and no more: the count floor applies to it, and noise + stands in for the dominance rule as it does for a cell. An un-noised copy + of the same cell gives the factor away. ## Validating published files @@ -149,14 +180,28 @@ python -m policyengine.disclosure validate public/data --reviewed notes.csv ``` It rejects private audit columns, unknown tiers, a value in a withheld row, -unrounded values and files it does not recognise. It cannot tell whether a -number rested on enough households; only the release that made the file can. +unrounded values, a direction that does not match its interval, and files it +does not recognise. It cannot tell whether a number rested on enough +households; only the release that made the file can. -## The private audit +## Private objects `audit_dir=` writes the exact values, household counts and shares behind every row. The directory must be outside any git repository. +`AreaEstimates`, the object that `estimate_change` returns, holds the same +exact values. With noise on it also holds the noise factor of every household +that gains or loses, and within one dataset those factors do the same job as +the key. Keep it in memory and treat it as you treat the key. It cannot be +pickled or copied. + +## What a reader can still derive + +`average_change` and `relative_change` together give the area's total +baseline income to within a few per cent (household count × average ÷ +relative). That total is a denominator: it carries no noise and is not +checked for dominance. + ## Arrays without a simulation The package works on arrays and imports nothing from the rest of the library: diff --git a/src/policyengine/disclosure/areas.py b/src/policyengine/disclosure/areas.py index 9ebe1152..e7058d9c 100644 --- a/src/policyengine/disclosure/areas.py +++ b/src/policyengine/disclosure/areas.py @@ -3,8 +3,8 @@ Two layouts hold the same thing. In the label layout each record sits in one area with one weight. In the matrix layout a record can carry weight in many areas. The rules read both through the same operations, so they are written -once. A matrix with at most one non-zero weight per record is stored in the -label layout. +once. Which rules apply depends on how the data was given, not on what it +holds: weights given as a matrix keep the weight-matrix rules. """ from collections.abc import Iterable, Sequence @@ -55,8 +55,10 @@ def __init__( label_index: Optional[np.ndarray] = None, weights: Optional[np.ndarray] = None, matrix: Optional[np.ndarray] = None, + from_matrix: bool = False, ) -> None: self.area_codes = tuple(area_codes) + self._from_matrix = from_matrix or matrix is not None self.n_records = int(n_records) self._position = {code: i for i, code in enumerate(self.area_codes)} # The label layout fills the first two; the matrix layout the third. @@ -89,6 +91,7 @@ def _from_codes( codes: np.ndarray, weights: np.ndarray, order: Optional[Sequence[str]] = None, + from_matrix: bool = False, ) -> "AreaWeights": weighted = {code for code in codes[weights > 0] if code != UNASSIGNED} ordered = order if order is not None else sorted(weighted) @@ -101,11 +104,17 @@ def _from_codes( n_records=len(codes), label_index=label_index, weights=weights, + from_matrix=from_matrix, ) @classmethod def from_matrix(cls, matrix: object, area_codes: Sequence[object]) -> "AreaWeights": - """An areas x records weight matrix. The matrix is not copied.""" + """An areas x records weight matrix. The matrix is not copied. + + Data given as a matrix keeps the weight-matrix rules whatever it + holds. One row of a dense matrix has one weight per record, yet every + record of the nation still sits in that area. + """ array = np.asarray(matrix) if array.ndim != 2: raise ValueError("A weight matrix must be two-dimensional.") @@ -119,14 +128,15 @@ def from_matrix(cls, matrix: object, area_codes: Sequence[object]) -> "AreaWeigh if (array < 0).any(): raise ValueError("Weights must be non-negative.") if array.shape[0] and np.count_nonzero(array, axis=0).max(initial=0) <= 1: - # One non-zero per record: the same data as one label per record. + # One non-zero per record: store it as labels, which is faster. + # The rules still treat it as a weight matrix. weights = array.max(axis=0).astype(float) row = array.argmax(axis=0) labels = np.array( [codes[r] if w > 0 else UNASSIGNED for r, w in zip(row, weights)], dtype=object, ) - return cls._from_codes(labels, weights, order=codes) + return cls._from_codes(labels, weights, order=codes, from_matrix=True) return cls(area_codes=codes, n_records=array.shape[1], matrix=array) @property @@ -135,7 +145,8 @@ def n_areas(self) -> int: @property def one_area_per_record(self) -> bool: - return self._matrix is None + """Whether the data was given as one area and one weight per record.""" + return not self._from_matrix @property def labels(self) -> np.ndarray: @@ -202,4 +213,4 @@ def totals( if self.n_areas == 0: return np.zeros((0, multipliers.shape[0])) by_area = multipliers[:, self._order] * weighted[self._order] - return np.add.reduceat(by_area, self._starts, axis=1).T + return np.asarray(np.add.reduceat(by_area, self._starts, axis=1).T) diff --git a/src/policyengine/disclosure/estimate.py b/src/policyengine/disclosure/estimate.py index 1244ff54..958e30fa 100644 --- a/src/policyengine/disclosure/estimate.py +++ b/src/policyengine/disclosure/estimate.py @@ -52,7 +52,13 @@ class AreaEstimates: - """Estimates and support for the areas of one level. Private.""" + """Estimates and support for the areas of one level. Private. + + The object holds exact values, household counts and, with noise on, the + noise factor of every household that gains or loses. Within one dataset + those factors do the same job as the key, so treat the object as you treat + the key: keep it in memory, and let a release write the audit. + """ def __init__( self, @@ -65,19 +71,17 @@ def __init__( settings: DisclosureSettings, noise: bool, contributions: dict[str, np.ndarray], - no_change: np.ndarray, table: pd.DataFrame, ) -> None: self.level = level self.reform = reform self.year = year self.areas = areas - self.source_index = source_index self.settings = settings self.noise = noise - self.contributions = contributions - self.no_change = no_change - self.table = table + self._source_index = source_index + self._contributions = contributions + self._table = table def __repr__(self) -> str: # The table holds exact values and counts, so it stays out of reprs. @@ -86,15 +90,21 @@ def __repr__(self) -> str: f"year={self.year!r}, areas={self.areas.n_areas}, private)" ) + def __getstate__(self) -> dict: + raise TypeError( + "AreaEstimates holds exact values and noise factors, so it cannot " + "be pickled or copied. Use Release.write_audit for a private record." + ) + def family_tiers(self, family: str) -> dict[str, str]: """The primary tier of each area's cell in one family.""" outcome = AVERAGE_CHANGE if family == FAMILY_CHANGE else family - rows = self.table[self.table["outcome"] == outcome] + rows = self._table[self._table["outcome"] == outcome] return dict(zip(rows["area_code"], rows["family_tier"])) def family_sources(self, family: str) -> dict[str, int]: outcome = AVERAGE_CHANGE if family == FAMILY_CHANGE else family - rows = self.table[self.table["outcome"] == outcome] + rows = self._table[self._table["outcome"] == outcome] return dict(zip(rows["area_code"], rows["n_sources"])) @@ -169,6 +179,23 @@ def _replicates( return {outcome: np.concatenate(chunks, axis=1) for outcome, chunks in out.items()} +def _interval( + draws: np.ndarray, exact: np.ndarray, alpha: float +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Percentile bounds and standard error of each area's replicates. + + A skewed set of replicates can leave the estimate outside its percentile + bounds, so the bounds are widened to include it. + """ + lower = np.nanpercentile(draws, 100.0 * alpha, axis=1) + upper = np.nanpercentile(draws, 100.0 * (1.0 - alpha), axis=1) + return ( + np.fmin(lower, exact), + np.fmax(upper, exact), + np.nanstd(draws, axis=1, ddof=1), + ) + + def _direction(value: float, lo: float, hi: float) -> str: if lo > 0: return GAIN @@ -197,7 +224,7 @@ def estimate_change( household net income. ``source_ids`` names the survey household each record was copied or cloned from. """ - settings = settings or DisclosureSettings() + settings = (settings or DisclosureSettings()).revalidated() # Resolve the key first, so a missing key fails before any computation. key = resolve_noise_key(noise_key) if settings.noise else None n_records = areas.n_records @@ -253,13 +280,7 @@ def estimate_change( with warnings.catch_warnings(): warnings.simplefilter("ignore", category=RuntimeWarning) for outcome, draws in replicates.items(): - lower = np.nanpercentile(draws, 100.0 * alpha, axis=1) - upper = np.nanpercentile(draws, 100.0 * (1.0 - alpha), axis=1) - bounds[outcome] = ( - np.fmin(lower, exact[outcome]), - np.fmax(upper, exact[outcome]), - np.nanstd(draws, axis=1, ddof=1), - ) + bounds[outcome] = _interval(draws, exact[outcome], alpha) one_area = areas.one_area_per_record rules = {"one_area_per_record": one_area, "noise": settings.noise} @@ -350,9 +371,15 @@ def estimate_change( tiers[family], reasons[family], binding_support(supports[family]), + # A cell with no contributor rests on the households that + # did not change, as a zero share rests on its other side. + binding_support((*supports[family], unchanged)), one_area, ) ) + third_cell = binding_support( + (*supports[FAMILY_GAINING], *supports[FAMILY_LOSING], unchanged) + ) rows.append( _row( level, @@ -367,13 +394,8 @@ def estimate_change( rounded[SHARE_NO_CHANGE], third, third_reason, - binding_support( - ( - *supports[FAMILY_GAINING], - *supports[FAMILY_LOSING], - unchanged, - ) - ), + third_cell, + third_cell, one_area, ) ) @@ -387,7 +409,6 @@ def estimate_change( settings=settings, noise=settings.noise, contributions=contributions, - no_change=no_change, table=pd.DataFrame(rows), ) @@ -406,6 +427,7 @@ def _row( family_tier: str, family_reason: Optional[str], cell: Support, + band_cell: Support, one_area_per_record: bool, ) -> dict: """One private row: exact and rounded values with the cell's support.""" @@ -442,5 +464,9 @@ def _row( "n_sources": cell.n_sources, "n_eff": cell.n_eff, "top_share": cell.top_share, - "support_band": support_band(cell, one_area_per_record=one_area_per_record), + "rest_share": cell.rest_share, + "support_band": support_band( + cell if cell.n_sources else band_cell, + one_area_per_record=one_area_per_record, + ), } diff --git a/src/policyengine/disclosure/release.py b/src/policyengine/disclosure/release.py index a9136b0f..464db136 100644 --- a/src/policyengine/disclosure/release.py +++ b/src/policyengine/disclosure/release.py @@ -65,7 +65,7 @@ def __init__( prior_public: Sequence[PriorPublic] = (), acknowledge: Iterable[str] = (), ) -> None: - self.settings = settings or DisclosureSettings() + self.settings = (settings or DisclosureSettings()).revalidated() self._acknowledged = frozenset(acknowledge) unknown = sorted(self._acknowledged - set(ACKNOWLEDGEABLE_FINDINGS)) if unknown: @@ -81,9 +81,10 @@ def __init__( def add(self, estimates: AreaEstimates, *, publish: bool = True) -> None: """Add the estimates of one level, reform and year. - With ``publish=False`` the estimates are a reference only: the release - does not emit their rows, and treats the rows that pass the rules as - public elsewhere. + With ``publish=False`` the estimates are a reference: the release does + not emit their rows. It treats the rows that pass the rules as public + elsewhere and the withheld rows as secret, and checks them like its + own. """ if estimates.settings != self.settings: raise DisclosureError( @@ -152,7 +153,7 @@ def check(self) -> list[Finding]: findings: list[Finding] = [] self._check_nesting(tiers, known, findings) self._check_packages(tiers, known, findings) - self._check_overlap(findings) + self._check_overlap(tiers, known, findings) self._check_prior(tiers, findings) self._findings = findings return findings @@ -221,7 +222,7 @@ def write_audit(self, private_dir: Union[str, Path]) -> Path: directory.mkdir(parents=True, exist_ok=True) frames = [] for key in sorted(self._estimates): - frame = self._estimates[key].table.drop(columns=["family_tier"]).copy() + frame = self._estimates[key]._table.drop(columns=["family_tier"]).copy() effective = list(self._effective_rows(key)) frame["tier"] = [tier for _, tier, _ in effective] frame["private_reason"] = [reason for _, _, reason in effective] @@ -251,7 +252,7 @@ def _effective_rows(self, key: _Key) -> Iterable[tuple[dict, str, Optional[str]] for level, reform, year, family, area in self._manual if (level, reform, year) == key } - for record in self._estimates[key].table.to_dict("records"): + for record in self._estimates[key]._table.to_dict("records"): outcome, area = record["outcome"], record["area_code"] families = ( (FAMILY_GAINING, FAMILY_LOSING) @@ -280,12 +281,18 @@ def _known_cells(self) -> set[tuple[str, str, str, int, str]]: def _passes( self, estimates: AreaEstimates, weights: np.ndarray, family: str ) -> tuple[bool, bool]: - """Whether a remainder may be recovered, and whether it is empty.""" + """Whether a remainder may be recovered, and whether it is empty. + + A remainder gets the rules of a cell, no more and no less. The count + floor always applies. With noise on, what subtraction recovers is the + noised sum of the withheld cells, so noise stands in for the dominance + rule here as it does for a cell. + """ records = np.flatnonzero(weights) supports = family_supports( weights[records], - estimates.contributions[family][records], - estimates.source_index[records], + estimates._contributions[family][records], + estimates._source_index[records], family, ) tier, _ = cell_tier( @@ -300,7 +307,7 @@ def _passes( @staticmethod def _require_same_records(first: AreaEstimates, second: AreaEstimates) -> None: same = first.areas.n_records == second.areas.n_records and np.array_equal( - first.source_index, second.source_index + first._source_index, second._source_index ) if not same: raise DisclosureError( @@ -372,7 +379,7 @@ def _check_nesting( for child_level, parent_level, parent_of in self._nestings: for key in sorted(self._estimates): level, reform, year = key - if level != child_level or not self._published[key]: + if level != child_level: continue child = self._estimates[key] parent_key = (parent_level, reform, year) @@ -382,6 +389,9 @@ def _check_nesting( f"Nesting {child_level} in {parent_level} needs " f"{parent_level} estimates for reform {reform!r}, {year}." ) + # A relation between two references is not this release's. + if not (self._published[key] or self._published[parent_key]): + continue self._require_same_records(child, parent) groups = self._groups(child, parent, parent_of) for family in FAMILIES: @@ -414,7 +424,7 @@ def _check_nesting( # There is no public parent to subtract from. The # risk is the reverse: published children that add # up to a parent this release withholds. - if empty and self._published[parent_key]: + if empty: findings.append( Finding( code=RECOVERABLE_PARENT, @@ -432,10 +442,14 @@ def _check_nesting( if not passes: # Withholding a cell with no contributor adds # nothing to the remainder. + # A cell that is public elsewhere cannot be + # withheld here either. shown = [ c for c in children - if child_tiers[c] != WITHHELD and sources[c] > 0 + if self._published[key] + and child_tiers[c] != WITHHELD + and sources[c] > 0 ] findings.append( Finding( @@ -469,7 +483,7 @@ def _check_packages( for package, components in self._packages: for key in sorted(self._estimates): level, reform, year = key - if reform != package or not self._published[key]: + if reform != package: continue whole = self._estimates[key] parts = [] @@ -487,12 +501,17 @@ def _check_packages( "A package and its components need the same areas." ) parts.append((part_key, part)) + # A relation between references only is not this release's. + if not ( + self._published[key] + or any(self._published[part_key] for part_key, _ in parts) + ): + continue for area in whole.areas.area_codes: hidden = [ (part_key, part) for part_key, part in parts - if self._published[part_key] - and tiers[(part_key, family)][area] == WITHHELD + if tiers[(part_key, family)][area] == WITHHELD and (level, area, part.reform, year, family) not in known ] whole_hidden = ( @@ -524,10 +543,10 @@ def _check_packages( elif hidden and not whole_hidden: records, weights = whole.areas.area_records(area) combined = sum( - part.contributions[family][records] for _, part in hidden + part._contributions[family][records] for _, part in hidden ) supports = family_supports( - weights, combined, whole.source_index[records], family + weights, combined, whole._source_index[records], family ) tier, _ = cell_tier( supports, @@ -540,7 +559,8 @@ def _check_packages( shown = [ part for part_key, part in parts - if tiers[(part_key, family)][area] != WITHHELD + if self._published[part_key] + and tiers[(part_key, family)][area] != WITHHELD and part.family_sources(family)[area] > 0 ] findings.append( @@ -569,8 +589,18 @@ def _check_packages( ) ) - def _check_overlap(self, findings: list[Finding]) -> None: - """Levels published together need a declared relation between them.""" + def _check_overlap( + self, + tiers: dict[tuple[_Key, str], dict[str, str]], + known: set[tuple[str, str, str, int, str]], + findings: list[Finding], + ) -> None: + """Levels released together with no declared nesting between them. + + Areas of one level that lie wholly inside an area of the other are + checked like a nesting. Other differences between overlapping areas + are not, and the finding says so. + """ related: dict[str, set[str]] = {} for child_level, parent_level, _ in self._nestings: related.setdefault(child_level, set()).add(parent_level) @@ -588,30 +618,131 @@ def nests(child: str, parent: str) -> bool: return False by_run: dict[tuple[str, int], set[str]] = {} - for (level, reform, year), publish in self._published.items(): - if publish: - by_run.setdefault((reform, year), set()).add(level) - pairs = { - pair - for levels in by_run.values() - for pair in combinations(sorted(levels), 2) - if not nests(*pair) and not nests(*reversed(pair)) - } - for first, second in sorted(pairs): + for level, reform, year in self._estimates: + by_run.setdefault((reform, year), set()).add(level) + pairs: dict[tuple[str, str], list[tuple[str, int]]] = {} + for (reform, year), levels in sorted(by_run.items()): + for first, second in combinations(sorted(levels), 2): + released = ( + self._published[(first, reform, year)] + or self._published[(second, reform, year)] + ) + if released and not nests(first, second) and not nests(second, first): + pairs.setdefault((first, second), []).append((reform, year)) + for (first, second), runs in sorted(pairs.items()): + contained = [ + self._check_containment( + (outer, reform, year), (inner, reform, year), tiers, known, findings + ) + for reform, year in runs + for outer, inner in ((first, second), (second, first)) + ] + scope = ( + "Areas that lie wholly inside an area of the other level are " + "checked. Other differences between overlapping areas are not." + if all(contained) + else "Differences between overlapping areas are not checked." + ) findings.append( Finding( code=UNCHECKED_OVERLAP, blocking=UNCHECKED_OVERLAP not in self._acknowledged, relation="release", message=( - f"{first} and {second} are published together with no " - "declared nesting. Differences between overlapping areas " - "are not checked." + f"{first} and {second} are released together with no " + f"declared nesting. {scope}" ), level=f"{first}|{second}", ) ) + def _check_containment( + self, + outer_key: _Key, + inner_key: _Key, + tiers: dict[tuple[_Key, str], dict[str, str]], + known: set[tuple[str, str, str, int, str]], + findings: list[Finding], + ) -> bool: + """Check an outer area against the inner areas wholly inside it. + + An outer area minus its public inner areas is a remainder. Returns + False when the two levels do not allow the check: they must hold the + same records with one area and one weight per record. + """ + outer, inner = self._estimates[outer_key], self._estimates[inner_key] + outer_level, reform, year = outer_key + inner_level = inner_key[0] + try: + self._require_same_records(outer, inner) + outer_labels, inner_labels = outer.areas.labels, inner.areas.labels + except (DisclosureError, ValueError): + return False + outer_weights = outer.areas.record_weights(outer.areas.area_codes) + inner_weights = inner.areas.record_weights(inner.areas.area_codes) + both = (outer_weights > 0) & (inner_weights > 0) + if not np.array_equal(outer_weights[both], inner_weights[both]): + return False + + pairs = pd.DataFrame({"outer": outer_labels, "inner": inner_labels}) + pairs = pairs[inner_weights > 0].drop_duplicates() + whole = pairs[~pairs["inner"].duplicated(keep=False) & (pairs["outer"] != "")] + inside = whole.groupby("outer")["inner"].apply(sorted).to_dict() + for family in FAMILIES: + outer_tiers = tiers[(outer_key, family)] + inner_tiers = tiers[(inner_key, family)] + sources = inner.family_sources(family) + for code, members in inside.items(): + if ( + outer_tiers.get(code) == WITHHELD + and (outer_level, code, reform, year, family) not in known + ): + continue + public = [ + c + for c in members + if inner_tiers[c] != WITHHELD + or (inner_level, c, reform, year, family) in known + ] + if not public: + continue + weights = outer.areas.record_weights([code]) + weights = np.clip(weights - inner.areas.record_weights(public), 0, None) + passes, _ = self._passes(outer, weights, family) + if passes: + continue + shown = [ + c + for c in public + if self._published[inner_key] + and inner_tiers[c] != WITHHELD + and sources[c] > 0 + ] + findings.append( + Finding( + code=RECOVERABLE_REMAINDER, + blocking=True, + relation="containment", + message=( + f"{outer_level} {code!r} minus the public " + f"{inner_level} cells inside it fails the rules, so " + "subtraction recovers the rest. Withhold one of " + f"those cells, or the {outer_level} cell." + ), + level=inner_level, + reform=reform, + year=year, + family=family, + group=code, + group_level=outer_level, + members=sorted(set(members) - set(public)), + suggested=min( + shown, key=lambda c: (sources[c], c), default=None + ), + ) + ) + return True + def _check_prior( self, tiers: dict[tuple[_Key, str], dict[str, str]], diff --git a/src/policyengine/disclosure/schema.py b/src/policyengine/disclosure/schema.py index 565558b7..380d9d2b 100644 --- a/src/policyengine/disclosure/schema.py +++ b/src/policyengine/disclosure/schema.py @@ -94,6 +94,7 @@ "n_sources", "n_eff", "top_share", + "rest_share", "private_reason", ) @@ -129,11 +130,16 @@ class DisclosureSettings(BaseModel): min_sources: int = Field(default=10, ge=10) max_top_share: float = Field(default=0.5, gt=0.0, le=0.5) + # The p% rule: what is left of a number without its largest household, or + # its two largest, must be at least this share of the largest one's part. + min_rest_share: float = Field(default=0.1, ge=0.1, le=1.0) min_n_eff: float = Field(default=10.0, ge=10.0) approximate_below_n_eff: float = Field(default=30.0, ge=0.0) # A household contributes, gains or loses when its change is larger than - # this. Smaller changes count as no change and are left out of every sum. - change_tolerance: float = Field(default=1.0, gt=0.0) + # this, in the outcome's unit (pounds a year). Smaller changes count as no + # change and are left out of every sum. The bound keeps floating-point + # residue from counting as a contributor. + change_tolerance: float = Field(default=1.0, ge=1.0) noise: bool = False # Distance of a noise factor from 1: uniform on this range, either side. noise_distance: tuple[float, float] = (0.08, 0.12) @@ -144,6 +150,14 @@ class DisclosureSettings(BaseModel): default_factory=lambda: dict(DEFAULT_MIN_ROUNDING_BASE) ) + def revalidated(self) -> "DisclosureSettings": + """These settings, validated again. + + ``model_copy(update=...)`` and ``model_construct`` skip validation, so + the layer validates the settings it is given before it uses them. + """ + return type(self).model_validate(self.model_dump()) + @field_validator("noise_distance") @classmethod def _noise_distance_keeps_clear_of_one( diff --git a/src/policyengine/disclosure/support.py b/src/policyengine/disclosure/support.py index 7d6b50c5..10c9a12f 100644 --- a/src/policyengine/disclosure/support.py +++ b/src/policyengine/disclosure/support.py @@ -5,6 +5,7 @@ count once: weights and contributions are summed by source first. """ +import math from typing import NamedTuple, Optional import numpy as np @@ -29,6 +30,7 @@ class Support(NamedTuple): n_sources: int n_eff: float top_share: float + rest_share: float = math.inf NO_SUPPORT = Support(0, 0.0, 0.0) @@ -45,14 +47,33 @@ def support( if not contributing.any(): return NO_SUPPORT _, by_source = np.unique(source_index[contributing], return_inverse=True) + weighted = weights[contributing] * contributions[contributing] source_weight = np.bincount(by_source, weights=weights[contributing]) - magnitude = np.bincount( - by_source, - weights=np.abs(weights[contributing] * contributions[contributing]), - ) + magnitude = np.bincount(by_source, weights=np.abs(weighted)) n_eff = float(source_weight.sum() ** 2 / (source_weight**2).sum()) top_share = float(magnitude.max() / magnitude.sum()) - return Support(len(source_weight), n_eff, top_share) + signed = np.bincount(by_source, weights=weighted) + return Support(len(source_weight), n_eff, top_share, _rest_share(signed)) + + +def _rest_share(signed: np.ndarray) -> float: + """The rest of a total next to its largest household: the p% measure. + + ``top_share`` works on magnitudes, so it misses a total whose other parts + cancel: the total is then one household's contribution. This measure + works on the signed total. It takes what is left without the largest + household, and what is left without the two largest, as a share of the + largest household's part. A small value means the number gives that part + away, to any reader or to the second-largest household. + """ + order = np.argsort(-np.abs(signed)) + largest = abs(float(signed[order[0]])) + if largest == 0: + return math.inf + total = float(signed.sum()) + without_one = abs(total - float(signed[order[0]])) + without_two = abs(total - float(signed[order[:2]].sum())) + return min(without_one, without_two) / largest def family_supports( @@ -83,7 +104,11 @@ def primary_tier( return WITHHELD, PRIVATE_BELOW_FLOOR if not one_area_per_record and cell.n_eff < settings.min_n_eff: return WITHHELD, PRIVATE_LOW_N_EFF - if not noise and cell.top_share > settings.max_top_share: + dominated = ( + cell.top_share > settings.max_top_share + or cell.rest_share < settings.min_rest_share + ) + if not noise and dominated: return WITHHELD, PRIVATE_DOMINANT if cell.n_eff < settings.approximate_below_n_eff: return APPROXIMATE, REASON_THIN_SUPPORT diff --git a/src/policyengine/disclosure/validate.py b/src/policyengine/disclosure/validate.py index e6294930..ef6cdd19 100644 --- a/src/policyengine/disclosure/validate.py +++ b/src/policyengine/disclosure/validate.py @@ -21,6 +21,9 @@ DEFAULT_MIN_ROUNDING_BASE, DIRECTIONS, ESTIMATOR_DIRECT, + GAIN, + LOSS, + NO_CHANGE, OPTIONAL_LABEL_COLUMNS, OUTCOMES, PUBLIC_COLUMNS, @@ -32,6 +35,7 @@ SHARE_OUTCOMES, SUPPORT_BANDS, TIERS, + TOO_CLOSE_TO_CALL, WITHHELD, PublicTable, ) @@ -121,6 +125,17 @@ def _check_frame( continue if not lo <= value <= hi: problems.append(f"{where}: the value is outside its interval") + if outcome in CHANGE_OUTCOMES and direction in DIRECTIONS: + if lo > 0: + expected = GAIN + elif hi < 0: + expected = LOSS + elif value == 0 and lo == 0 and hi == 0: + expected = NO_CHANGE + else: + expected = TOO_CLOSE_TO_CALL + if direction != expected: + problems.append(f"{where}: the direction does not match the interval") if row["support_band"] == BAND_UNDER_10 and value != 0: problems.append(f"{where}: a non-zero value has a support band under 10") if outcome in SHARE_OUTCOMES: diff --git a/src/policyengine/outputs/uk_geography_disclosure.py b/src/policyengine/outputs/uk_geography_disclosure.py index 48be972c..511eacfe 100644 --- a/src/policyengine/outputs/uk_geography_disclosure.py +++ b/src/policyengine/outputs/uk_geography_disclosure.py @@ -170,19 +170,26 @@ def areas_for(level: str) -> Optional[AreaWeights]: publish=level in self.levels, ) + regions = estimated.get(UK_REGION_LEVEL) for level, areas in estimated.items(): if level == UK_NATIONAL_LEVEL: continue - under_region = ( - level != UK_REGION_LEVEL - and UK_REGION_LEVEL in estimated - and estimated[UK_REGION_LEVEL].one_area_per_record - and (areas.one_area_per_record or level in parents) - ) - if under_region: - release.add_nesting(level, UK_REGION_LEVEL, parents.get(level)) - else: + if ( + level == UK_REGION_LEVEL + or regions is None + or not regions.one_area_per_record + ): release.add_nesting(level, UK_NATIONAL_LEVEL) + continue + if not areas.one_area_per_record and level not in parents: + # Against the nation alone, every withheld area goes into one + # remainder, and that nearly always passes. + raise ValueError( + f"Level {level!r} is a weight matrix. Pass " + f"parent_of={{{level!r}: ...}} with the region code of each " + "area, so its areas are checked against their region." + ) + release.add_nesting(level, UK_REGION_LEVEL, parents.get(level)) self._release = release self.refresh() @@ -236,7 +243,8 @@ def compute_uk_geography_disclosure( area_weights: Weights for a level, in place of the level's area code column, for example an area x household weight matrix. parent_of: For a level given as a weight matrix, the region code of - each of its areas. + each of its areas. Required for such a level when the household + output carries regions. acknowledge: Findings the caller accepts after review. prior_public: Earlier publications that count as already public. audit_dir: Private directory, outside any repository, for the exact diff --git a/src/policyengine/outputs/uk_geography_impact.py b/src/policyengine/outputs/uk_geography_impact.py index 2f48f252..6a283c63 100644 --- a/src/policyengine/outputs/uk_geography_impact.py +++ b/src/policyengine/outputs/uk_geography_impact.py @@ -1,5 +1,6 @@ """Longwise UK geography impact helpers.""" +import logging import warnings from pathlib import Path from typing import Optional @@ -12,16 +13,22 @@ default_local_search_dirs, ) +logger = logging.getLogger(__name__) + def warn_unprotected_uk_geography_output(name: str, level: str) -> None: - """Warn that an output without disclosure control is deprecated.""" - warnings.warn( + """Warn that an output without disclosure control is deprecated. + + policyengine-core silences Python warnings in every process that loads a + country model, so the message also goes to the log, which it leaves alone. + """ + message = ( f"{name} is deprecated. Its output is not disclosure-controlled, so it " "must not be published when the data is licensed survey microdata. Use " - f"compute_uk_geography_disclosure(levels=({level!r},)) instead.", - DeprecationWarning, - stacklevel=3, + f"compute_uk_geography_disclosure(levels=({level!r},)) instead." ) + logger.warning(message) + warnings.warn(message, DeprecationWarning, stacklevel=3) def _normalise_code(value) -> str: diff --git a/tests/test_constituency_impact.py b/tests/test_constituency_impact.py index c0dcdc92..d7b9b29c 100644 --- a/tests/test_constituency_impact.py +++ b/tests/test_constituency_impact.py @@ -257,7 +257,7 @@ def test_compute_constituency_impacts_requires_longwise_geography_column(): ) -def test_constituency_impacts_warn_that_they_are_not_disclosure_controlled(): +def test_constituency_impacts_warn_that_they_are_not_disclosure_controlled(caplog): household = { "household_net_income": [100.0], "household_weight": [1.0], @@ -273,3 +273,7 @@ def test_constituency_impacts_warn_that_they_are_not_disclosure_controlled(): ) with pytest.warns(DeprecationWarning, match=message): impact.run() + # policyengine-core silences Python warnings once a country model loads, + # so the message must also reach the log. + logged = [record.getMessage() for record in caplog.records] + assert sum("not disclosure-controlled" in line for line in logged) == 2 diff --git a/tests/test_disclosure_areas.py b/tests/test_disclosure_areas.py index 883b2d20..c4a06d68 100644 --- a/tests/test_disclosure_areas.py +++ b/tests/test_disclosure_areas.py @@ -24,7 +24,7 @@ def test_labels_and_a_one_nonzero_matrix_are_the_same_data(): from_labels = AreaWeights.from_labels(codes, weights) from_matrix = AreaWeights.from_matrix(matrix, order) - assert from_matrix.one_area_per_record + assert from_labels.one_area_per_record assert from_labels.area_codes == from_matrix.area_codes == ("A", "B", "C") assert np.array_equal(from_labels.labels, from_matrix.labels) assert np.allclose(from_labels.totals(values), from_matrix.totals(values)) @@ -34,6 +34,15 @@ def test_labels_and_a_one_nonzero_matrix_are_the_same_data(): ) +def test_weights_given_as_a_matrix_keep_the_matrix_rules(): + """One row of a dense matrix looks like labels, but is not label data.""" + one_row = AreaWeights.from_matrix(np.array([[5.0, 0.0, 2.0, 1.0]]), ["A"]) + + assert not one_row.one_area_per_record + assert one_row.labels.tolist() == ["A", UNASSIGNED, "A", "A"] + assert one_row.totals(np.ones(4)).tolist() == [8.0] + + def test_totals_equal_a_direct_weighted_sum(): codes, weights, _, _ = _labels_and_matrix() values = np.arange(1.0, 8.0) diff --git a/tests/test_disclosure_estimate.py b/tests/test_disclosure_estimate.py index 6c6a21ce..270e2cde 100644 --- a/tests/test_disclosure_estimate.py +++ b/tests/test_disclosure_estimate.py @@ -1,6 +1,8 @@ """Private estimates: values, support, intervals and rounding. Synthetic data.""" +import copy import math +import pickle import numpy as np import pandas as pd @@ -72,7 +74,7 @@ def _estimate(data, *, settings=None, noise_key=None, areas=None): def _cell(estimates, area, outcome): - table = estimates.table + table = estimates._table rows = table[(table["area_code"] == area) & (table["outcome"] == outcome)] assert len(rows) == 1 return rows.iloc[0] @@ -159,7 +161,10 @@ def test_a_cell_with_no_contributing_household_publishes_zero(): assert (cell["value"], cell["lo"], cell["hi"]) == (0.0, 0.0, 0.0) assert (cell["tier"], cell["direction"]) == (REPORTED, "no_change") - assert cell["support_band"] == BAND_UNDER_10 + # A zero rests on the households that did not change, as a zero share + # rests on its other side. + assert cell["n_sources"] == 0 + assert cell["support_band"] == BAND_30_PLUS assert (losing["value"], losing["tier"]) == (0.0, REPORTED) assert losing["support_band"] == BAND_30_PLUS @@ -218,6 +223,29 @@ def test_one_heavy_record_in_many_areas_fails_the_effective_sample_gate(): assert cell["support_band"] == BAND_UNDER_10 +def test_one_row_of_a_weight_matrix_keeps_the_effective_sample_gate(): + """A per-seat slice of a dense matrix has one weight per record.""" + n = 50 + weights = np.array([10_000.0] * 3 + [1.0] * 47) + data = _population({"a": (n, n)}) + data["weights"] = weights + + as_matrix = _estimate(data, areas=AreaWeights.from_matrix(weights[None, :], ["a"])) + as_labels = _estimate(data) + + matrix_cell = _cell(as_matrix, "a", "average_change") + label_cell = _cell(as_labels, "a", "average_change") + assert matrix_cell["n_eff"] == pytest.approx(label_cell["n_eff"]) + assert matrix_cell["n_eff"] < 10 + assert (matrix_cell["tier"], matrix_cell["private_reason"]) == ( + WITHHELD, + PRIVATE_LOW_N_EFF, + ) + assert matrix_cell["support_band"] == BAND_UNDER_10 + assert label_cell["tier"] == APPROXIMATE + assert label_cell["support_band"] == BAND_30_PLUS + + def test_an_even_weight_matrix_passes_and_bands_by_effective_sample(): n = 300 data = _population({"X": (n, n)}) @@ -287,7 +315,7 @@ def test_noise_reads_the_key_from_the_environment(monkeypatch): from_argument = _estimate(data, settings=settings, noise_key=TEST_KEY) assert from_environment.noise - pd.testing.assert_frame_equal(from_environment.table, from_argument.table) + pd.testing.assert_frame_equal(from_environment._table, from_argument._table) def test_shares_stay_between_zero_and_one_under_noise(): @@ -406,11 +434,11 @@ def test_intervals_are_reproducible_and_do_not_depend_on_the_chunk(monkeypatch): data["reformed"] != data["baseline"] ) - first = _estimate(data).table - again = _estimate(data).table + first = _estimate(data)._table + again = _estimate(data)._table monkeypatch.setattr(estimate_module, "_REPLICATE_CHUNK", 37) - chunked = _estimate(data).table - other_seed = _estimate(data, settings=DisclosureSettings(bootstrap_seed=1)).table + chunked = _estimate(data)._table + other_seed = _estimate(data, settings=DisclosureSettings(bootstrap_seed=1))._table pd.testing.assert_frame_equal(first, again) pd.testing.assert_frame_equal(first, chunked) @@ -425,7 +453,7 @@ def test_values_are_rounded_to_the_order_of_their_standard_error(): ) settings = DisclosureSettings() - table = _estimate(data).table + table = _estimate(data)._table for row in table[table["tier"] != WITHHELD].to_dict("records"): base = row["rounding_base"] @@ -434,6 +462,9 @@ def test_values_are_rounded_to_the_order_of_their_standard_error(): for column in ("value", "lo", "hi"): assert row[column] / base == pytest.approx(round(row[column] / base)) assert row["lo"] <= row["value"] <= row["hi"] + # The interval is rounded outward, never to the nearest step. + assert row["lo"] <= row["exact_lo"] + 1e-9 + assert row["hi"] >= row["exact_hi"] - 1e-9 if row["outcome"] != "share_no_change" and row["standard_error"] > base: assert base == 10.0 ** math.floor(math.log10(row["standard_error"])) average = _cell(_estimate(data), "A", "average_change") @@ -441,6 +472,27 @@ def test_values_are_rounded_to_the_order_of_their_standard_error(): assert average["direction"] == "gain" +def test_rounding_goes_down_up_or_to_the_nearest_step(): + assert estimate_module._round(27.0, 10.0) == 30.0 + assert estimate_module._round(27.0, 10.0, "down") == 20.0 + assert estimate_module._round(23.0, 10.0, "up") == 30.0 + assert estimate_module._round(-23.0, 10.0, "down") == -30.0 + assert estimate_module._round(0.3, 0.1, "down") == 0.3 + assert estimate_module._round(0.0437, 0.001, "up") == 0.044 + + +def test_an_interval_is_widened_to_include_its_estimate(): + """Skewed replicates can leave the estimate outside their percentiles.""" + draws = np.array([[2.0, 3.0, 4.0, 5.0], [2.0, 3.0, 4.0, 5.0]]) + + lo, hi, error = estimate_module._interval(draws, np.array([10.0, 0.0]), 0.05) + + assert hi[0] == 10.0 + assert lo[1] == 0.0 + assert lo[0] == pytest.approx(2.15) and hi[1] == pytest.approx(4.85) + assert error.tolist() == pytest.approx([np.std([2, 3, 4, 5], ddof=1)] * 2) + + def test_an_interval_across_zero_is_too_close_to_call(): rng = np.random.default_rng(6) data = _population({"A": (200, 0)}) @@ -468,6 +520,18 @@ def test_relative_change_is_withheld_where_it_is_undefined(): assert _cell(estimates, "A", "average_change")["tier"] != WITHHELD +def test_estimates_cannot_be_pickled_or_copied(): + """They hold exact values and, with noise on, the noise factors.""" + estimates = _estimate(_population({"A": (50, 30)})) + + with pytest.raises(TypeError, match="cannot be pickled"): + pickle.dumps(estimates) + with pytest.raises(TypeError, match="cannot be pickled"): + copy.deepcopy(estimates) + for name in ("table", "contributions", "source_index"): + assert not hasattr(estimates, name) + + def test_repr_holds_no_values(): estimates = _estimate(_population({"A": (50, 30)})) diff --git a/tests/test_disclosure_release.py b/tests/test_disclosure_release.py index 2b970dad..b1b719fe 100644 --- a/tests/test_disclosure_release.py +++ b/tests/test_disclosure_release.py @@ -220,7 +220,7 @@ def test_the_suggestion_skips_cells_that_add_nothing(): assert {finding.group_level for finding in findings} == {"region"} -def test_noise_never_protects_a_remainder(): +def test_noise_does_not_lift_the_count_floor_for_a_remainder(): settings = DisclosureSettings(noise=True) release = _release({**STRONG, "S1": 4}, settings=settings, key=TEST_KEY) @@ -279,7 +279,7 @@ def test_a_reference_level_is_checked_against_but_not_published(): assert set(release.to_public().to_frame()["level"]) == {"constituency"} -def _package_release(): +def _package_release(reference=()): data = _records(STRONG) seat_one = data["seats"] == "S1" broad = data["change"] @@ -287,7 +287,10 @@ def _package_release(): narrow = np.where(~seat_one & (np.arange(len(seat_one)) % 120 < 60), 500.0, narrow) release = Release() for reform, change in (("A", broad), ("B", narrow), ("P", broad + narrow)): - release.add(_estimates(data, "constituency", reform=reform, change=change)) + release.add( + _estimates(data, "constituency", reform=reform, change=change), + publish=reform not in reference, + ) release.add_package("P", ["A", "B"]) return release @@ -312,6 +315,46 @@ def test_a_package_gives_away_a_withheld_component(): assert release.check() == [] +def test_a_withheld_reference_component_is_hidden_too(): + """Public elsewhere means public where it passes; the rest is secret.""" + release = _package_release(reference=("B",)) + + findings = _blocking(release.check(), RECOVERABLE_REMAINDER) + + assert [(f.relation, f.group, f.members, f.suggested) for f in findings] == [ + ("package", "S1", ["B"], "A") + ] + + +def test_a_relation_between_references_only_is_not_checked(): + assert _package_release(reference=("A", "B", "P")).check() == [] + + +def test_withheld_cells_of_a_reference_child_level_are_hidden_too(): + data = _records({**STRONG, "S1": 4}) + release = Release() + release.add(_estimates(data, "region")) + release.add(_estimates(data, "constituency"), publish=False) + release.add_nesting("constituency", "region") + + findings = _blocking(release.check(), RECOVERABLE_REMAINDER) + + assert {(f.group, tuple(f.members), f.suggested) for f in findings} == { + ("R1", ("S1",), None) + } + # The seats are public elsewhere, so the fix is on the published side. + for finding in findings: + release.withhold( + level=finding.group_level, + area_code=finding.group, + reform="reform", + year=2026, + outcome="average_change" if finding.family == "change" else finding.family, + ) + assert release.check() == [] + assert set(release.to_public().to_frame()["level"]) == {"region"} + + def test_a_withheld_package_with_all_components_published_is_recoverable(): release = _package_release() for reform in ("A", "P"): @@ -401,6 +444,43 @@ def build(**options): Release(acknowledge=[RECOVERABLE_REMAINDER]) +def test_an_area_minus_the_public_areas_inside_it_is_checked(): + """A local authority that is one seat plus a few households of another.""" + data = _records(STRONG) + authorities = np.where(data["seats"] == "S1", "L1", "L2").astype(object) + authorities[np.flatnonzero(data["seats"] == "S2")[:5]] = "L1" + data["local_authoritys"] = authorities + release = Release(acknowledge=[UNCHECKED_OVERLAP]) + release.add(_estimates(data, "constituency")) + release.add(_estimates(data, "local_authority")) + + findings = release.check() + + blocking = _blocking(findings, RECOVERABLE_REMAINDER) + assert { + (f.relation, f.level, f.group, f.group_level, f.suggested) for f in blocking + } == {("containment", "constituency", "L1", "local_authority", "S1")} + # The rest is five households. Every family fails on it: also the share + # losing, because "none of these five lose" is a statement about five. + assert {f.family for f in blocking} == {"change", "share_gaining", "share_losing"} + overlap = [f for f in findings if f.code == UNCHECKED_OVERLAP] + assert [f.blocking for f in overlap] == [False] + assert "wholly inside" in overlap[0].message + with pytest.raises(DisclosureError, match=RECOVERABLE_REMAINDER): + release.to_public() + + for outcome in ("average_change", "share_gaining", "share_losing"): + release.withhold( + level="constituency", + area_code="S1", + reform="reform", + year=2026, + outcome=outcome, + ) + assert not any(finding.blocking for finding in release.check()) + release.to_public() + + def test_a_cell_that_is_already_public_is_not_a_secret(): cells = [ ("constituency", "S1", "reform", 2026, "average_change"), diff --git a/tests/test_disclosure_schema.py b/tests/test_disclosure_schema.py index 782dc3ff..417eb614 100644 --- a/tests/test_disclosure_schema.py +++ b/tests/test_disclosure_schema.py @@ -8,10 +8,13 @@ from pydantic import ValidationError from policyengine.disclosure import ( + AreaWeights, DisclosureSettings, PriorPublic, PublicRow, PublicTable, + Release, + estimate_change, ) from policyengine.disclosure.schema import ( BAND_10_TO_29, @@ -31,6 +34,8 @@ {"max_top_share": 0.6}, {"min_n_eff": 9.9}, {"change_tolerance": 0.0}, + {"change_tolerance": 1e-9}, + {"min_rest_share": 0.01}, {"n_replicates": 50}, ], ) @@ -39,6 +44,31 @@ def test_settings_cannot_be_loosened(loosened): DisclosureSettings(**loosened) +def test_settings_are_validated_again_before_use(): + """``model_copy`` and ``model_construct`` skip validation; the layer does not.""" + copied = DisclosureSettings().model_copy(update={"min_sources": 1}) + constructed = DisclosureSettings.model_construct(min_sources=1) + areas = AreaWeights.from_labels(["A"] * 20, [1.0] * 20) + + assert copied.min_sources == 1 + for settings in (copied, constructed): + with pytest.raises(ValidationError): + settings.revalidated() + with pytest.raises(ValidationError): + Release(settings) + with pytest.raises(ValidationError): + estimate_change( + level="region", + reform="r", + year=2026, + areas=areas, + source_ids=range(20), + baseline=[100.0] * 20, + reformed=[200.0] * 20, + settings=settings, + ) + + def test_settings_can_be_raised_for_special_conditions(): settings = DisclosureSettings(min_sources=30, min_n_eff=30.0) diff --git a/tests/test_disclosure_support.py b/tests/test_disclosure_support.py index 8bbb4e19..83a20321 100644 --- a/tests/test_disclosure_support.py +++ b/tests/test_disclosure_support.py @@ -113,6 +113,43 @@ def test_dominance_withholds_only_while_noise_is_off(): assert with_noise == (REPORTED, None) +def test_a_total_that_is_one_households_change_is_dominated(): + """Offsetting gains and losses leave the total equal to one household's.""" + contributions = np.array([1000.0] + [110.0] * 5 + [-110.0] * 5) + cell = support(np.ones(11), contributions, np.arange(11)) + + assert cell.top_share == pytest.approx(1000.0 / 2100.0) + assert cell.top_share < SETTINGS.max_top_share + assert cell.rest_share == pytest.approx(0.0, abs=1e-12) + assert primary_tier(cell, SETTINGS, **LABELS) == (WITHHELD, PRIVATE_DOMINANT) + with_noise = primary_tier(cell, SETTINGS, one_area_per_record=True, noise=True) + assert with_noise[0] != WITHHELD + + +def test_two_households_that_carry_a_total_are_dominated(): + contributions = np.array([490.0, 490.0] + [2.5] * 8) + cell = support(np.ones(10), contributions, np.arange(10)) + + assert cell.top_share == pytest.approx(0.49) + assert cell.rest_share == pytest.approx(20.0 / 490.0) + assert primary_tier(cell, SETTINGS, **LABELS) == (WITHHELD, PRIVATE_DOMINANT) + + +def test_the_rest_share_measures_what_is_left_without_the_largest(): + even = support(np.ones(40), np.full(40, 50.0), np.arange(40)) + skewed = support(np.ones(3), np.array([100.0, 50.0, 50.0]), np.arange(3)) + clones = support( + np.ones(4), np.array([60.0, 40.0, 50.0, 50.0]), np.array([0, 0, 1, 2]) + ) + cancelled = support(np.ones(2), np.array([70.0, -70.0]), np.array([0, 0])) + + assert even.rest_share == pytest.approx(38.0) + assert skewed.rest_share == pytest.approx(0.5) + assert clones.rest_share == pytest.approx(0.5) + assert cancelled.rest_share == float("inf") + assert primary_tier(even, SETTINGS, **LABELS)[0] == REPORTED + + def test_the_count_rule_comes_before_dominance(): cell = Support(3, 3.0, 0.9) diff --git a/tests/test_disclosure_validate.py b/tests/test_disclosure_validate.py index a3cac99c..fa3e5318 100644 --- a/tests/test_disclosure_validate.py +++ b/tests/test_disclosure_validate.py @@ -145,9 +145,22 @@ def test_a_non_zero_value_cannot_have_a_band_under_ten(frame): frame.loc[row, "support_band"] = "<10" assert any("under 10" in problem for problem in validate_public(frame)) + + +def test_a_zero_row_shows_the_households_it_was_measured_on(frame): zero = _first(frame, outcome="average_change", area_code="D") - assert frame.loc[zero, "support_band"] == "<10" + assert frame.loc[zero, "value"] == "0.0" + assert frame.loc[zero, "direction"] == "no_change" + assert frame.loc[zero, "support_band"] == "30+" + + +def test_a_direction_must_match_its_interval(frame): + row = _first(frame, tier="reported", outcome="average_change", area_code="C") + assert frame.loc[row, "direction"] == "gain" + frame.loc[row, "direction"] = "too_close_to_call" + + assert any("direction does not match" in p for p in validate_public(frame)) def test_shares_must_add_to_one_and_stay_in_range(frame): diff --git a/tests/test_local_authority_impact.py b/tests/test_local_authority_impact.py index aa32d154..3c536502 100644 --- a/tests/test_local_authority_impact.py +++ b/tests/test_local_authority_impact.py @@ -259,7 +259,7 @@ def test_compute_local_authority_impacts_requires_longwise_geography_column(): ) -def test_local_authority_impacts_warn_that_they_are_not_disclosure_controlled(): +def test_local_authority_impacts_warn_that_they_are_not_disclosure_controlled(caplog): household = { "household_net_income": [100.0], "household_weight": [1.0], @@ -277,3 +277,7 @@ def test_local_authority_impacts_warn_that_they_are_not_disclosure_controlled(): ) with pytest.warns(DeprecationWarning, match=message): impact.run() + # policyengine-core silences Python warnings once a country model loads, + # so the message must also reach the log. + logged = [record.getMessage() for record in caplog.records] + assert sum("not disclosure-controlled" in line for line in logged) == 2 diff --git a/tests/test_uk_geography_disclosure.py b/tests/test_uk_geography_disclosure.py index ff89d906..322163b9 100644 --- a/tests/test_uk_geography_disclosure.py +++ b/tests/test_uk_geography_disclosure.py @@ -237,13 +237,17 @@ def test_a_weight_matrix_can_replace_the_area_codes(): n = len(baseline["household_id"]) matrix = AreaWeights.from_matrix(np.full((2, n), 50.0), ["M1", "M2"]) + options = {"levels": ("constituency",), "area_weights": {"constituency": matrix}} + + with pytest.raises(ValueError, match="parent_of"): + _compute(baseline, reform, **options) + impact = _compute( baseline, reform, - levels=("constituency",), - area_weights={"constituency": matrix}, + parent_of={"constituency": {"M1": "NORTH_EAST", "M2": "LONDON"}}, + **options, ) - assert set(impact.public.to_frame()["area_code"]) == {"M1", "M2"} From 5ffa150a6ef05dc3c10e9dbee56cee3875779a81 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mar=C3=ADa=20Juaristi?= <127882282+juaristi22@users.noreply.github.com> Date: Wed, 14 Oct 2026 14:37:46 +0100 Subject: [PATCH 3/3] Address second review of the disclosure layer Apply the p% rule with every household as a target and any reader or any one other household as the attacker. Check each area against the areas of another level that meet it, and flag public areas that add up to a withheld area. Count a withheld parent as recoverable only when the remainder holds no records. Keep noise at least 8% from 1. Say in the changelog that the old functions log a warning. Co-Authored-By: Claude Opus 5.5 --- changelog.d/589.changed.md | 2 +- docs/disclosure.md | 33 +++-- src/policyengine/disclosure/release.py | 163 +++++++++++++++++-------- src/policyengine/disclosure/schema.py | 13 +- src/policyengine/disclosure/support.py | 38 ++++-- tests/test_disclosure_release.py | 110 +++++++++++++++-- tests/test_disclosure_schema.py | 2 + tests/test_disclosure_support.py | 48 ++++++++ 8 files changed, 317 insertions(+), 92 deletions(-) diff --git a/changelog.d/589.changed.md b/changelog.d/589.changed.md index c06e4f1e..2606e59e 100644 --- a/changelog.d/589.changed.md +++ b/changelog.d/589.changed.md @@ -1 +1 @@ -Deprecated `compute_uk_constituency_impacts`, `compute_uk_local_authority_impacts` and their `ConstituencyImpact` and `LocalAuthorityImpact` classes. Their output is not disclosure-controlled; use `compute_uk_geography_disclosure` instead. They now emit a `DeprecationWarning` and otherwise behave as before. +Deprecated `compute_uk_constituency_impacts`, `compute_uk_local_authority_impacts` and their `ConstituencyImpact` and `LocalAuthorityImpact` classes. Their output is not disclosure-controlled; use `compute_uk_geography_disclosure` instead. They now log a warning through the `policyengine` logger and emit a `DeprecationWarning`, and otherwise behave as before. diff --git a/docs/disclosure.md b/docs/disclosure.md index cb61fcd2..7cbffe84 100644 --- a/docs/disclosure.md +++ b/docs/disclosure.md @@ -57,13 +57,15 @@ must pass: the households that gain and those that do not. **Dominance.** While noise is off, a number is withheld in two cases: - one survey household carries more than half of it, by size; -- what is left of it without its largest household, or without its two - largest, is less than 10% of the largest household's part (the p% rule). - This catches a total whose other parts cancel, and a total that two - households carry between them. +- the p% rule: a reader of the number, or one household that knows its own + part, could estimate some household's part to within 10%. Every household + is a possible target and every other household a possible reader, so the + rule also catches a total whose other parts cancel. -These rules protect one or two households. They do not stop a number that -three or more households carry almost alone, as long as 10 contribute. +These rules protect each household against any reader and against any one +other household. They do not stop two or more households that pool what they +know: a number that three households carry almost alone passes, as long as 10 +contribute. **Effective sample.** Weights that differ a lot make an estimate less precise than its household count suggests. Where the data gives one area and one @@ -105,8 +107,9 @@ looks for such remainders: - child areas against their parent area (constituencies in regions, regions in the nation); - the components of a package against the package; -- an area against the areas of another level that lie wholly inside it, when - two levels that do not nest are released together; +- when two levels that do not nest are released together, an area against + the areas of the other level that lie wholly inside it, and against the + areas that meet it (their sum minus the area is what they hold outside it); - earlier publications, passed as `prior_public`. A remainder gets the rules of a cell. Estimates added with `publish=False` @@ -143,9 +146,9 @@ Two findings can be accepted after review, with `acknowledge=[...]`: - `unchecked_overlap`: levels that overlap without nesting, such as constituencies and local authorities, are released together. The release - checks each area against the areas wholly inside it. It does not check - other differences between overlapping areas, and on weight-matrix data it - checks none. + checks each area against the areas wholly inside it and against the areas + that meet it. It does not check other differences between overlapping + areas, and on weight-matrix data it checks none. - `noise_reference_public`: an earlier table shows the same cells under another noise setting. @@ -155,7 +158,13 @@ Noise is off by default. With `DisclosureSettings(noise=True)`, each survey household's change is scaled by a fixed factor between 8% and 12% away from 1, derived from its id and a secret key. The same household has the same factor in every table and every rebuild. Noise replaces the dominance rule; the count -floor still applies. +floor still applies. The range cannot start closer to 1 than 8%. + +A number that one household dominates, published with noise, gives that +household's part to within 8 to 12%. The p% rule treats 10% as the line, so +noise puts such a number right at that line, not safely past it. Whether the +UK Data Service accepts noise in place of the dominance rule is a licence +question, which is why noise stays off by default. - **Size.** The noise in a total is about a tenth of that total's sampling error. diff --git a/src/policyengine/disclosure/release.py b/src/policyengine/disclosure/release.py index 464db136..6e45895c 100644 --- a/src/policyengine/disclosure/release.py +++ b/src/policyengine/disclosure/release.py @@ -281,7 +281,11 @@ def _known_cells(self) -> set[tuple[str, str, str, int, str]]: def _passes( self, estimates: AreaEstimates, weights: np.ndarray, family: str ) -> tuple[bool, bool]: - """Whether a remainder may be recovered, and whether it is empty. + """Whether a remainder passes, and whether it holds no records at all. + + Only a remainder with no records gives a withheld parent away by + addition. A remainder whose records happen not to change still hides + them: a reader cannot know that they did not. A remainder gets the rules of a cell, no more and no less. The count floor always applies. With noise on, what subtraction recovers is the @@ -301,8 +305,7 @@ def _passes( one_area_per_record=estimates.areas.one_area_per_record, noise=estimates.noise, ) - empty = supports[0].n_sources == 0 - return tier != WITHHELD, empty + return tier != WITHHELD, not len(records) @staticmethod def _require_same_records(first: AreaEstimates, second: AreaEstimates) -> None: @@ -638,8 +641,9 @@ def nests(child: str, parent: str) -> bool: for outer, inner in ((first, second), (second, first)) ] scope = ( - "Areas that lie wholly inside an area of the other level are " - "checked. Other differences between overlapping areas are not." + "Each area is checked against the areas of the other level that " + "lie wholly inside it and against those that meet it. Other " + "differences between overlapping areas are not." if all(contained) else "Differences between overlapping areas are not checked." ) @@ -664,11 +668,14 @@ def _check_containment( known: set[tuple[str, str, str, int, str]], findings: list[Finding], ) -> bool: - """Check an outer area against the inner areas wholly inside it. - - An outer area minus its public inner areas is a remainder. Returns - False when the two levels do not allow the check: they must hold the - same records with one area and one weight per record. + """Check each outer area against the inner areas inside it and around it. + + Two remainders come from one outer area: the area minus the public + inner areas wholly inside it, and the public inner areas that meet it + minus the area. The second is the households those inner areas hold + outside it. Returns False when the two levels do not allow the check: + they must hold the same records with one area and one weight per + record. """ outer, inner = self._estimates[outer_key], self._estimates[inner_key] outer_level, reform, year = outer_key @@ -688,59 +695,111 @@ def _check_containment( pairs = pairs[inner_weights > 0].drop_duplicates() whole = pairs[~pairs["inner"].duplicated(keep=False) & (pairs["outer"] != "")] inside = whole.groupby("outer")["inner"].apply(sorted).to_dict() + meeting = ( + pairs[pairs["outer"] != ""] + .groupby("outer")["inner"] + .apply(sorted) + .to_dict() + ) for family in FAMILIES: outer_tiers = tiers[(outer_key, family)] inner_tiers = tiers[(inner_key, family)] sources = inner.family_sources(family) - for code, members in inside.items(): - if ( - outer_tiers.get(code) == WITHHELD - and (outer_level, code, reform, year, family) not in known - ): - continue - public = [ - c - for c in members - if inner_tiers[c] != WITHHELD - or (inner_level, c, reform, year, family) in known - ] - if not public: - continue - weights = outer.areas.record_weights([code]) - weights = np.clip(weights - inner.areas.record_weights(public), 0, None) - passes, _ = self._passes(outer, weights, family) - if passes: - continue + + def public(cell: str) -> bool: + return ( + inner_tiers[cell] != WITHHELD + or (inner_level, cell, reform, year, family) in known + ) + + def suggest(cells: list[str]) -> Optional[str]: + # A cell that is public elsewhere cannot be withheld here. shown = [ c - for c in public + for c in cells if self._published[inner_key] and inner_tiers[c] != WITHHELD and sources[c] > 0 ] - findings.append( - Finding( - code=RECOVERABLE_REMAINDER, - blocking=True, - relation="containment", - message=( - f"{outer_level} {code!r} minus the public " - f"{inner_level} cells inside it fails the rules, so " - "subtraction recovers the rest. Withhold one of " - f"those cells, or the {outer_level} cell." - ), - level=inner_level, - reform=reform, - year=year, - family=family, - group=code, - group_level=outer_level, - members=sorted(set(members) - set(public)), - suggested=min( - shown, key=lambda c: (sources[c], c), default=None - ), - ) + return min(shown, key=lambda c: (sources[c], c), default=None) + + location: dict[str, Any] = { + "relation": "containment", + "level": inner_level, + "reform": reform, + "year": year, + "family": family, + "group_level": outer_level, + } + for code in outer.areas.area_codes: + outer_public = ( + outer_tiers.get(code) != WITHHELD + or (outer_level, code, reform, year, family) in known ) + members = inside.get(code, []) + shown_inside = [c for c in members if public(c)] + if shown_inside: + weights = outer.areas.record_weights([code]) + weights -= inner.areas.record_weights(shown_inside) + weights = np.clip(weights, 0, None) + weights[np.abs(weights) < 1e-9] = 0.0 + passes, empty = self._passes(outer, weights, family) + if not outer_public and empty: + findings.append( + Finding( + code=RECOVERABLE_PARENT, + blocking=True, + message=( + f"{outer_level} {code!r} is withheld but the " + f"public {inner_level} cells inside it add up " + "to it, so addition recovers it." + ), + group=code, + **location, + ) + ) + elif outer_public and not passes: + findings.append( + Finding( + code=RECOVERABLE_REMAINDER, + blocking=True, + message=( + f"{outer_level} {code!r} minus the public " + f"{inner_level} cells inside it fails the " + "rules, so subtraction recovers the rest. " + "Withhold one of those cells, or the " + f"{outer_level} cell." + ), + group=code, + members=sorted(set(members) - set(shown_inside)), + suggested=suggest(shown_inside), + **location, + ) + ) + around = meeting.get(code, []) + if not outer_public or not around or not all(map(public, around)): + continue + weights = inner.areas.record_weights(around) + weights -= outer.areas.record_weights([code]) + passes, _ = self._passes(outer, np.clip(weights, 0, None), family) + if not passes: + findings.append( + Finding( + code=RECOVERABLE_REMAINDER, + blocking=True, + message=( + f"The public {inner_level} cells that meet " + f"{outer_level} {code!r} cover it. Their sum " + "minus it fails the rules, so subtraction " + "recovers the households they hold outside it. " + "Withhold one of those cells, or the " + f"{outer_level} cell." + ), + group=code, + suggested=suggest(around), + **location, + ) + ) return True def _check_prior( diff --git a/src/policyengine/disclosure/schema.py b/src/policyengine/disclosure/schema.py index 380d9d2b..af8792cf 100644 --- a/src/policyengine/disclosure/schema.py +++ b/src/policyengine/disclosure/schema.py @@ -116,6 +116,12 @@ ACKNOWLEDGEABLE_FINDINGS = (UNCHECKED_OVERLAP, NOISE_REFERENCE_PUBLIC) +# The lowest distance of a noise factor from 1. Noise stands in for the +# dominance rule, so it may not come closer to 1 than the p% rule's 10% +# allows, give or take the noise range. +MIN_NOISE_DISTANCE = 0.08 + + def is_power_of_ten(value: float) -> bool: if not math.isfinite(value) or value <= 0: return False @@ -142,6 +148,7 @@ class DisclosureSettings(BaseModel): change_tolerance: float = Field(default=1.0, ge=1.0) noise: bool = False # Distance of a noise factor from 1: uniform on this range, either side. + # Noise replaces the dominance rule, so the range cannot start below 8%. noise_distance: tuple[float, float] = (0.08, 0.12) n_replicates: int = Field(default=1000, ge=200) interval: float = Field(default=0.90, gt=0.5, lt=1.0) @@ -164,8 +171,10 @@ def _noise_distance_keeps_clear_of_one( cls, value: tuple[float, float] ) -> tuple[float, float]: low, high = value - if not 0.0 < low < high < 1.0: - raise ValueError("noise_distance must satisfy 0 < low < high < 1") + if not MIN_NOISE_DISTANCE <= low < high < 1.0: + raise ValueError( + f"noise_distance must satisfy {MIN_NOISE_DISTANCE} <= low < high < 1" + ) return value @field_validator("min_rounding_base") diff --git a/src/policyengine/disclosure/support.py b/src/policyengine/disclosure/support.py index 10c9a12f..6b7e967e 100644 --- a/src/policyengine/disclosure/support.py +++ b/src/policyengine/disclosure/support.py @@ -57,23 +57,35 @@ def support( def _rest_share(signed: np.ndarray) -> float: - """The rest of a total next to its largest household: the p% measure. + """The p% measure on signed parts, for every household and every attacker. ``top_share`` works on magnitudes, so it misses a total whose other parts - cancel: the total is then one household's contribution. This measure - works on the signed total. It takes what is left without the largest - household, and what is left without the two largest, as a share of the - largest household's part. A small value means the number gives that part - away, to any reader or to the second-largest household. + cancel. This measure works on the signed parts, one per survey household. + For a household *i*, the rest of the total without it is ``R = T - c_i``. + Any reader of the total estimates ``c_i`` with error ``R``. A household + *j* that knows its own part estimates it with error ``R - c_j``. The + measure is the smallest such error as a share of ``|c_i|``, over every + *i* and every *j*. A small value means the number gives one household's + part away, to any reader or to one other household. """ - order = np.argsort(-np.abs(signed)) - largest = abs(float(signed[order[0]])) - if largest == 0: + parts = signed[signed != 0] + if not len(parts): return math.inf - total = float(signed.sum()) - without_one = abs(total - float(signed[order[0]])) - without_two = abs(total - float(signed[order[:2]].sum())) - return min(without_one, without_two) / largest + rest = float(signed.sum()) - parts + order = np.argsort(parts, kind="stable") + ranked = parts[order] + rank = np.empty(len(parts), dtype=int) + rank[order] = np.arange(len(parts)) + # The part nearest to each R is next to where R would sort. If that part + # is the household's own, the nearest other part is one step further. + slot = np.searchsorted(ranked, rest) + error = np.abs(rest) + for step in (-2, -1, 0, 1): + position = slot + step + valid = (position >= 0) & (position < len(ranked)) & (position != rank) + nearest = ranked[np.clip(position, 0, len(ranked) - 1)] + error = np.where(valid, np.minimum(error, np.abs(rest - nearest)), error) + return float(np.min(error / np.abs(parts))) def family_supports( diff --git a/tests/test_disclosure_release.py b/tests/test_disclosure_release.py index b1b719fe..d90cdbd5 100644 --- a/tests/test_disclosure_release.py +++ b/tests/test_disclosure_release.py @@ -444,6 +444,28 @@ def build(**options): Release(acknowledge=[RECOVERABLE_REMAINDER]) +def _withhold_suggestions(release, rounds=20): + """Resolve blocking findings one at a time, as a person would.""" + withheld = set() + for _ in range(rounds): + blocking = [finding for finding in release.check() if finding.blocking] + if not blocking: + return withheld + finding = blocking[0] + level = finding.level if finding.suggested else finding.group_level + area = finding.suggested or finding.group + withheld.add((level, area)) + family = finding.family + release.withhold( + level=level, + area_code=area, + reform="reform", + year=2026, + outcome="average_change" if family == "change" else family, + ) + raise AssertionError("findings did not resolve") + + def test_an_area_minus_the_public_areas_inside_it_is_checked(): """A local authority that is one seat plus a few households of another.""" data = _records(STRONG) @@ -457,11 +479,16 @@ def test_an_area_minus_the_public_areas_inside_it_is_checked(): findings = release.check() blocking = _blocking(findings, RECOVERABLE_REMAINDER) + # The five households are recovered three ways: L1 minus S1, S1's + # authorities minus S1, and S2 to S6 minus L2. assert { (f.relation, f.level, f.group, f.group_level, f.suggested) for f in blocking - } == {("containment", "constituency", "L1", "local_authority", "S1")} - # The rest is five households. Every family fails on it: also the share - # losing, because "none of these five lose" is a statement about five. + } == { + ("containment", "constituency", "L1", "local_authority", "S1"), + ("containment", "local_authority", "S1", "constituency", "L1"), + ("containment", "constituency", "L2", "local_authority", "S2"), + } + # Also the share losing: "none of these five lose" is a statement about five. assert {f.family for f in blocking} == {"change", "share_gaining", "share_losing"} overlap = [f for f in findings if f.code == UNCHECKED_OVERLAP] assert [f.blocking for f in overlap] == [False] @@ -469,18 +496,77 @@ def test_an_area_minus_the_public_areas_inside_it_is_checked(): with pytest.raises(DisclosureError, match=RECOVERABLE_REMAINDER): release.to_public() - for outcome in ("average_change", "share_gaining", "share_losing"): - release.withhold( - level="constituency", - area_code="S1", - reform="reform", - year=2026, - outcome=outcome, - ) - assert not any(finding.blocking for finding in release.check()) + withheld = _withhold_suggestions(release) + + # Two cells must go: one that breaks L1 minus S1, one that breaks the + # seats around L2. + assert len(withheld) == 2 + assert ("constituency", "S2") in withheld release.to_public() +def test_an_area_covered_by_the_areas_that_meet_it_is_checked(): + """Two seats that cover an authority, each with a few households outside it.""" + seats = ["S1"] * 200 + ["S2"] * 200 + ["S3"] * 200 + ["S4"] * 200 + authorities = ( + ["L1"] * 197 + + ["L2"] * 3 + + ["L1"] * 198 + + ["L2"] * 2 + + ["L2"] * 100 + + ["L3"] * 100 + + ["L2"] * 100 + + ["L3"] * 100 + ) + n = len(seats) + baseline = np.random.default_rng(1).uniform(20_000.0, 60_000.0, n) + data = { + "seats": np.array(seats, dtype=object), + "local_authoritys": np.array(authorities, dtype=object), + "weights": np.full(n, 100.0), + "source_ids": np.arange(n), + "baseline": baseline, + "change": np.full(n, 200.0), + } + release = Release(acknowledge=[UNCHECKED_OVERLAP]) + release.add(_estimates(data, "constituency")) + release.add(_estimates(data, "local_authority")) + + blocking = _blocking(release.check(), RECOVERABLE_REMAINDER) + + # No area of either level lies wholly inside an area of the other, yet + # S1 + S2 - L1 gives the five households that the two seats hold in L2. + assert {(f.level, f.group, f.group_level) for f in blocking} == { + ("constituency", "L1", "local_authority") + } + assert {f.suggested for f in blocking} <= {"S1", "S2"} + _withhold_suggestions(release) + release.to_public() + + +def test_public_areas_inside_a_withheld_area_cannot_add_up_to_it(): + data = _records(STRONG) + data["local_authoritys"] = np.where( + np.isin(data["seats"], ["S1", "S2"]), "L1", "L2" + ).astype(object) + release = Release(acknowledge=[UNCHECKED_OVERLAP]) + release.add(_estimates(data, "constituency")) + release.add(_estimates(data, "local_authority")) + release.withhold( + level="local_authority", + area_code="L1", + reform="reform", + year=2026, + outcome="average_change", + ) + + findings = _blocking(release.check(), RECOVERABLE_PARENT) + + assert [(f.relation, f.group, f.group_level, f.family) for f in findings] == [ + ("containment", "L1", "local_authority", "change") + ] + + def test_a_cell_that_is_already_public_is_not_a_secret(): cells = [ ("constituency", "S1", "reform", 2026, "average_change"), diff --git a/tests/test_disclosure_schema.py b/tests/test_disclosure_schema.py index 417eb614..5db9552d 100644 --- a/tests/test_disclosure_schema.py +++ b/tests/test_disclosure_schema.py @@ -36,6 +36,8 @@ {"change_tolerance": 0.0}, {"change_tolerance": 1e-9}, {"min_rest_share": 0.01}, + {"noise_distance": (1e-6, 2e-6)}, + {"noise_distance": (0.05, 0.1)}, {"n_replicates": 50}, ], ) diff --git a/tests/test_disclosure_support.py b/tests/test_disclosure_support.py index 83a20321..4ad32809 100644 --- a/tests/test_disclosure_support.py +++ b/tests/test_disclosure_support.py @@ -135,6 +135,54 @@ def test_two_households_that_carry_a_total_are_dominated(): assert primary_tier(cell, SETTINGS, **LABELS) == (WITHHELD, PRIVATE_DOMINANT) +def test_any_household_can_be_the_attacker(): + """The -400 household knows its own part and learns the +1,000.""" + contributions = np.array([1000.0, 500.0, -400.0] + [-495.0 / 17] * 17) + cell = support(np.ones(20), contributions, np.arange(20)) + + assert cell.top_share < SETTINGS.max_top_share + assert cell.rest_share == pytest.approx(5.0 / 1000.0) + assert primary_tier(cell, SETTINGS, **LABELS) == (WITHHELD, PRIVATE_DOMINANT) + + +def test_any_household_can_be_the_target(): + """Two large parts cancel, and the total gives the +600 away.""" + contributions = np.array([1000.0, -990.0, 600.0] + [2.0, -2.0] * 8 + [2.0]) + cell = support(np.ones(20), contributions, np.arange(20)) + + # Any reader is 12 off; a +2 household that knows its own part is 10 off. + assert cell.rest_share == pytest.approx(10.0 / 600.0) + assert primary_tier(cell, SETTINGS, **LABELS) == (WITHHELD, PRIVATE_DOMINANT) + + +def test_the_rest_share_agrees_with_a_direct_search(): + """Every household as target, every other household or none as attacker.""" + rng = np.random.default_rng(0) + for _ in range(500): + n = int(rng.integers(1, 25)) + parts = rng.normal(0.0, 100.0, n) * rng.choice([1.0, 0.01], n) + cell = support(np.ones(n), parts, np.arange(n)) + total = parts.sum() + direct = min( + min( + [abs(total - part)] + + [abs(total - part - other) for j, other in enumerate(parts) if j != i] + ) + / abs(part) + for i, part in enumerate(parts) + ) + assert cell.rest_share == pytest.approx(direct, rel=1e-9, abs=1e-12) + + +def test_three_households_that_carry_a_total_are_not_caught(): + """The rules protect against one reader or one household, not a coalition.""" + contributions = np.array([500.0] * 3 + [2.0] * 17) + cell = support(np.ones(20), contributions, np.arange(20)) + + assert cell.rest_share > SETTINGS.min_rest_share + assert primary_tier(cell, SETTINGS, **LABELS)[0] != WITHHELD + + def test_the_rest_share_measures_what_is_left_without_the_largest(): even = support(np.ones(40), np.full(40, 50.0), np.arange(40)) skewed = support(np.ones(3), np.array([100.0, 50.0, 50.0]), np.arange(3))