Improve cross-machine reproducibility "fingerprint" test for area weights

[donboyd5]

@martinholmer raised a concern about whether the current "fingerprint" test for reproducibility across runs and machines is sufficient. It is not. This issue and proposed solution, prepared in collaboration with Claude, reflects that conversation with @martinholmer.

Goals

Area weight fingerprint tests should:

1. Catch legitimate differences — any real change to code, data, targets, or solver parameters that alters area weights should be detected, even if the change is small (e.g., 1% shift in a small area's weight sum)
2. No false positives from numerical noise — different hardware (x86 vs ARM, different SIMD), different BLAS implementations (MKL vs OpenBLAS vs Accelerate), or different operation ordering in algebraically equivalent code should not cause test failures
3. Scale to thousands of areas — if we add county weighting, we may have 3,000+ areas; the test must not become slow or produce spurious failures due to sheer number of comparisons
4. Be fast — the fingerprint check itself should take seconds, not minutes; I/O (reading weight files) will dominate, not computation

Current approach and its weaknesses

The current fingerprint (in _compute_fingerprint) works as follows:

• For each area, read the weight file, round each weight to the nearest integer, sum the rounded values → one integer per area
• Sort areas alphabetically, join as AK:12345|AL:67890|..., take first 16 chars of SHA-256

Problems:

• Insensitive to small areas: Area weights are often in the range 0.1–8.0. Rounding to integers destroys most of the information. A weight shifting from 1.5 to 2.4 (60% change) is invisible if both round to 2.
• Distribution-blind: Two completely different weight vectors with the same integer sum produce the same fingerprint. Redistribution of weight across records goes undetected.
• Rounding boundary false positives: A weight of 2.4999 on one machine and 2.5001 on another (a ~1e-5 relative difference — pure numerical noise) rounds to 2 vs 3, changing the sum and failing the test.
• Hash provides no diagnostics: When the SHA-256 check fails, you know something changed but not what or where. The per-area integer sum comparison helps, but inherits the insensitivity and boundary problems above.

Context: solver numerical properties

The area weight solver uses Clarabel (interior-point QP) with convergence tolerances of 1e-7 on gap and feasibility. Clarabel uses its own QDLDL factorization (not BLAS), so given identical inputs it aims for cross-machine determinism. However, the inputs to Clarabel (the constraint matrix from Tax-Calculator output) depend on numpy/BLAS, so independent make data runs on different machines can produce slightly different inputs, leading to solver path divergence.

Expected cross-machine differences in aggregate weight statistics: O(1e-6) to O(1e-4) relative. Expected differences from a real code/data change: typically > 1% relative. There is a clean separation between noise and signal — the test tolerance just needs to sit between them.

[donboyd5]

Alternatives considered and recommendation

Alternatives

A. Full-vector hash at double precision
SHA-256 of the entire weight vector at full floating-point precision. Maximum sensitivity — catches any change to any weight. But too sensitive: even 1 bit of difference in 1 weight (inevitable across BLAS implementations) changes the hash. Will produce false positives on any cross-machine comparison. No diagnostics on failure. Reject.

B. Full-vector hash at rounded precision (e.g., 3 significant figures)
Round each weight to N significant figures, then hash. Absorbs small differences. But suffers from the rounding boundary problem: a weight of 2.40500001 rounds up on one machine, 2.40499999 rounds down on another. With 215K weights per area and potentially 3,000+ areas, boundary collisions become likely. Reject.

C. Relative tolerance on weight sums only
Compare unrounded weight sums per area using np.isclose(rtol=1e-3). Solves the boundary problem and handles small areas correctly. But still distribution-blind — a complete redistribution of weight across records is undetected if the sum is preserved. Insufficient alone.

D. Multiple aggregate statistics per area with relative tolerance ✅
Compare several statistics per area, each with np.isclose(rtol=1e-3). No rounding boundaries, distribution-sensitive, diagnostic on failure, fast. Recommended.

Recommended design

Statistics per area (computed on the WT{TAXYEAR} column):

Statistic	What it catches	Match type
n_records	Data pipeline changes	Exact integer
weight_sum	Overall level shift (1st moment)	rtol=1e-3
weight_std	Distribution shape change (2nd moment)	rtol=1e-3
n_positive (weight > 0.01)	Sparsity pattern changes	Exact integer
max_weight	Outlier / extreme weight behavior	rtol=1e-3

Why these statistics:

• Most of the ~215K records per area have near-zero weight (only records "in" the area carry meaningful weight), so percentiles (p25, p50, p75) would be uninformative for most areas
• weight_std captures distribution shape and is more interpretable in failure messages than sum_of_squares (they carry equivalent information given weight_sum and n_records)
• n_positive catches the solver zeroing out different records across machines — a structural difference that sum/std might miss
• max_weight catches pathological solver behavior at the tail

Why rtol=1e-3:

• Cross-machine noise floor: O(1e-6) to O(1e-4) relative in aggregate statistics
• Real code/data changes: typically > 1% (O(1e-2)) relative
• 1e-3 sits cleanly between noise and signal with an order of magnitude margin on each side
• With 3,000 areas × 3 float statistics = 9,000 comparisons, the probability of a false positive remains negligible at this tolerance

Storage format: One JSON file per area type (states, CDs, future counties):

{
  "metadata": {"taxyear": 2022, "n_areas": 51, "generated": "2025-03-30"},
  "areas": {
    "AK": {"n_records": 215494, "weight_sum": 733583.2, "weight_std": 3.41, "n_positive": 4521, "max_weight": 8.34},
    "AL": {}
  }
}

Test structure:

Overall: all expected areas present, n_records matches exactly
Per-area float stats: compare with rtol=1e-3, report which areas and statistics diverged and by how much
Scalability: Adding county weights (3,000+ areas) just means a larger JSON file and more iterations in the comparison loop. Computation per area is O(n) with a small constant (one pass through the weight vector). The dominant cost is I/O (decompressing gzipped CSVs), not the fingerprint math.

Note on the existing national weight fingerprint

test_weights.py already uses a multi-statistic approach (n, total, mean, sdev, min, p25, p50, p75, max, sum_sq) with np.allclose defaults (rtol=1e-5). This is sound for the national weights (single test, well-conditioned). The area fingerprint design above adapts this pattern with a relaxed tolerance appropriate for many areas and cross-machine use.

[donboyd5]

@martinholmer, if you think the proposed improvement above makes sense, I can prepare a draft PR that you can check out on your machine. Does that make sense?

[martinholmer]

@donboyd5, yes I think your approach makes sense. Thanks for the careful consideration of this issue.