R_d rises with temperature, and is settable (#41)

[dfalster]

R_d was rd_to_vcmax_ratio_ * vcmax_(T), so it inherited Vcmax's PEAKED
Arrhenius and FELL above the thermal optimum, where real dark respiration
rises. No value of the ratio fixes a function of the wrong shape. Now
Tjoelker's declining Q10, with R_d_25 a plain trait -- default 1.44, the
same double 0.015 * 96 gave, so the defaults are unchanged.

Respiration no longer follows vcmax_25, and results move at every
temperature but 25 C, where agreement is exact by construction. plant's SCM
is bit-identical at its defaults. The golden grid gains a 40 C block because
it was blind to all of this. Measurements and traps in the comments.

[dfalster]

The measurements

Blast radius of the shape change, at the defaults

T	R_d old → new	A old → new
15 °C	0.669 → 0.646	6.1505 → 6.1572	+0.11%
25 °C	1.440 → 1.440	5.5992 → 5.5992	exact
35 °C	1.610 → 2.592	3.4924 → 3.1550	−9.66%
40 °C	1.012 → 3.171	1.8265 → 0.8555	−53.16%
45 °C	0.507 → 3.618	0.6339 → −3.6176	shut-down

The 25 °C row is exact rather than approximate: the reference value defaults to
rd_to_vcmax_ratio_ * vcmax_25, and vcmax_(25) == vcmax_25 by construction of the
peaked Arrhenius.

Three arms, over the 1152-point golden grid, against origin/master

Harness: tests/validate/rd_shape_grid.cpp (one tree, three binaries, per hazard 5).

arm	cells differing	what it isolates
default, 25 °C only	0 of 2592	the control
rd_tracks_vcmax_ = TRUE, all four temperatures	0 of 10368	the shape change alone
default, all four temperatures	6018 of 10368	the blast radius

The middle arm is the result of this PR. With the old shape restored, everything
else in #41 — the cache key, the settable R_d_25, the shut-down exit, the whole R
API — is jointly bit-identical to master. That is proved rather than argued, and it
is what the escape hatch was worth adding for.

Assimilation, relative change by temperature:

T	min	median	max
15 °C	+0.05%	+0.13%	+3.51%
25 °C	0	0	0
35 °C	−100%	−10.86%	−5.67%
40 °C	−213%	−62.35%	−34.95%

Per field, worst / median relative: assim 1 / 0.109, profit 1.91 / 0.082,
ci 0.615 / 0.078, psi_stem 0.484 / 0.022, the four flux fields 1 / 0.054.

No row changed feasibility — the 48 shut-down rows per temperature are the same
rows. Of the 864 rows whose A moves, 269 keep their exact operating point: 144
shut down (E = 0), and 125 sit pinned to a bound that temperature does not shift, so
only their carbon balance moves.

plant revalidation, deliberately not pinned at 25 °C

Harness: tests/validate/rd_shape_arms.R. TF24 SCM offspring production, one
species, max_patch_lifetime = 5 — the setup plant's own regression test uses.

plant configuration	phylloptim master	this branch	change
default (leaf_temp = 25, PM off)	81.857087216483691	81.857087216483691	bit-identical
leaf_temp driver = 30 °C	24.283073719078399	22.369992929645893	−7.9%
leaf_temp driver = 35 °C	2.1302e-04	6.4472e-07	÷331
Penman-Monteith on	2.2034981948510728	1.1638497249192270	−47.2%

⚠️ plant pins 25 °C itself, which is what makes the first row exact rather than
lucky: TF24_Environment sets leaf_temp as a constant extrinsic driver at 25 and
use_energy_balance defaults to 0. So a default plant run is unaffected, by
construction, and the exposure is entirely in the two configurations that leave 25 °C.

⚠️ The PM row is the one to carry forward. Leaf temperature on that path runs
19.1 to 62.0 °C — up to +22 K above air (over PAR 400–2000 × Tair 20–40 × VPD
1–3), so a PM leaf at 25 °C air already sits at 32.9 °C and loses 23% of its mean
assimilation. Tleaf itself moved on 22 of 90 cells: the energy balance couples
back through transpiration, so on that path leaf temperature is an output, not an
input.

plant needs no source change — it never calls set_traits, so the fourteenth
argument is invisible to it, and it binds only R_d_ of the respiration fields. Its
leaf/TF24 test files give 486 checks, 0 failures and one error, and that error is
present identically on the control, so it is pre-existing on plant develop.

The gradient reallocation, which is the one thing a caller must know

Making R_d_25 a parameter turns dY/dvcmax_25 from a total derivative into a
partial at fixed respiration. It used to carry ratio * dY/dR_d_25 as well.

Of the 60 cells in tests/testthat/gradient_golden.tsv, 14 moved and all 14 are
vcmax_25, by 16% to 250%; every stem_b, psi_crit,
leaf_specific_conductance_max and resistance cell is bit-identical, as is every
solved value. The 250% is the five-layer row, where the two terms nearly cancelled
(8.6e-06 total against 3.0e-05 partial), so dropping one is not a small edit.

The chain rule reconstructs all 14 old values — exactly where both arms are exact,
and to 2.8e-04 at worst where they are not. That residual is the direct term's own
finite-difference round-off, not a modelling discrepancy: it falls to 2.5e-06 at
step = 1e-4 and rises to 0.030 at step = 1e-8, i.e. coarser is better there.
Asserted as the identity rather than left as two unrelated recorded numbers.

Cleanest case, a shut-down point where A = -R_d exactly:

	before	now
dA/dvcmax_25	−0.015	exactly 0
dA/dR_d_25	(no such column)	−1, to 6e-11

Verified against a finite difference of the solve: dA/dR_d_25 = -0.24874 against
−0.24868 at one observation. ⚠️ Not ≈ −1, because the optimiser re-optimises the
operating point, so the total derivative is well below the direct −R_d term.

Cost

Cache key 3.03–3.04 µs/solve on both arms, interleaved ×3 from one tree, with
identical checksums — inside the ±0.01 within-process noise. The key is compared once
per set_physiology(), not per solve iteration.

What was rejected, and why

A constant Q10 = 2 was implemented first and rejected on measurement: R_d 4.07
at 40 °C (2.8× its 25 °C value over 15 K), −73% on assimilation, and no operating
point at all by 45 °C. That is not the physics; it is a constant Q10 used far
outside where it holds, which is what Tjoelker measured. Slope zero recovers a
constant Q10 exactly, and that is asserted — it is the escape hatch for anyone who
wants the conventional form.

⚠️ Slope zero does NOT recover the pre-#41 shape, which was peaked. That is why
rd_tracks_vcmax_ is a separate flag and not a parameter value.

Extending the generated Leaf constructor to take R_d_25, instead of assigning
the field after construction. Not taken: its argument list is positional, seventeen
long, and pinned by the raw-versus-friendly comparison in test-surface.R. The cost
is recorded — R_d_25 is now listed in the constructor-partition test as a trait the
constructor does not take, with a test that the assignment is really there.

Traps, each of which cost real time

⚠️ NaN == NaN is false. The temperature cache key is compared with ==, so the
raw sentinel must never go in it — it keys on rd_reference(), the resolved value.
Putting R_d_25 in directly silently disables the cache for every default caller, and
nothing notices.

⚠️ The parameter enumeration deliberately does not follow set_traits' argument
order. R_d_25 is set_traits' fourteenth argument and the enumeration's
sixteenth, appended after the two non-traits, because slotting it in at 13 would
displace par_kmax and par_resistance — the reorder gradient.hpp warns is unsafe.
Appending is safe; reordering differentiates the wrong parameter and returns plausible
numbers. test-gradient-batch.R asserts the divergence in both directions.

⚠️ .gradient_setter's positional trait call is where this broke. Its own comment
predicted that "if a trait is ever added, this call is the one place that will not fail
to compile" — and that is exactly how it went, at run time, with argument R_d_25 is missing raised inside the generated binding, naming neither the call nor the count.

⚠️ The batch path read R_d_25 as a missing trait, and the message was worse than
useless: `traits` is missing: NA, where the NA was the name of a trait that
does not exist. .gradient_theta_matrix() assumed the traits were the leading block
of par_names and sliced off "the last two", which now takes
leaf_specific_conductance_max. Both ends are placed by match() now.

⚠️ A theta column cannot hold the sentinel — and not only because of the validator.
A parameter whose value is "derive it" is not differentiable: perturbing NA gives NA.
Resolving it where theta is built is what makes R_d_25 differentiable at the
default, which is what a calibration starting from the defaults needs.

⚠️ The CI consumer is a fourth C++ program, in a heredoc in cpp-tests.yml, that no
local build covers. Its theta is brace-initialised, so appending to the enumeration
silently zero-filled the new element — a leaf with no respiration — and its arbitrated
dA/dvcmax_25 still carried the pre-#41 total derivative. Both fixed and verified by
extracting the two heredocs and building it locally: dA/dvcmax = 0.020934892290554417,
profit unchanged, exit 0, ctest 2/2.

⚠️ R_LIBS falls back to the site library in silence. It happened once here,
mid-install, and produced a plausible file describing the wrong package. Both new
harnesses hard-fail on the wrong library rather than warning, because every number
they print is plausible for either arm.

⚠️ use_energy_balance lives in plant's strategy$pars. Assigning it at the
strategy's top level is accepted and silently ignored, and the run comes back exactly
equal to the default arm — which reads as "PM changes nothing" rather than as a broken
switch. Asserted after setting.

⚠️ Doxygen 1.9 on CI rejects an indented block inside a function body, so the Q10
comparison table lives in the doc block above
update_temperature_dependent_params, not inline.

Corrections to the record

• The C++ suite's blast-radius table was labelled "old form vs Q10 = 2" — the
  rejected first implementation. Its numbers were always the shipped declining-Q10
  form. A printed label is asserted by nothing, so it outlived what it named. The
  "old" column is now the rd_tracks_vcmax_ branch itself, so the table is evidence
  the escape hatch reproduces what plant used to compute; the printed numbers are
  unchanged.
• vignettes/phylloptim.Rmd said "all 16" parameters when the answer was 14. That
  label is now computed from leaf_gradient() rather than written down.
• compare_with_plant.R's header blames the constructor signature for the harness
  being unrunnable. That is no longer the obstacle: against plant at 76df7169 both
  Leaf() and set_physiology() work. What actually stops it is root_collar_psi_
  versus opt_root_psi_ — the One representation for water potential: positive magnitudes throughout (supersedes the two-type approach in #23) #25 rename, which flipped the sign — on top of results
  deliberately moved by API cleanup: renames, a purely-intensive input set, the hidden hard-coded pressure, and the shutdown fix #15, 11a, 11b, One representation for water potential: positive magnitudes throughout (supersedes the two-type approach in #23) #25, Bound the root vulnerability curves past the grid #77 and Solve for leaf temperature inside the collar FOC #84. That is tests/validate/ has not compiled or run since #15 #64, and it needs a
  decision about the reference build, not a rename. Recorded so the next attempt
  starts from the real obstacle.
• bench_solve.cpp claimed to cover "the same state space the regression baseline
  covers". It deliberately stays at 288 points at 25 °C: a timing baseline is only
  useful against its own history (tools/cost-baseline.tsv,
  tools/bench_history.sh), so quadrupling it would end that history to measure
  nothing new about the solve.

The golden grid, and why it needed the axis

⚠️ The golden file was blind to every temperature response in the model, and had
been all along. 288 points, all at 25 °C — where the reference values of Vcmax,
Jmax and R_d are defined, so a change to any response curve is inert there by
construction. #41 changed R_d's shape at every temperature except 25 °C and the file
did not move a bit. The same blindness covered the thirteen temperature-response
parameters bound in the first commit.

Now 15 / 25 / 35 / 40 °C, 1152 rows. 40 °C is deliberate: it is where the response
bites hardest, and the shut-down rows at the hot end pin the temperature at which the
model stops having an operating point, which is a limit that moved here.

The regeneration is checkable as an addition rather than taken on trust.
Temperature is the outermost loop, so the file is four contiguous copies of the
original grid and the 25 °C block is byte-identical to the pre-#41 file in all nine
output columns. 864 rows added, none moved.

The classification is counted per temperature now, because it moves:

T	interior	pinned wet	pinned dry	shutdown
15 °C	198	24	18	48
25 °C	198	24	18	48
35 °C	194	43	3	48
40 °C	160	80	0	48

A single total would have let points move between branches and still add up. The
direction is physical: a hot leaf assimilates less and respires more, so it has less
to gain from water and its optimum presses against the wet bound instead of the
dry one. The 48 shut-down rows are temperature-independent, which is also right —
those are the ψ_soil = 6 MPa rows, where hydraulics forbid transpiration whatever the
temperature.

Two things the hot end does not reach, so the zeros are not read as coverage:
shade-death (reached by light, not by heat) and the compensation-point exit, which
needs about 45 °C at the defaults. test_rd_temperature_response covers the latter.

Still not addressed, and arguably the bigger question

R_d now rises with temperature, which was the ask. The calibration finding that
motivated it is not explained: leaf_calibration_test's ManyPeaksRange misfit is
reducible (a within-site random forest gets R² 0.45 on A where phylloptim gets
−0.27), and the carbon-shape defaults were refuted as the cause by a control.
R_d was the last untested carbon-side candidate. That study can free it now —
sabot-2022/R/ LEAF_OPTIONAL_PARS is where it would go — and that is the test of
whether any of this actually helps.

Verification

• C++ suite 496 checks, 0 failures; golden 1152 points, all bit-identical
• R suite 940 checks, 0 failures, 1 skip
• R CMD check: 1 NOTE, pre-existing (tibble not declared for a vignette)
• make -C tests/cpp bench builds (not in all, but CI builds it)
• CMake no-R build + ctest 2/2, and the installed-package consumer exit 0

[dfalster]

Follow-up: the recorded gradients needed their own tolerance (e714703)

The first CI run caught something worth writing down, on Linux only:

gradient_golden.tsv: worst relative difference 0.0013 (compared with tolerance)
── Failure: the recorded gradients have not moved ──
Expected interior-1layer R_d_25 A to equal `expected`.
  `actual`: -0.2484      `expected`: -0.2487

Not a regression — the tolerance was borrowed and the borrowing was wrong. This
file used golden_tolerance() (1e-03 for the argmax class), which is the
solved-output baseline's policy. But a value in gradient_golden.tsv is a finite
difference of the solve, so cross-platform it carries the solve's ~1e-9 floor
divided by the step — one amplification more than the solved outputs it is built
from. The step is relative, so the amplification is set by the differentiated
parameter's own magnitude, and the smallest parameter in the file sets the tolerance
for all of it. That used to be stem_b at 3.9; since this PR it is R_d_25 at
1.44, against vcmax_25 at 96 — a 67× smaller absolute step.

The arithmetic closes on the floor exactly: 1.3e-03 × 0.2487 × 2h = 9.3e-10 with
h = 1e-6 × 1.44. A step sweep agrees independently — at step = 1e-7 the R_d_25
column has already moved 1.3% while vcmax_25 has moved 0.2%, so it is the column
that goes noise-dominated first:

step	dA/dvcmax_25	dA/dstem_b	dA/dR_d_25
1e-04	0.020939751	2.8465363	−0.24873735
1e-05	0.020939195	2.8465389	−0.24873711
1e-06 (default)	0.020934892	2.8465358	−0.24873740
1e-07	0.020894746	2.8474104	−0.25197702

gradient_golden_tolerance() is 5e-03, 3.8× the one observation there is, and it
cannot hide a real change: this PR's own reallocation moved these cells by 16–250%,
two orders above it. The values stay recorded at the default step, so the file
still pins what a default caller gets — recording them at a coarser step would have
been more portable and would have stopped pinning the default.

⚠️ macOS never takes the tolerant branch, so the fix was verified by forcing it:
the CI disagreement passes at 5e-03, fails under the old 1e-03, and a real 5%
regression still fails.

Three stale figures found while reading the CI summary line

The guide says "if you need a magnitude, read the summary line, not the FAIL lines".
Doing that turned up three comments that nothing asserts, so all three had rotted:

• test_golden.cpp's own header carried profit 1.85e-06 / 5.87e-07 and the
  argmax class 5.53e-04 / 2.73e-04 — the pre-11a figures, two generations old,
  still showing a gcc-versus-clang split that 11a removed.
• The developer guide's profit row said 1.82e-07; CI now reports 2.14e-09,
  both compilers identical, over 1152 points with 7702 of 10368 values differing.
  ⚠️ Deliberately not attributed: a bigger grid can only move a maximum upward
  and the 25 °C block is byte-identical, so the improvement is a change on master
  since the figure was recorded (Bound the root vulnerability curves past the grid #77 and Solve for leaf temperature inside the collar FOC #84 are the candidates) — and confirming
  that needs a cross-platform run of master, which no macOS box can give you.
• The guide's gradient figure said 1.3e-4, which is the solved-output class, not
  this file's.

The sqrt(worst profit) ≈ worst argmax sanity check no longer holds numerically
(sqrt(2.14e-09) ≈ 4.6e-05 against 1.4e-04 observed) and the guide now says why that
is expected rather than alarming: the two column maxima fall at different operating
points, and the argmax column now has 864 hot rows the profit column's worst case is
not among. It was never a global identity.

[dfalster]

Good implementation, but we have put in too much detail ensuring and documenting back compatibility. back compatibility is ensured by using the same values as was previously used. I don't think we need extra infrastructure. We should treat R_d_25 as a first class trait with no special fallbacks.

Also too much text in the comments, particularly related to past history. We want to simplify for the future.

[dfalster]

Simplify

[dfalster]

Let's fix it properly with reorder

[dfalster]

remove. Always require respiration to be set

[dfalster]

why no default?

[dfalster]

Is this needed?

[dfalster]

is this needed?

[dfalster]

Done — all of it. c074647.

R_d_25 is a plain trait with no fallbacks. Default 1.44, which is
0.015 * 96 and — checked — the same double (0x1.70a3d70a3d70ap+0), so
back-compatibility comes from the value as you said. Deleted: the NaN sentinel,
rd_to_vcmax_ratio_ (C++ member and R binding; plant never used it), and
rd_tracks_vcmax_. An unset or negative R_d_25 now fails rather than deriving
anything. The special validator in leaf_traits() is gone too — .check_scalars()
covers it like every other trait.

Enumeration reordered properly. R_d_25 is the fourteenth trait, in set_traits'
own argument order, with leaf_specific_conductance_max and resistance after it.
par_kmax → 14, par_resistance → 15. All indices are named now
(par_vcmax_25 … par_resistance), so apply() no longer mixes theta[0] with
theta[par_R_d_25], and the R side is back to c(unlist(traits), kmax, resistance)
with no set difference to reconcile. The recorded gradients came out bit-identical
across the reorder — it is purely positional.

Golden grid: two temperatures, 25 and 40 °C. 1152 rows → 576. The 25 °C block is
still byte-identical, so the earlier regeneration is still checkable as an addition.

Both validate harnesses removed, and their numbers kept in PLAN 41. rd_shape_grid.cpp
had to go anyway — it existed to drive the escape hatch.

Comments cut throughout: −1965 lines against +372. The history is out of
leaf_model.hpp, gradient.hpp, the YAML, the workflow, the tests, CLAUDE.md and the
vignettes; NEWS is five #41 sections down to three. Two stale comments turned up while
cutting: the Doxygen block above update_temperature_dependent_params still said R_d
"falls above the thermal optimum … recorded rather than changed", and the cache-key
comment still cited rd_to_vcmax_ratio_ as its example. Both fixed.

---

One consequence to flag, measured rather than argued

⚠️ Respiration no longer follows vcmax_25. That is the point of a trait, and it
is the one place "same values ⇒ back compatible" has an exception — a run that
varies vcmax_25 used to get proportionally more or less respiration. plant's TF24
SCM offspring production, rebuilt both sides:

	phylloptim master	this branch
default (vcmax_25 = 96)	81.857087216483691	81.857087216483691
vcmax_25 = 60	0.40692	0.34336 (−15.6%)
vcmax_25 = 150	73.703	92.026 (+24.9%)

The default arm is bit-identical to all 17 digits, so your premise holds exactly where
it matters: plant's vcmax_25 is 96 and it pins leaf_temp at 25 °C. The direction of
the other two is right — a high-Vcmax strategy no longer pays proportionally more
respiration. If a study wants that coupling, TF24_hyperpar is where plant already
derives its other parameters from traits, and that seems the right home for it. Flagged
rather than acted on.

And a correction to my own earlier justification

My argument for the escape hatch was that "no published plant result under the old
shape could be regenerated". The revalidation I then ran undercut it: plant pins
leaf_temp at a constant 25 °C with PM off, and the two shapes agree exactly there
— so those results are regenerable without any hatch. Removing it is right, and my
reason for adding it was wrong.

Verification

486 C++ checks, golden 576 points bit-identical, 933 R checks, ctest 2/2, the
installed-package consumer exit 0, R CMD check 1 pre-existing NOTE (tibble).
⚠️ make bench caught one thing make all could not: bench_gradient.cpp was
passing the old NaN sentinel and now aborts on the validation — fixed.

[dfalster]

No — you caught a real miss. 22d412b fixes it.

Three edits from the last round never landed. They ran from a drifted working
directory, so the file was written but somewhere else, and I reported success from the
tool's "ok" rather than from the file. The developer guide still described a
1152-point grid at four temperatures, and both vignettes still described R_d_25
as an appended parameter with a NaN sentinel — including a whole ⚠️ paragraph telling
readers to resolve a sentinel that no longer exists.

So this time I swept for every count instead of the ones I remembered:

	was	now
golden grid	1152 points / 10368 cells / "4 leaf temperatures"	576 / 5184 / 25 and 40 °C
"thirteen traits"	in nine places	fourteen
gradient_par_names() roxygen	"a character vector of fifteen names"	sixteen
bench_gradient's total line	"all 13 traits"	"13 of the 14" — which is what it measures
the YAML's temperature-response block	"these thirteen"	fifteen bound fields

The cross-platform figures in the table were not stale — CI reports the same
profit 2.14e-09 / argmax 1.4e-04 on the smaller grid — but the prose around them
quoted the old cell counts, so those came from CI's current summary line.

Two judgement calls in there:

• bench_gradient stays at 13 traits. R_d_25 needs no spline rebuild, so it lands
  in the cheap bucket the split already has twelve examples of — adding it would change
  the recorded totals tools/bench_history.sh compares without measuring anything new.
  Labelled accurately rather than grown.
• test-surface.R's "13 traits and 4 tolerances onto 17 positional slots" is
  correct and left alone: that is about the constructor, which still takes 13.

And the history you flagged is out of the golden-file section: the "it was 288
points until #41" paragraph, the 1.82e-07 retraction and the sqrt-story archaeology
are gone. What is left is the mechanism (flat maximum, sqrt(dp/k), why the two
classes sit five orders apart), the current figures, and the one instruction that
matters — read the summary line, not this table, because nothing asserts it.

486 C++ checks, golden bit-identical, 933 R checks, R CMD check 1 pre-existing NOTE,
bench builds, all six CI jobs green.