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

.claude/CLAUDE.md

@@ -72,7 +72,8 @@ tests/cpp/                     plain-C++ suite, no R, no framework
 tests/cpp/root_network.hpp     the suite's root-architecture fixture: the two
                                ex-Leaf-default beta_R_* constants, in ONE place
                                because the golden file's bit-exactness depends on them
-tests/cpp/golden/              bit-exact regression baseline, 288 operating points
+tests/cpp/golden/              bit-exact regression baseline, 576 operating points
+                               -- one 288-point state grid at 25 and 40 C
 tests/cpp/bench_solve.cpp      timing harness for the collar solve (hazard 5)
 tests/cpp/bench_gradient.cpp   timing harness for a TRAIT GRADIENT: the IFT
                                composite against differencing the solve, with
@@ -87,7 +88,10 @@ tests/testthat/gradient_golden.tsv
                                on that platform only, like tests/cpp/golden/ --
                                these are derivatives of argmax-evaluated outputs,
                                so they inherit its sqrt-amplified class and
-                               disagree cross-platform by up to 1.3e-4
+                               disagree cross-platform by up to 1.3e-3 -- ten
+                               times the solved outputs, because a finite
+                               difference divides the solver floor by the step
+                               and `R_d_25`'s step is the smallest here
 tests/validate/                R scripts comparing against plant (needs R)
 CMakeLists.txt                 the no-R build: C++ and Python consumers, and the
                                thing that makes "does not need R" runnable
@@ -276,12 +280,20 @@ while looking like history. Two rows were produced that way before it was notice
 
 ## The golden file is the safety net — treat it that way
 
-`tests/cpp/golden/operating_points.tsv` records 288 operating points and is
-compared **bit-exactly** (`%.17g` round-trips a double). It is what makes a large
-refactor of this code checkable rather than hopeful, and it earned that role
-repeatedly: it proved three "surely dead" `set_physiology` arguments really were
-dead, confined the shutdown fix to exactly 48 rows × 5 fields, and showed the
-`area_leaf` change was 2 ULP.
+`tests/cpp/golden/operating_points.tsv` records 576 operating points — one
+288-point state grid at 25 and 40 °C — and is compared **bit-exactly** (`%.17g`
+round-trips a double). It is what makes a large refactor of this code checkable
+rather than hopeful, and it earned that role repeatedly: it proved three "surely
+dead" `set_physiology` arguments really were dead, confined the shutdown fix to
+exactly 48 rows × 5 fields, and showed the `area_leaf` change was 2 ULP.
+
+⚠️ **A grid at one leaf temperature is blind to every temperature response in the
+model**, because every reference value is *defined* at 25 °C and a change to any
+response curve is inert there **by construction**. Hence the second temperature, and
+hence the classification being printed per temperature: points move between branches
+as the leaf warms (dry-pinned 18 → 0, wet-pinned 24 → 80). Temperature is the
+OUTERMOST loop, which is what lets a regeneration that adds one be checked as an
+addition — the 25 °C block must come out byte-identical.
 
 **Only run `make golden` deliberately.** Running it after an accidental change
 rubber-stamps the change. If a diff is intended, regenerate and say so in the
@@ -336,33 +348,26 @@ magnitude apart:
 
 | field | gcc | clang | why |
 |---|---|---|---|
-| `profit` | **1.82e-07** | **1.82e-07** | it is the maximum itself — well-conditioned |
+| `profit` | **2.14e-09** | **2.14e-09** | it is the maximum itself — well-conditioned |
 | the other eight | **1.4e-04** | **1.4e-04** | evaluated at the **argmax** — sqrt-amplified |
 
-**These figures changed when PLAN 11a replaced the collar solver, and how they
-changed is informative.** They were `profit` 1.85e-06 / 5.87e-07 and the other
-eight 5.53e-04 / 2.73e-04. Two things to take from the move:
-
-- **The gcc-versus-clang split was itself a golden-section artefact.** The two
-  compilers now report *identical* figures, where they used to differ by 2–3×. What
-  differed between them was which way a golden-section comparison fell; what is
-  left is libm's `exp`/`pow`, which is a property of the platform and not of the
-  compiler. So do not expect a compiler-dependent column here any more — and if one
-  reappears, something has reintroduced a discrete decision into the solve.
-- **The sqrt story now fits better than it did.** `sqrt(1.82e-07)` ≈ 4.3e-04
-  against the 1.4e-04 observed, where before it was `sqrt(1.85e-06)` ≈ 1.4e-03
-  against 5.5e-04. Same order in both cases, but the residual factor shrank, which
-  is what you would expect once the argmax stopped carrying an extra `GSS_tol_abs`
-  of arbitrary displacement on top of the flat-maximum amplification.
-
-The maximum is *flat*: curvature measured directly at the two worst points gives
-k ≈ 1.0 and 0.9 in `profit ≈ p* − k(psi_stem−x*)²`. For a flat maximum an error
-`dp` in the profit **value** displaces its **location** by `sqrt(dp/k)`, which is
-why the well-conditioned column sits three orders below the other. Checked
-pointwise — at those points the residuals imply k = 1.01 and 1.70 against the 1.0
-and 0.9 measured from the curvature. It is *not* a global identity: the two column
-maxima fall at different operating points, so `sqrt(worst profit)` is not meant to
-reproduce `worst argmax`.
+Read off CI's summary line: 3757 of 5184 values differ, against tolerances of 1e-05
+and 5e-03. ⚠️ **This table is a CI reading and nothing asserts it, so it goes stale
+silently** — it has been wrong three times. Read the summary line, never the FAIL
+lines, which are truncated at 20 and biased toward whichever rows come first.
+
+gcc and clang report *identical* figures. They used to differ by 2–3×, and what
+differed was which way a golden-section comparison fell; since PLAN 11a removed that
+search, what is left is libm's `exp`/`pow` — a property of the platform, not the
+compiler. If a compiler-dependent column reappears, something has reintroduced a
+discrete decision into the solve.
+
+**Why the two classes are five orders apart:** the maximum is *flat*, with curvature
+k ≈ 1.0 measured directly at the worst points in `profit ≈ p* − k(psi_stem−x*)²`, so
+an error `dp` in the profit **value** displaces its **location** by `sqrt(dp/k)`.
+That is a mechanism, not an identity — `sqrt(2.14e-09)` ≈ 4.6e-05 against the 1.4e-04
+observed, because the two column maxima fall at *different* operating points. Do not
+expect `sqrt(worst profit)` to predict `worst argmax`.
 
 Two things follow that matter beyond this file:
 

.github/workflows/cpp-tests.yml

@@ -188,9 +188,13 @@ jobs:
             d.PPFD = 900; d.psi_soil = {2.0}; d.soil_depth = soil_depth;
             d.atm_vpd = 2.0; d.ca = 40.0; d.leaf_temp = 25.0;
             d.atm_o2_kpa = 21.0; d.atm_kpa = 101.3;
+            // ⚠️ ALL n_pars ELEMENTS, WRITTEN OUT: a short initialiser compiles
+            // and zero-fills the rest, so adding a parameter would silently pass
+            // 0.0 for it here. The traits in set_traits' order, then the two
+            // non-traits.
             double theta[phylloptim::gradient::n_pars] = {
                 96.0, 2.680147, 3.898245, 5.870283, 2.680147, 3.898245,
-                5.870283, 1.5, 157.44, 0.30, 0.7, 0.99, 7.5,
+                5.870283, 1.5, 157.44, 0.30, 0.7, 0.99, 7.5, 1.44,
                 1.0 * 0.000157 / 5.0, 0.0};
             const int pars[1] = {0};  // vcmax_25
             const std::vector<phylloptim::gradient::Drivers> obs{d, d};
@@ -207,7 +211,7 @@ jobs:
             }
             // The R layer's arbitrated reference for this operating point, loosely
             // compared for the same reason profit is below.
-            if (std::fabs(g[0].grad[0] - 0.0172086) > 1e-5) {
+            if (std::fabs(g[0].grad[0] - 0.0209349) > 1e-5) {
               return 1;
             }
             l.set_physiology(roots, 900, {2.0}, soil_depth,

NEWS.md

@@ -1,5 +1,167 @@
 # phylloptim 0.2.1
 
+## `R_d_25` is a trait, and respiration rises with temperature (#41)
+
+Dark respiration was `rd_to_vcmax_ratio_ * vcmax_(T)`: a fixed fraction of Vcmax,
+unreachable from R, and — because it inherited Vcmax's **peaked** Arrhenius — it
+**fell** above the thermal optimum, where real dark respiration rises. No value of
+the fraction repairs a function of the wrong shape.
+
+Now `R_d_25` is a `leaf_traits()` trait like any other, and the response is
+Tjoelker's declining Q10:
+
+```
+Q10(T) = rd_q10_intercept_ - rd_q10_slope_ * (T + 25) / 2      3.09, 0.0430
+R_d(T) = R_d_25 * Q10(T)^((T - 25) / 10)
+```
+
+The Q10 is evaluated at the mean of T and the reference, which is the form the
+land-surface literature implements. Set the slope to zero and the intercept is a
+constant Q10, for anyone who wants the conventional form; a constant Q10 of 2 was
+tried first and measured too aggressive (R_d 2.8× its 25 °C value over 15 K, −73% on
+assimilation at 40 °C, no operating point at all by 45 °C).
+
+**There is no fallback and no sentinel.** `R_d_25` defaults to **1.44** — which is
+`0.015 * 96`, the value the old derivation gave at the default `vcmax_25`, and the
+same double — so a caller who does not touch it gets exactly what they got before.
+An unset or negative value fails rather than deriving something.
+
+⚠️ **RESPIRATION NO LONGER FOLLOWS `vcmax_25`.** That is the point of making it a
+trait, and it is the one place "same defaults ⇒ same results" has an exception:
+a run that *varies* `vcmax_25` used to get proportionally more or less respiration
+and now does not. Measured on plant's TF24 SCM, offspring production:
+
+| | phylloptim master | this branch |
+|---|---|---|
+| default (`vcmax_25` = 96) | 81.857087216483691 | 81.857087216483691 |
+| `vcmax_25` = 60 | 0.40692 | 0.34336 |
+| `vcmax_25` = 150 | 73.703 | 92.026 |
+
+The default arm is bit-identical to all 17 digits. If a study wants respiration to
+track Vcmax, that coupling belongs in the caller's parameterisation — plant's
+`TF24_hyperpar` is where every other derived parameter is already computed.
+
+`R_d_25` is also a `pars` entry for `leaf_gradient()` and `leaf_gradient_batch()`:
+`dA/dR_d_25 = -0.24874` against a central difference of the solve at −0.24868. ⚠️ Not
+≈ −1, because the optimiser moves the operating point in response. `dY/dvcmax_25` is
+correspondingly a **partial** at fixed respiration; a fit that moves both adds the
+two columns.
+
+**Blast radius**, and it is zero at 25 °C by construction — the trait *is* the value
+there:
+
+| 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 |
+
+Two further consequences:
+
+- **A leaf too hot to gain carbon now shuts down instead of throwing.** A higher R_d
+  can drive net assimilation negative across the whole `[gamma*, ca]` bracket, and
+  then `psi_stem_to_ci` has no supply-equals-demand root. The model already knew that
+  case — `ci_at_compensation_point_` — but only on the energy-balance path; it now
+  applies on the default path too, and the temperature at which the model stops
+  having an operating point is pinned as a measurement.
+- **On plant's Penman-Monteith path the change is large**, because leaf temperature
+  there reaches **62 °C, up to +22 K above air**: offspring production −47.2%
+  (2.2035 → 1.1638). ⚠️ plant's *default* configuration is unaffected — it sets
+  `leaf_temp` as a constant 25 °C driver with PM off — so the exposure is entirely in
+  configurations that leave 25 °C.
+
+`R_d_25` is the fourteenth entry of the gradient enumeration, i.e. in `set_traits()`'
+own argument order, with `leaf_specific_conductance_max` and `resistance` after it.
+Every index in `gradient.hpp` now has a name (`par_vcmax_25`, `par_R_d_25`, ...) so
+nothing indexes `theta` with a bare integer.
+
+**API removal:** `rd_to_vcmax_ratio_` is gone, from both C++ and the R bindings. It
+existed only as the fraction, and there is no fraction now. plant never used it.
+
+## The golden grid gains a temperature axis: 288 points become 576 (#41)
+
+⚠️ **A grid at one leaf temperature is blind to every temperature response in the
+model**, and this one was, at 25 °C. The reference values of Vcmax, Jmax and R_d are
+*defined* there, so a change to any response curve is inert **by construction** — the
+respiration change above moved results at every temperature except 25 °C and this
+file did not move a bit.
+
+The grid is now the same 288 operating points at **25 and 40 °C**, with a `leaf_temp`
+column. 40 °C is where the response bites hardest, and the extra pinned rows at the
+hot end are a feature.
+
+**Regenerated deliberately, and checkable as an addition**: temperature is the
+outermost loop, so the 25 °C block is **byte-identical** to the previous file in all
+nine output columns. 288 rows added, none moved.
+
+The classification is counted per temperature now, because it moves:
+
+| T | interior | pinned wet | pinned dry | shutdown |
+|---|---|---|---|---|
+| 25 °C | 198 | 24 | 18 | 48 |
+| 40 °C | 160 | 80 | 0 | 48 |
+
+A single total would 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 do not move, because there it is hydraulics rather than heat
+that forbids transpiration.
+
+`tests/cpp/bench_solve.cpp` 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`).
+
+## The recorded gradients have their own tolerance
+
+`gradient_golden.tsv` used the solved-output baseline's tolerance, and that was
+wrong. A value there is a **finite difference**, so cross-platform it carries the
+solve's ~1e-9 floor **divided by the step** — one amplification more than the outputs
+it is built from. The step is relative, so the smallest-magnitude parameter sets the
+tolerance for the whole file: `R_d_25` is 1.44 where `vcmax_25` is 96, a 67× smaller
+absolute step, and it disagrees 1.3e-03 on Linux against 1.3e-04 for the other
+columns. `gradient_golden_tolerance()` is 5e-03, with the arithmetic.
+
+## The temperature cache now keys on everything it reads, so R_d is genuinely settable (#41)
+
+Binding the temperature-response parameters to R made them *writable* but not
+*effective*. `set_physiology()`'s temperature block was cached on
+`(leaf_temp_, atm_o2_kpa_)` alone, justified by "same inputs → bit-identical
+outputs, so reusing is exact" — a statement that was true while those parameters
+were unreachable C++ members and **became false the moment they were bound**. The
+block's outputs depend on fourteen further inputs the key never mentioned.
+
+The failure was silent. Setting `rd_to_vcmax_ratio_ <- 0.03` on a solved leaf and
+re-supplying the same drivers left `R_d` at **1.44** where a freshly built leaf gave
+**2.88**, and assimilation unchanged to every digit. Which is the outcome #41 cared
+about: a calibration that cannot move respiration absorbs the mismatch into whatever
+it *can* move.
+
+The cache key now covers every scalar `update_temperature_dependent_params()` reads.
+Two consequences beyond the reported bug:
+
+- **It also covers `vcmax_25` and `jmax_25`**, which closes the third and least
+  visible part of hazard 10 — a bare `l$vcmax_25 <- x` no longer leaves
+  `vcmax_`/`jmax_`/`R_d_` describing the old value. ⚠️ `set_traits()` is still the
+  correct way to change a trait; the vulnerability splines and the solved operating
+  point need clearing too, and a cache key cannot do that.
+- ⚠️ **The existing test passed while the feature was broken**, because it pushed
+  each change through by calling `update_temperature_dependent_params()` directly —
+  a route no caller has. `test_temperature_params_invalidate_cache` uses the real
+  one, and fails four ways without this fix.
+
+**No behaviour change at the defaults**: the golden file is bit-identical, and the
+solve is 2.99 against 3.00 µs/solve interleaved ×3, inside the ±0.01 within-process
+noise.
+
+⚠️ **What this does NOT fix**, recorded because the numbers invite it: `R_d` still
+inherits Vcmax's *peaked* Arrhenius and therefore **falls** above the thermal
+optimum, where real dark respiration rises. No value of the ratio repairs a function
+of the wrong shape — see the note in `update_temperature_dependent_params()`. And
+`rd_to_vcmax_ratio` is still not a `leaf_traits()` member, so it cannot yet be
+fitted or differentiated; it is set as a field.
+
 ## An out-of-domain transport lookup says which spline, and which caller
 
 The stem curve is the only interpolator here built with extrapolation disabled, so