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..2606e59e --- /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 log a warning through the `policyengine` logger and 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..7cbffe84 --- /dev/null +++ b/docs/disclosure.md @@ -0,0 +1,239 @@ +--- +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, 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 in two cases: + +- one survey household carries more than half of it, by size; +- 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 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 +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 +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; +- 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` +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 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 +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 released together. The release + 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. + +## 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. 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. +- **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. 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 + +`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, 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. + +## 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: + +```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..e7058d9c --- /dev/null +++ b/src/policyengine/disclosure/areas.py @@ -0,0 +1,216 @@ +"""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. 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 +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, + 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. + 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, + 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) + 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, + 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. + + 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.") + 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: 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, from_matrix=True) + 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: + """Whether the data was given as one area and one weight per record.""" + return not self._from_matrix + + @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.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 new file mode 100644 index 00000000..958e30fa --- /dev/null +++ b/src/policyengine/disclosure/estimate.py @@ -0,0 +1,472 @@ +"""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. + + 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, + *, + level: str, + reform: str, + year: int, + areas: AreaWeights, + source_index: np.ndarray, + settings: DisclosureSettings, + noise: bool, + contributions: dict[str, np.ndarray], + table: pd.DataFrame, + ) -> None: + self.level = level + self.reform = reform + self.year = year + self.areas = areas + self.settings = settings + self.noise = noise + 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. + return ( + f"AreaEstimates(level={self.level!r}, reform={self.reform!r}, " + 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] + 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 _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 + 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()).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 + 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(): + bounds[outcome] = _interval(draws, exact[outcome], alpha) + + 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]), + # 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, + code, + reform, + year, + SHARE_NO_CHANGE, + SHARE_NO_CHANGE, + position, + exact, + bounds, + rounded[SHARE_NO_CHANGE], + third, + third_reason, + third_cell, + third_cell, + one_area, + ) + ) + + return AreaEstimates( + level=level, + reform=reform, + year=int(year), + areas=areas, + source_index=source_index, + settings=settings, + noise=settings.noise, + contributions=contributions, + 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, + band_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, + "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/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..6e45895c --- /dev/null +++ b/src/policyengine/disclosure/release.py @@ -0,0 +1,858 @@ +"""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()).revalidated() + 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: 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( + "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(tiers, known, 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 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 + 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], + family, + ) + tier, _ = cell_tier( + supports, + self.settings, + one_area_per_record=estimates.areas.one_area_per_record, + noise=estimates.noise, + ) + return tier != WITHHELD, not len(records) + + @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: + 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}." + ) + # 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: + 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: + 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. + # A cell that is public elsewhere cannot be + # withheld here either. + shown = [ + c + for c in children + if self._published[key] + and 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: + 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)) + # 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 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 self._published[part_key] + and 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, + 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) + + 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 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 = ( + "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." + ) + findings.append( + Finding( + code=UNCHECKED_OVERLAP, + blocking=UNCHECKED_OVERLAP not in self._acknowledged, + relation="release", + message=( + 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 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 + 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() + 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) + + 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 cells + if self._published[inner_key] + and inner_tiers[c] != WITHHELD + and sources[c] > 0 + ] + 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( + 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..af8792cf --- /dev/null +++ b/src/policyengine/disclosure/schema.py @@ -0,0 +1,278 @@ +"""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", + "rest_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) + + +# 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 + 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) + # 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, 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 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) + bootstrap_seed: int = 0 + min_rounding_base: dict[str, float] = Field( + 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( + cls, value: tuple[float, float] + ) -> tuple[float, float]: + low, high = value + 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") + @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..6b7e967e --- /dev/null +++ b/src/policyengine/disclosure/support.py @@ -0,0 +1,166 @@ +"""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. +""" + +import math +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 + rest_share: float = math.inf + + +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) + weighted = weights[contributing] * contributions[contributing] + source_weight = np.bincount(by_source, weights=weights[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()) + 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 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. 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. + """ + parts = signed[signed != 0] + if not len(parts): + return math.inf + 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( + 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 + 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 + 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..ef6cdd19 --- /dev/null +++ b/src/policyengine/disclosure/validate.py @@ -0,0 +1,211 @@ +"""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, + GAIN, + LOSS, + NO_CHANGE, + OPTIONAL_LABEL_COLUMNS, + OUTCOMES, + PUBLIC_COLUMNS, + REASON_DISCLOSURE_CONTROL, + REASON_THIN_SUPPORT, + SHARE_GAINING, + SHARE_LOSING, + SHARE_NO_CHANGE, + SHARE_OUTCOMES, + SUPPORT_BANDS, + TIERS, + TOO_CLOSE_TO_CALL, + 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 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: + 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..511eacfe --- /dev/null +++ b/src/policyengine/outputs/uk_geography_disclosure.py @@ -0,0 +1,275 @@ +"""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, + ) + + regions = estimated.get(UK_REGION_LEVEL) + for level, areas in estimated.items(): + if level == UK_NATIONAL_LEVEL: + continue + 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() + 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. 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 + 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..6a283c63 100644 --- a/src/policyengine/outputs/uk_geography_impact.py +++ b/src/policyengine/outputs/uk_geography_impact.py @@ -1,5 +1,7 @@ """Longwise UK geography impact helpers.""" +import logging +import warnings from pathlib import Path from typing import Optional @@ -11,6 +13,23 @@ 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. + + 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." + ) + logger.warning(message) + warnings.warn(message, DeprecationWarning, stacklevel=3) + def _normalise_code(value) -> str: if value is None or pd.isna(value): 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..d7b9b29c 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,25 @@ 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(caplog): + 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() + # 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 new file mode 100644 index 00000000..c4a06d68 --- /dev/null +++ b/tests/test_disclosure_areas.py @@ -0,0 +1,132 @@ +"""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_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)) + assert np.allclose( + from_labels.totals(values, multipliers), + from_matrix.totals(values, multipliers), + ) + + +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) + 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..270e2cde --- /dev/null +++ b/tests/test_disclosure_estimate.py @@ -0,0 +1,558 @@ +"""Private estimates: values, support, intervals and rounding. Synthetic data.""" + +import copy +import math +import pickle + +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") + # 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 + + +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_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)}) + 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"] + # 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") + assert average["exact_lo"] < average["exact_value"] < average["exact_hi"] + 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)}) + 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_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)})) + + 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..d90cdbd5 --- /dev/null +++ b/tests/test_disclosure_release.py @@ -0,0 +1,678 @@ +"""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_does_not_lift_the_count_floor_for_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(reference=()): + 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), + publish=reform not in reference, + ) + 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_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"): + 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 _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) + 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) + # 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"), + ("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] + assert "wholly inside" in overlap[0].message + with pytest.raises(DisclosureError, match=RECOVERABLE_REMAINDER): + release.to_public() + + 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"), + ("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..5db9552d --- /dev/null +++ b/tests/test_disclosure_schema.py @@ -0,0 +1,186 @@ +"""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 ( + AreaWeights, + DisclosureSettings, + PriorPublic, + PublicRow, + PublicTable, + Release, + estimate_change, +) +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}, + {"change_tolerance": 1e-9}, + {"min_rest_share": 0.01}, + {"noise_distance": (1e-6, 2e-6)}, + {"noise_distance": (0.05, 0.1)}, + {"n_replicates": 50}, + ], +) +def test_settings_cannot_be_loosened(loosened): + with pytest.raises(ValidationError): + 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) + + 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..4ad32809 --- /dev/null +++ b/tests/test_disclosure_support.py @@ -0,0 +1,242 @@ +"""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_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_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)) + 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) + + 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..fa3e5318 --- /dev/null +++ b/tests/test_disclosure_validate.py @@ -0,0 +1,217 @@ +"""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)) + + +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, "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): + 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..3c536502 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,27 @@ 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(caplog): + 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() + # 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 new file mode 100644 index 00000000..322163b9 --- /dev/null +++ b/tests/test_uk_geography_disclosure.py @@ -0,0 +1,277 @@ +"""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"]) + + options = {"levels": ("constituency",), "area_weights": {"constituency": matrix}} + + with pytest.raises(ValueError, match="parent_of"): + _compute(baseline, reform, **options) + + impact = _compute( + baseline, + reform, + parent_of={"constituency": {"M1": "NORTH_EAST", "M2": "LONDON"}}, + **options, + ) + 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