Reach the vulnerability domain the splines claim

[dfalster]

The knot grid accumulated psi += step, so it yielded resolution or resolution - 1 knots depending on b and c, and its upper bound was neither reproducible nor equal to vulnerability_psi_max. Both splines built from it disable extrapolation, so that edge decides whether a dry-end lookup throws. At the package defaults it fell one full step (1.0%) short. Indexed now (#92).

Moves results: 3496 of 5184 golden cells, median 2.0e-15, largest absolute move 1.97e-10; recorded gradients 63 of 100, worst 5.1e-04. Nothing reaches the ~1e-4 band. Split by cause in the first comment.

With the edge well defined, psi_crit is checked against it by name rather than failing inside the interpolator (#38) — it is P95 of the stem curve, not a free trait, and that was written down nowhere. The thirteen traits are readable from R now, read-only (#95). Closes the stale note on the temperature cache key (#55).

Closes #92, closes #38, closes #95, closes #55.

[dfalster]

The defect, reproduced

for (double psi = step; psi <= psi_max; psi += step) over b in [1, 6] and c in {0.4, 1, 2.04, 2.680147, 6, 12}:

requested resolution	knots obtained
50	49 or 50
100	99 or 100
200	199 or 200
1000	999 or 1000

⚠️ One correction to #38's own figures. The issue reported the plant TF24 pair (b = 3.457267907, c = 2.04) as dropping a full knot, with psi_max = 7.3089243750824018. On this machine that pair keeps all 100 knots and psi_max = 7.3089243744058408 — the issue's psi_max differs in the 10th digit, so it was computed from a slightly different b. The pair still ends 1.6e-16 relative short of psi_max, so a lookup at exactly vulnerability_psi_max still throws; the reported 1.0% shortfall belongs to a different pair.

The pair that does drop a knot is this package's own default:

b = 3.898245, c = 2.680147, resolution = 100
old:  99 knots, last at 6.822923224755713
new: 100 knots, last at 6.8918416411673782 == vulnerability_psi_max(b, c)
     one full step, 1.0000% of the intended domain, on BOTH the stem and the root curve

Worst interior knot displacement at the defaults is 1.22e-15 relative, at i = 95 — so the results move below comes almost entirely from the added knot, not from the reindexing.

Blast radius: tests/cpp/golden/operating_points.tsv

3496 of 5184 finite cells move. Median relative move 2.0e-15 — reassociation band.

relative move	cells
< 1e-15	483
1e-15 – 1e-13	2177
1e-13 – 1e-11	427
1e-11 – 1e-9	284
1e-9 – 1e-7	54
1e-7 – 1e-5	71
> 1e-5	0

The 71 cells above 1e-7 are the ones worth explaining, and they are a small-denominator artefact rather than a model change: all 71 sit at psi_soil = 3 with 5 layers, where the quantities themselves run from 1e-11 to 1e-3. The largest absolute move anywhere in the file is 1.97e-10; the median absolute move among those 71 is 1.0e-15. Worst relative is 5.65e-06, on assim = 1.53e-06.

Per field:

field	moved	worst rel	median rel (moved)
psi_stem	421	3.19e-09	2.70e-15
opt_root_psi	372	5.02e-10	2.25e-15
ci	352	1.08e-09	3.37e-15
assim	342	5.65e-06	1.46e-14
transpiration	363	3.26e-06	1.39e-14
gc	363	3.26e-06	1.39e-14
profit	408	1.82e-07	2.01e-15
e_up	437	3.01e-06	1.37e-14
uptake	438	3.01e-06	1.32e-14

Point classification is unchanged: 198 interior / 42 pinned / 48 shutdown at 25 °C, 160 / 80 / 48 at 40 °C. No point changed branch.

Blast radius: tests/testthat/gradient_golden.tsv

63 of 100 cells, median 1.1e-08, worst 5.15e-04. Larger than the solved outputs, and it has to be: these are finite differences of the same solve, so they carry its floor divided by a ~1e-6 relative step. The file's own cross-platform tolerance is 5e-3, so the worst move sits an order inside the noise band it is already compared with elsewhere.

By case: interior-1layer 5.15e-04, pinned-dry-3layer 1.20e-04, single-potential 9.0e-08, interior-5layer 5.3e-08, shutdown-1layer 0 (it never reaches the collar solve).

The largest single cell, in absolute terms: interior-1layer / vcmax_25 / profit, 0.01052672793 → 0.01053214776, i.e. 5.42e-06 absolute.

The four R-side hex rows in test-golden.R moved 27 of 36 values, scaling with how dry the point is — 1e-15 at psi_soil 0.5 and 2, 1e-11 to 3e-11 at 4, and the shut-down row's nine unchanged.

⚠️ Two things this went wrong on first, recorded so they are not repeated

1. tools/gradient_golden.R regenerated the file from the SITE build. It does library(phylloptim), which resolves against the library path and not the tree you are standing in. Run from this worktree with uninstalled changes, it reproduced the old values exactly and reported "0 of 100 cells moved" against a suite that was failing on 63 of them. Same trap tools/bench_history.sh guards with --expect-version, and the same one the guide records for plant.

The script now refuses to run unless the loaded build's headers match inst/include/ byte-for-byte — R CMD INSTALL copies them into the package directory, so that is a direct answer to "is the loaded build this tree?", where comparing paths cannot work (an installed package is never in the tree) and comparing versions cannot either (0.2.1 both sides). Verified in both directions: refuses against the site build, runs against the installed tree and regenerates the committed file byte-identically.

2. A magnitude I asserted in a comment before measuring it. The test-golden.R header note first said the 27 moved values were "all at or below 3.5e-13 relative except assim/gc/transpiration at the 4-layer row". They are not — the split is by row, not by column. Measured and corrected before commit; flagging it because a number in a comment beside a regenerated file is exactly the kind of claim that is never re-checked.

What #38's check found inside this repo

The check rejects only combinations that already threw, so nothing that worked stops working — but two places in the tree held such a combination and had never solved through it:

• tests/testthat/test-surface.R built leaf_model(leaf_traits(stem_b = 2.0)), leaving the default psi_crit of 5.870283 against a P99 of 3.5359. It passed because that test only calls proportion_of_conductivity, a closed form that never reads the spline. psi_crit now moves with stem_b there.
• ?leaf_traits's example was leaf_traits(vcmax_25 = 120, stem_b = 2.5) — P99 = 4.42, so a reader copying it into leaf_model() got the error the issue is about. Fixed, and ?leaf_traits now states the P95 relationship.

test-gradient.R already worked around this at line 653 ("stem_b is raised, not lowered. #38: …"), which is corroboration that the constraint was real and being managed by hand.

The one test threshold that moved, and why it is not a loosening

test-gradient.R's expect_lt(abs(exact), 1e-12) now reads 1e-9. exact went from 2.4e-15 to 1.2e-10 on macOS/arm64 — the knot grid moved where the collar root-find lands inside its own tolerance, which is precisely the mechanism the comment above that assertion already documents for Linux (it records exact = 6.6e-11 there). 1e-12 was pinning a landing point; the guide names ~1e-9 as the solver floor, which is what the bound should have been.

What carries the test is the gap, and it is intact: exact 1.2e-10 against fd 2.1e-04 is a factor of 1.7e+06, with both bounds three or more orders inside it. Removing the assignment the test guards still fails it.

What #38 assumed and I did not do

#38 says "⚠️ Check the root curve too — MultiLayerRoots::setup_vulnerability uses the same builder, so root_psi_crit presumably has the same latent constraint."

It does not, and #77 is why. root_vuln_from_psi has extrapolation disabled like the stem pair, but root_vuln_at clamps its argument to the last knot rather than letting the read fail. So a root_psi_crit past the root P99 silently reports the floor conductivity instead of throwing — a different defect, and the conductivity half of it is #85's open question. Adding a hard error there would reject configurations that work today, so the check is stem-only and says so at the declaration.

Verification

• make -C tests/cpp — 503 checks, 0 failures; golden bit-identical after the deliberate regeneration
• make -C tests/cpp bench_solve bench_gradient — both build (the targets CI builds and make all does not)
• make test-cmake — 2/2 ctest, the no-R route a C++ consumer takes
• devtools::test() — 1096 pass, 0 fail, 0 skip; gradient_golden.tsv: worst relative difference 0
• The new tests were run against the old builder to confirm they have teeth: 712 of 1800 grids wrong count, 1577 ending short (worst 0.100000 relative, at resolution 10 where one dropped knot is 10%).
• CI's inline consumer program constructs a default phylloptim::Leaf, so psi_crit above the spline domain fails with an error naming neither psi_crit nor the curve #38's check cannot break it.

The knot-count sweep asserts x.back() == vulnerability_psi_max(b, c) bit-exactly, not with a tolerance. That is deliberate: the whole defect was an endpoint that was close and not equal, and with extrapolation disabled one ULP short is the difference between a lookup at vulnerability_psi_max working and throwing. A tolerance there would pass on the code the test exists to reject.

[dfalster]

CI found one thing, and it was not the change

All seven jobs are green now, on both compilers and all three platforms. Getting there took one fix worth recording, because it looked like a results problem and was a measurement problem.

ubuntu-latest (release) failed test-surface.R's check that Leaf::operating_point_values() — one flat vector across the R boundary — agrees with reading the twelve outputs one binding at a time. macOS passed. Eleven of the twelve elements were bit-identical; element [10] differed by 1 ULP:

[10] 0.00000405328878168362639  -  0.00000405328878168362724

Element 10 is uptake, and it is the only one of the twelve that is not a binding — it is a sum over soil layers, so the R side of the comparison had to derive it, and it used sum().

R's sum() accumulates in LDOUBLE, which is 80-bit on x86-64 and 64-bit on arm64. operating_point_values() accumulates in double. So sum() agrees with it bit-for-bit on arm64 and can differ by an ULP on x86-64, for any sum of more than one term — which is why only the three-layer case failed and only on ubuntu. The one-layer cases have nothing to accumulate.

So the assertion had always been platform-dependent and had happened to hold. This PR moved the values and it stopped holding. Nothing about the code under test changed platform behaviour.

Fixed by making the comparison exact everywhere rather than relaxing it: Reduce("+", finite, 0) accumulates left to right in plain double addition, from the same zero the C++ loop starts at. Verified locally that sum(), Reduce and the C++ value are all three identical on arm64 — i.e. the change is a no-op on the platform that generated the golden files — and ubuntu now passes. A tolerance would also have made it green, and would have blunted the one thing the test is for: a shifted column, which is worth catching at the last bit.

No shipped code is affected. operating_point() reads the C++ value rather than recomputing the sum, so there was never a two-route disagreement in the package itself — only in the test's independent reading. I checked that before reaching for a fix, because "two implementations of one sum" would have been the more serious finding.

Recorded in the guide next to the existing note on R's decimal parser, since it is the same shape: as.numeric there, sum here, and both times the model was right and the instrument was not. sum, mean, cumsum and prod all use LDOUBLE; + on doubles does not.