Stop formulas writing into arrays they read from other variables

[MaxGhenis]

Summary

policyengine-core returns the cached array itself from tax_unit("x", period) (Simulation._calculate returns holder.get_array(...)), and add(entity, period, sources) returns that same array when sources holds a single same-entity variable (for_each_variable starts from the first variable's values). A formula that writes into it (x += y) therefore changes variable x's cached value for every later reader in the simulation, and inside a branch it changes the branch's copy. Five formulas did this:

formula	variable it overwrote	when
reforms/crfb/agi_surtax.py	adjusted_gross_income (+ the increased-base sources)	CRFB surtax with increased_base.in_effect
ga_deductions.py	federal itemized_taxable_income_deductions (− the Georgia adjustment)	Georgia, adjustment ≠ 0
ut_personal_exemption.py	ut_total_dependents (+ newborn dependents)	Utah, from 2023
ut_dependent_exemption_reform.py	ut_total_dependents	reform installed but not in effect that year
basic_income_phase_in.py	tax_unit_earned_income (+ Social Security)	include_ss_benefits_as_earnings

Each now builds a new array (x = x + y). What changes is every variable calculated after them from the overwritten one. The formulas' own outputs are unchanged with one exception: agi_surtax now sums its base in float64. The old agi += … rounded AGI plus the extra sources to float32 before applying the rates (float32 spacing is $1 at $10 million and grows with AGI). In the 2026 run this moves 3 records by at most $0.02. For a single filer the gap grows with AGI: the review's NumPy emulation of the schedule puts it at up to $0.125 for AGI of $10m–$100m and $2 for $100m–$1bn.

nj_gross_income and ma_ccfa_countable_income_person wrote into the result of add() over a parameter list. Those lists hold 13 and 14+ variables in every dated value, so the write never reached the cache. They now build new arrays too, so a later edit of either list to a single variable cannot turn them into cache writes. With the current lists their outputs don't change (measured below). With a one-variable MA earned list, the old subtraction also overwrote the formula's own earned local (the same cached array). The new code gives the formula's result instead. A test checks that the source variable is no longer overwritten; it does not pin the formula's result. Separately, if a reform turned on the minor-earnings exclusion with parent-only counting off, the dependent path would subtract a minor's earnings twice. Baseline never combines the two (the exclusion starts in 2026; parent-only counting has applied since October 2023), and this PR doesn't change it: #9759.

Tests

• tests/test_formulas_do_not_write_into_cached_arrays.py (17 tests). One test per site calculates the consumer, then checks the source variable: AGI stays $400,000 after income_tax under the CRFB reform; New York state income tax under the reform equals baseline; Georgia's federal itemized deductions are unchanged after household_net_income; ut_total_dependents stays 2 for a newborn plus a child, in baseline and under the dependent-exemption reform; earned income stays $10,000 under the basic-income phase-in; MA's self_employment_income stays put with a one-variable earned list. A final test calculates three kinds of household in every state and DC under baseline, the CRFB reform, the Utah reform and the basic income. Every cached array is read-only, both as stored and as read, which also covers the arrays branches inherit. Any formula on that path that writes into a cached array then raises at that line. A further test checks that this fixture catches a write into a branch's inherited array. The per-site tests are parametrized over each formula's branches: Utah's newborn credit on and off, Social Security counted or not, per-person phase-in, and MA's parent-only rule and minor exclusion. That way every line this PR changes runs, along with every branch of ut_personal_exemption, basic_income_phase_in and ma_ccfa_countable_income_person. The Utah reform's in-effect branch, which this PR doesn't touch, is not exercised here. 13 of the 17 fail on main; all pass here. The 4 that pass on main are the three paths where the old code never wrote, plus the fixture test.
• tests/code_health/test_no_in_place_writes_to_cached_arrays.py: an AST scan of variables/ and reforms/ that fails CI on new writes into arrays read from variables. It catches augmented and subscript assignment, np.place/putmask/copyto, ufunc.at, out=, and in-place methods such as .fill() and .sort(), including writes through views such as x.reshape(-1)[:] += 1. It tracks the result of entity calls (with a literal or variable name), .members(...), .calculate(...), and add() over anything but a literal list of two or more variables, through aliases, views, conditional expressions and tuple unpacking. The sides of an if, match cases and try blocks start from the same state and merge, and loop bodies are walked as if they ran zero, one or two times. So an optional copy on one path doesn't hide a write on another. On main it reports the five sites plus the four NJ lines and the MA line; here it reports none, across the whole package. It is not exhaustive. Positional out, alias chains that need three or more loop passes, walrus bindings, closures and some views escape it. None of those appears in the codebase; Harden the in-place cached-array write scan (positional out, loop fixed point, walrus, more views, tools/) #9760 tracks them. It has its own positive and negative self-tests. The scan from the original audit also flagged three_digit_zip_code.py, which writes into a new pandas Series; this one doesn't.
• With every cached array read-only, as stored and as read (so branch-inherited arrays count), 2,195 existing YAML tests pass and no write was logged: GA 413, UT 393, NJ 643, MA 651, contrib CRFB 35, contrib Utah 40, UBI Center 20. The patch was a sitecustomize.py on InMemoryStorage.put and .get that logs any read-only ValueError raised inside a formula.

What changes (real microsimulation runs)

Default dataset (populace_us_2024.h5@populace-us-2024-spm-20260909), 2026, one Microsimulation per run. Each run calculates household_net_income first, then reads the other variables from the same simulation, comparing main (20ccd5ac7c) with this branch. Scripts and outputs: impact.py, compare.py, out/compare.md.

Baseline (current law). Every tax, benefit and net-income total is identical. The only output that moves is the reported ut_total_dependents: 1.165m on main vs 1.127m here. On main, newborn dependents were added into it a second time; the Utah exemption itself was always right. 41 records (0.04m weighted tax units) differ. Georgia's adjustment is an input the dataset does not carry (zero for every record), so no Georgia output moves in microsimulation. The household case is fixed and tested.

CRFB AGI surtax with the increased base (reform minus baseline):

	on main	fixed	correction
household net income	−$111.89bn	−$80.49bn	+$31.40bn
federal income tax	+$82.31bn	+$80.52bn	−$1.79bn
income tax before credits (≈ the surtax)	+$80.52bn	+$80.52bn	0
state income tax	+$29.28bn	−$0.03bn	−$29.31bn
premium tax credit	−$6.55bn	$0.00bn	+$6.55bn

The surtax itself is unchanged. On main, everything calculated from AGI after agi_surtax used AGI plus the increased-base sources (AGI totalled $20,400bn on main vs $18,931bn). That covers state income taxes, federal credits calculated after the surtax, and the premium tax credit. The state correction is largest in CA ($7.66bn), IL ($3.09bn), MD ($1.84bn), MN ($1.83bn) and NC ($1.73bn). On main the reform appeared to cut net income in the bottom half of the distribution ($21 to $224 per household in deciles 1–5). Fixed, it is about zero there, as the surtax starts at $100,000 / $200,000.

Utah dependent exemption reform, installed but not in effect in 2026 (in_effect from 2027, so the 2026 run uses the reform's formula, checked): same as baseline. Only ut_total_dependents moves (1.165m → 1.127m); no tax or income output changes.

Basic income ($12,000 flat) phased in with Social Security counted as earnings: the basic income itself is unchanged. On main the reported tax_unit_earned_income total was $1,677.9bn too high (39.95m weighted tax units differ). One later reader picked it up: state income tax falls $0.011bn here (8 records), and household net income rises by the same amount.

Runs of the CRFB increased-base option before this fix (it was added 2026-01-23) had the same bug: everything calculated from AGI after the surtax used the inflated AGI, as in the main-branch column above.

Invariants

• Calculating any variable leaves every other variable's cached value unchanged, so results don't depend on calculation order. The read-only-cache test checks this on the paths it runs, and the code-health scan checks it statically within its limits (Harden the in-place cached-array write scan (positional out, loop fixed point, walrus, more views, tools/) #9760).
• Each changed formula returns the same value as before for every input with the current parameters, except agi_surtax's float32-to-float64 rounding (see above). An independent review executed old and new formula versions on vectors including negatives, cents, large float32 values and int32 counts, with identical outputs for GA, both Utah sites, the basic income, NJ and MA. The per-site tests check the formula outputs (agi_surtax, ga_deductions, ut_personal_exemption, basic_income_phase_in) alongside the sources.

axiom: n/a: fixes cache aliasing in formulas; no rule change

Follow-ups: #9759 (MA CCFA double subtraction, pre-existing), #9760 (scan gaps).

Reviews: two rounds by an independent agent (GPT-6.1 Sol via Subfleet). Round 1 requested changes; they were fixed in fd20f55 or answered in the response. Round 2 approved the head 95d9736.

Related: PolicyEngine/policyengine-core#556 (branches share cached arrays and copy on first read) deliberately keeps today's behaviour for these formulas, so this fix is independent of it.

🤖 Generated with Claude Code

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (909176a) to head (95d9736).
⚠️ Report is 118 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff            @@
##              main     #9739   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files            4         5    +1     
  Lines           76        90   +14     
  Branches         2         5    +3     
=========================================
+ Hits            76        90   +14     

Flag	Coverage Δ
unittests	100.00% <100.00%> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[MaxGhenis]

Merged with --admin --match-head-commit 95d97361 on Max's ruling (decision d791: merge on gates). Gates at that head: every check passing (gh pr checks exit 0, 37/37), mergeable, not a draft, no changes requested. Independent review: two rounds by GPT-6.1 Sol via Subfleet. Round 1 (jobs 20261001-113258) requested changes, fixed in fd20f55 and 95d9736. Round 2 (20261001-135140) approved head 95d9736; its PR-text findings are applied above, and its scan gaps are #9760. Pre-existing MA issue found in review: #9759.