Differentiate profit, the output plant bills

[dfalster]

The four differentiated outputs are what a gas-exchange calibration observes.
plant bills carbon from leaf.profit_ and water from soil_consumption_, so
the two sets were disjoint and no trait gradient here reached a demographic
model. Adds profit as a fifth column.

It is the one output the envelope theorem reaches: at an interior optimum
dprofit/dpsi = 0, so its gradient is the direct partial alone. The dropped
term is noise, not truncation — a central difference of a flat maximum — and
keeping it moves the answer by up to 8.0e-05.

The four existing columns are bit-identical. Closes #87.

[dfalster]

What changed

n_outputs 4 → 5, output_names() gains "profit", outputs() gains
y[4] = l.profit_, and R's .gradient_output_names matches. Both sides append
rather than reorder, for the reason par_names already records — R indexes these
positions by integer.

The envelope is applied as one assignment at the site that knows status,
not as a special case inside the per-parameter loop:

dY_dpsi <- (hi - lo) / (2 * h_psi)
if (identical(status, "interior")) dY_dpsi[["profit"]] <- 0

gradient_ift then computes direct + rounded(0.0 * dpsi_dtheta), which is
bit-identically direct, so no new product and no new FMA hazard. It also puts
the fact where it belongs: dY_dpsi[["profit"]] is an estimate of
∂profit/∂ψ, and the envelope theorem is the statement that it is zero.

The measurement, which is the reason for the assignment

Over the golden grid's 136 interior rows (192 points scanned: 136 interior, 24
pinned, 32 no-gradient), comparing the two available estimates of ∂profit/∂ψ
at ψ*:

	min	median	max
|dprofit/dpsi|, exact forward AD	1.7e-18	4.9e-15	5.4e-10
|dprofit/dpsi|, central difference	0	7.8e-10	2.1e-04
ratio fd / exact	0	1.4e+05	8.7e+10
relative move in dprofit/dθ if kept	0	3.3e-10	8.0e-05

The dropped term is noise, not an h² truncation, and that distinction is
the whole argument. profit is the maximum, so it is flat; a central difference
of it divides the solve's ~1e-09 floor by a ~1e-06 step, which lands at ~1e-04.
This is hazard "profit is the wrong instrument for checking a collar-solve
change" from the dev guide, arriving from a new direction — and here the package
already owns the right instrument, so the comparison is available for free.

The worst row (8.0e-05) sits in the band this repo calls a real difference rather
than rounding, so this is worth doing and not merely tidy.

Conditional on status, not on use_ift

The identity comes from dprofit/dpsi == 0, which is what Interior means. At a
pinned optimum ψ* is a trait-dependent bound, dprofit/dpsi is not zero, and
the indirect term is real — so profit takes the same finite-difference fallback
as the other four columns. Someone who forces method = "ift" at a pinned point
already gets a confidently wrong number; they should not get a differently
wrong one for this column alone.

Visible in the regenerated golden file: psi_crit's profit gradient is exactly
0x0p+0 on the two interior rows and 0x1.03cb5d0db7589p-1 on the pinned one,
which is the same sharp statement the file already made for A.

Verification

• tests/testthat/gradient_golden.tsv: the four existing columns are
  bit-identical. Generated to a temp file and diffed on cut -f1-8 before
  installing the new one — no existing cell moved.
• R suite: 0 failures, 1 CRAN-only skip. gradient_golden.tsv: worst relative difference 0 (bit-exact platform) — so the R and C++ composites agree
  bit-for-bit on all five columns.
• C++ suite: 463 checks, 0 failures; 288 golden operating points bit-identical
  (198 interior / 42 pinned / 48 shutdown, unchanged). make bench also builds.

⚠️ The new test was written twice, and the first version was worthless

test_that("profit's gradient is the direct term alone at an interior optimum")
originally used the suite's usual grid_drivers(2.0). At that point the dropped
term is 2.9e-10 relative — so the test passed whether or not the zeroing was
there, which is the only property a test like this must not have.

It now runs at psi_soil = 0.5, vpd = 2, 3 layers, the worst of the 136 interior
rows. Measured there: reported 0.0102951008, with the term kept
0.01029462274, a relative difference of 4.6e-05 against a tolerance of
1e-7 — so removing the zeroing fails it by 464×. The comment records both the
chosen point and the rejected one.

(The first draft also had a genuine bug: it read hi/lo off a leaf left
carrying a perturbed trait from the preceding difference, which is why it
reported 1.8e-4 at a point where the true figure is 1.4e-9. Both instruments are
now read on a leaf at base traits.)

What was considered and is not here

uptake. gradient.hpp's outputs() reads members directly while R goes
through operating_point_values(), and the header notes that the three columns R
computes rather than copies are not among the outputs. uptake is one of those
three — a sum over the finite soil layers — so adding it means reproducing that
summation, and its order, on the C++ side. Small, but a different decision from
this one, and plant's water budget is the motivation rather than an afterthought.

lambda (dA/dE). Same class, also a computed column, and it is #3's and the
companion manuscript's quantity. It needs its own motivation.

Per-layer soil_consumption_. A different shape; not a column.

Follow-ups: #88 (a gradient at a ψ the caller sets) and traitecoevo/plant#614
(TF24f's dψ/dθ companion state) — this PR is what gives either of them
something on plant's side of the boundary to propagate into.

[aornugent]

The profit column is the right thing to add and the envelope argument is the right
way to get it — carbon reaching plant's mass budget is leaf.profit_, so a gradient
set that stops at assim_colimited_ reaches no demographic model. The measurement
justifying the zeroing (eleven orders between the exact and differenced instruments
at the median) is the kind of evidence this decision needed.

Requesting changes on five points. The first two are rebase blockers rather than
design objections.

1. The golden file is stale, and it is not a text conflict

26574f4 (#86, R_d rises with temperature) landed after this branched and moved every
recorded cell. interior-1layer / vcmax_25 / A:

master   0x1.56ff4d6eecf96p-6
this PR  0x1.19f24c79f383p-6

So "the four existing columns are bit-identical" holds against the branch point and
not against master. Rebasing means regenerating tests/testthat/gradient_golden.tsv
in full via tools/gradient_golden.R — and per .claude/CLAUDE.md that has to happen
on macOS/arm64, so it cannot be done as part of a conflict resolution on any other
box. Worth doing before the rest of the review, since the regenerated numbers are what
the profit column's own values will be measured against.

2. helper-golden.R reopens a question master has already closed, at a number master's own CI says is too tight

Master added gradient_golden_tolerance() <- 5.0e-3 and a tolerance = parameter on
expect_golden(). This PR adds golden_tolerance(field, kind) returning 1e-3 for the
gradient file. The two are the same fix by different routes, but the numbers disagree
and master's is backed by the tighter measurement: it is not that the gradients are
argmax-evaluated, it is that the tolerance is set by the smallest differentiated
parameter's relative step. R_d_25 is 1.44 where vcmax_25 is 96, and master
records dA/dR_d_25 disagreeing by 1.3e-03 on Linux CI. A rebased version of this
PR carrying kind = "gradient" → 1e-3 would fail on that column.

Suggest dropping the kind argument entirely and using master's
tolerance = gradient_golden_tolerance() at the call site. It is also the more general
shape — per-call-site tolerance rather than a two-valued enum that has to grow a member
per file.

3. Two hand-maintained copies of the output names

.gradient_output_names in R/gradient.R and output_names() in gradient.hpp are
now both literals, each carrying a ⚠️ comment saying appending is safe and reordering
is not. That comment is describing a hazard that does not have to exist: par_names
has the same problem and solves it by having test-gradient-batch.R read the C++ list
back and compare, and the output list can go one better by having R read it rather than
restate it — export a gradient_output_names() alongside the existing
gradient_par_names() and make .gradient_output_names a first-call cache over it,
like .gradient_outputs_idx() already is. Then adding a sixth output is one edit
instead of two that can disagree.

4. direct + rounded(0 * dpsi_dtheta) rather than assigning the column

In gradient_ift, profit's answer is reached by zeroing dY_dpsi[out_profit] and
letting the generic line multiply it by dpsi_dtheta. collar in the same loop is
handled the other way — assigned explicitly, with a comment saying that beats relying
on a difference of two identical numbers. The same argument applies here and buys more:
0 * d_psi is NaN for any non-finite d_psi, and #88 makes d_psi caller-supplied
(and .gradient_dpsi_dtheta checks is.numeric && !anyNA, not is.finite), so a
caller passing Inf gets ±Inf in four columns and NaN in profit alone. Setting
out[k * n_outputs + out_profit] = direct[out_profit] is immune and says the thing the
comment above it is already saying.

5. The shut-down rows record a number nothing asserts

shutdown-1layer / stem_b records dprofit/dstem_b ≈ 2.13 where A, gc, psi_stem
and collar are all exactly zero, at a row whose own status is no-gradient. As a
drift detector that is fine; as a baseline it pins a large value with nothing saying
what it ought to be, in the one regime where profit_ comes from a branch that sets it
without the other outputs (hazard 8). Worth an explicit assertion of the shut-down
profit column against its closed form, rather than only the recorded hex.

Minor

The envelope paragraph and its three-row table now appear in NEWS.md, the roxygen
block, the generated .Rd, README.md, two vignettes, and inline in both R/gradient.R
and gradient.hpp — seven coordinated edits if the decision ever changes. CLAUDE.md
makes this argument about commit bodies ("that is the habit to drop"); it applies to
source comments too. The inline C++ and R comments could be the decision plus a pointer
to ?leaf_gradient without losing anything.

[dfalster]

All five addressed, pushed as e1d8249. Taking them in order.

1. Stale golden file — done, and the check is now against master

Rebased onto 26574f4 and regenerated gradient_golden.tsv in full on
macOS/arm64. The bit-identical claim is now made the way you asked for it:

Rscript tools/gradient_golden.R > /tmp/new.tsv
cut -f1-8 /tmp/new.tsv | diff - <(git show origin/master:tests/testthat/gradient_golden.tsv)
→ no differences

So the four pre-existing columns are bit-identical to master, not to the
branch point, across all twenty rows including the R_d_25 ones #86 added.

⚠️ I also re-took the envelope measurement, which the rebase invalidated.
#84 and #86 both moved the solve, and the table in ?leaf_gradient was measured
before either. Re-run on the rebased build: 136 interior rows still, exact
median 4.9e-15 → 4.8e-15, relative move median 3.3e-10 → 2.7e-10, and
max 8.0e-05 and the worst operating point both unchanged. The two medians are
updated. Had I not re-measured, the docs would have carried numbers from a
different solver — which is the failure mode this repo already has a name for.

2. helper-golden.R — dropped entirely, master's design adopted

You are right that these are the same fix and that master's number is the
better-founded one. Mine reasoned from the argmax amplification, which the
gradient columns do inherit but which is not what binds; yours reasons from the
smallest differentiated parameter's step, which is, and it has dA/dR_d_25 at
1.3e-03 behind it. My 1e-3 would indeed have failed on that column after the
rebase.

golden_tolerance(field, kind) is gone. The call site is
expect_golden(..., tolerance = gradient_golden_tolerance()), and I agree
per-call-site beats an enum that grows a member per file.

3. Two copies of the output names — now one, and R reads it

gradient_output_names() is exported alongside gradient_par_names(), and
.gradient_output_names is a first-call cache over it.

Worth recording why the two enumerations differ, since the asymmetry now looks
arbitrary: par_names has to stay a second copy because R builds theta from
leaf_traits() before any C++ call, so a test comparing them is the best
available. The output list has no such constraint. The hazard comment is deleted
rather than reworded — there is no longer a hazard to describe.

4. direct + rounded(0 * dpsi_dtheta) — assigned instead

gradient_ift / .gradient_ift() take an envelope flag and assign
out[k * n_outputs + out_profit] = direct[out_profit], exactly as collar is
assigned. dY_dpsi[out_profit] is no longer zeroed at all; it now carries the
exact dprofit/dpsi (the AD resid, already computed), which is strictly
better than a difference of the same quantity and is also what #88 wants to
return.

The NaN case you identified is real and this closes it — with d_psi = Inf the
profit column now comes back as the direct term rather than NaN. Bit-identical
on the existing grid, since direct + rounded(0 * finite) was already exactly
direct.

5. Shut-down row — asserted, and it has two closed forms rather than none

Good catch, and it turned out better than a spot check. At that row E = 0, so:

• A = -R_d exactly and R_d(25 °C) = R_d_25, and the hydraulic cost does not
  depend on R_d_25 — so dprofit/dR_d_25 = -1. Recorded value
  -0x1.00000000636f2p+0, i.e. −1 + 5.8e-09, the solver floor.
• the shut-down collar is pinned at psi_crit, so dcollar/dpsi_crit = 1
  (recorded 0x1.000000001aa4dp+0) — which is the whole explanation of why
  psi_crit alone carries a non-zero profit gradient there. It moves the
  collar, and the collar sets the cost.

Both are now assertions in the existing shut-down test, alongside
profit == A - hydraulic_cost read off a solved leaf. The stem_b value you
quoted (2.13 then, 1.129 after the rebase) is a branch change under
perturbation and stays a recorded number — it is not a closed form, and I have
not pretended otherwise.

Minor — trimmed

The table now appears in three places, each with a different job: NEWS.md
(release note), the roxygen (canonical, ?leaf_gradient), and the test (which
needs only the one figure that justifies its operating point, so it keeps
8.0e-05 and points at the roxygen for the rest). Both inline comments are the
decision plus a pointer, as you suggested. README and the two vignettes never
had the table.

---

Local: R suite 0 failures, gradient_golden.tsv: worst relative difference 0,
C++ 486 checks / 0 failures, 576 golden points bit-identical, make bench
builds.

Thanks, requested changes have been implemented