diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index 6013b17..187b819 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -252,6 +252,15 @@ byte-for-byte. If you add a header, you need do nothing. If you want a literal Doxygen command, write a `///` comment and the filter will leave it alone. Publishing is off until someone sets the repo variable `PUBLISH_DOCS=true`. +⚠️ **One indented display per comment block.** An indented run becomes +`\verbatim`, and Doxygen 1.9 — which is the version CI installs, where local +Homebrew is 1.17 — drops the second one's OPEN and then reports `unexpected +command endverbatim` at a line in the *filtered* stream that lands in unrelated +code, 120 lines away in the case that found it. **A local `doxygen` run renders +it in silence**, so this is only ever visible in CI and only cryptically. Put the +equations in one display rather than one on each side of a paragraph; `docs.yml` +now asserts the shape and names the block, so it fails legibly. + ## Cost: the one thing the golden file cannot see `tests/cpp/bench_solve.cpp` and `bench_gradient.cpp` cover the C++ side. @@ -815,6 +824,31 @@ the per-cause split and the tolerance bands go in the first PR comment — see removed the same composite wins 4.4× on the eleven traits that touch no spline. +11. **The single-layer optimisers search a LOCAL maximum, and the profit is + neither unimodal nor interior.** `brent_fmin` steps in from the bounds, so it + can return neither an endpoint nor the global maximum of a multi-modal + objective. Both happen: at a leaf hot enough that net assimilation is negative + across the whole supply stream, the ProfitMax profit is highest at **full + closure** and carries a local maximum out in the interior. Measured at Tair + 50 °C with the thermal cost on — profit −1.5314 at `psi_soil`, −1.5510 at 1.19, + −1.5459 at 1.88 — and the solver returned **1.643**, reporting an open stoma + where the objective says shut. + + `optimise_psi_stem_ProfitMax` now evaluates the objective on the scan + `prepare_profitmax` already runs, takes the grid argmax, and refines with Brent + only when it is interior. No extra model evaluations: `A` and `Tleaf` are + stored per grid point and `HC`/`TC` are analytic. + + ⚠️ **`optimise_psi_stem_TF` and `optimise_psi_stem_Sperry` still have it**, and + are documented rather than fixed because neither has a scan to reuse. The + collar solve is unaffected — `maximise_profit_over_collar` handles a pinned + optimum explicitly, which is why 42 of 240 feasible golden rows are pinned and + correct. + + The general form: **a bracketing optimiser answers "where is the interior + maximum", and that is not the same question as "where is the maximum".** If an + objective can be maximised at a constraint, the search has to be told. + ## Validating against plant ⚠️ **There is no longer a harness that compares against plant, and there cannot be diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index eeac8b8..e81bb74 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -68,6 +68,53 @@ jobs: done exit $status + # ⚠️ Doxygen 1.9 -- the version this runner installs -- handles the FIRST + # `\verbatim` in a comment block and drops the second one's OPEN, then + # reports `unexpected command endverbatim` at its close. The line it names + # is in the FILTERED stream and lands in unrelated code: in the case that + # found this, 120 lines past the comment responsible. Doxygen 1.17 renders + # the same input in silence, so a local run says nothing and the only + # oracle is this job. Assert the shape here, where the message can name the + # block, rather than leaving the next person to bisect the render. + # + # ⚠️ It is the `///` blocks that break, NOT the `/*! \file */` one, and the + # difference is load-bearing rather than untested: closed_form.hpp's file + # block has carried two runs since long before this check, on a green + # master. Counting it too fails the build on code Doxygen renders happily. + - name: Check no comment block opens two verbatim runs + run: | + set -eu + status=0 + for f in $(find inst/include -name '*.hpp' | sort); do + awk -f tools/doxygen_filter.awk "$f" \ + | awk -v f="$f" ' + function note(t) { + if (first == "") { + sub(/^[ \t]*(\/\/\/|\*)?[ \t]*/, "", t) + first = substr(t, 1, 60) + } + } + function flush() { + if (opens > 1) { + printf "::error file=%s::comment block [%s] opens %d verbatim runs; Doxygen 1.9 drops all but the first. One indented display per block -- see tools/doxygen_filter.awk\n", f, first, opens + bad = 1 + } + opens = 0; first = "" + } + /^[ \t]*\/\*![ \t]*\\file/ { infile = 1; flush(); next } + infile { + if ($0 ~ /^[ \t]*\*\/[ \t]*$/) infile = 0 + next + } + /^[ \t]*\/\// { + if ($0 ~ /\\verbatim[ \t]*$/) opens++ + note($0); next + } + { flush() } + END { flush(); exit (bad ? 1 : 0) }' || status=1 + done + exit $status + # Overriding on stdin is Doxygen's documented way to change one setting # without a second Doxyfile. Warnings stay non-fatal for a local `doxygen` # run, where they are informative rather than blocking. diff --git a/COMPARISON.md b/COMPARISON.md index a43d855..1366d24 100644 --- a/COMPARISON.md +++ b/COMPARISON.md @@ -28,10 +28,10 @@ source read at 1.0.6.1). All three are **pure R with no compiled code**. | Soil water potential | **yes, multi-layer** (separable — PLAN 7b) | prescribed scalar (Tuzet only) | no | no | | Root resistance | **yes, per layer** | explicitly not implemented | no | no | | Profit / gain-risk optimisation | **yes** | Cowan-Farquhar with fixed λ | no | no | -| Thermal acclimation or damage | yes (in plant's TF24t) | no | no | no | +| Thermal acclimation or damage | yes — an instantaneous PSII cost here (Sicangco et al. 2026, default off), and a lasting damage ratchet in plant's TF24t | no | no | no | | Fits to measured data | **no** | yes (A-Ci, Ball-Berry) | no | yes (that's the point) | | Swappable *empirical* gs schemes | Medlyn present, not dispatched | **yes — 4, plus Tuzet** | n/a | n/a | -| Swappable *hydraulic optimality* schemes | **planned — TF24, Sperry, Prentice14** | no (none are hydraulic) | no | no | +| Swappable *hydraulic optimality* schemes | **TF24 and Sperry ProfitMax runnable at identical drivers; Prentice14 planned** | no (none are hydraulic) | no | no | | Exact derivatives | **yes — forward-mode AD (XAD)** | no | no | no | | Language | C++ header-only | R | R | R | | Cost per solve | ~4 µs | not measured; `mapply` over scalars | one `uniroot` per leaf | vectorised over time series | diff --git a/NEWS.md b/NEWS.md index f629010..c90520d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -570,6 +570,126 @@ of the wrong shape — see the note in `update_temperature_dependent_params()`. `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. +## Fick's law uses the LEAF-to-air vapour deficit, not the air's (#7) + +Transpiration and stomatal conductance divided by `atm_vpd_` however hot the leaf +got. Diffusion out of a stoma is driven by the deficit *at the leaf*, and on the +energy-balance path the leaf runs above air temperature, so that deficit is the +larger one. Measured at the package defaults it is **3.4× / 3.2× / 4.0× the air's +at Tair 25 / 35 / 45 °C** — the factor by which `gs` and `A` were overstated +([#7](https://github.com/traitecoevo/phylloptim/issues/7), PLAN 13.1). + +⚠️ **Off the energy balance nothing moves, and that is a property of the grid +rather than evidence the change is inert.** At air temperature +`vpd_leaf_ == atm_vpd_` exactly, so the prescribed-temperature path is untouched +and `tests/cpp/golden/` is **bit-identical across all five commits** — because +the grid runs with `use_energy_balance_` off. A bit-identical golden run says +nothing here. + +⚠️ **The decoupling signature moves from conductance to transpiration.** Over the +acceptance window `E` rises ×1.034 against `D_leaf` ×1.089, so `gs` falls ×0.950. +A test written against the old arithmetic would read the same physics off the +wrong variable. + +`constants::vpd_leaf_min = 0.01` kPa floors it. A leaf cooled below the dew point +has a *negative* deficit, and dividing a positive transpiration by it reports a +negative conductance; the floor keeps such points finite and unattractive instead +of sign-inverted. The driver defaults are 1.5–2 kPa, so it binds only where a +one-way diffusion equation has already stopped describing the leaf. + +## Sperry's ProfitMax, on the same footing as TF24 + +`$optimise_psi_stem_ProfitMax()` maximises Sperry (2017)'s profit with **both** +terms normalised, and reports the λ that makes the existing entry point agree. +`$profitmax_curve()` returns the whole cost, gain and profit curve in one call, +and `$optimise_psi_stem_TF()` runs the TF24 cost through the same solver — so the +two formulations can be compared at identical drivers. An optional instantaneous +thermal cost is included, **default off**. + +⚠️ **These write `$carbon_gain_`, `$hydraulic_cost_norm_` and `$thermal_cost_`, +which are unitless and are NOT `$hydraulic_cost_`.** Reading the normalised cost +as the TF24 one is a units error the names are chosen to prevent. + +The normalised hydraulic cost is invariant to `kmax` by construction, and it is +asserted rather than assumed: worst difference **1.1e-16** across a 3× change. + +Motivated by `leaf_calibration_test/sicangco-2026`, a replication of Sicangco et +al. (2026), which needed ProfitMax measured against TF24 rather than described. + +## ⚠️ Four fixes on the single-layer optimisers, and the fourth is a solver bug + +These paths are reachable only from `optimise_psi_stem_*`, not from +`find_root_collar_psi()` — so **`plant` is unaffected, provably rather than +probably**: `prepare_collar_solve` has its own `assim_max_ < 0` exit and returns +false before any candidate potential is evaluated. + +1. **`optimise_psi_stem_Sperry` searched a NaN objective silently.** `lambda_` is + an input with no default, cleared by `setup_clean_leaf` and never set by + `set_physiology`, so a caller who drove the leaf and called this got a + plausible potential (2.551266 MPa at the defaults) beside `profit_ = NaN`. It + now refuses. +2. **Hazard 8, live on this path.** `opt_root_psi_`, `E_up_` and + `soil_consumption_` survived from an earlier collar solve, so the object + reported `E = 9.216e-5` beside `E_up = 2.626e-5`. +3. **`set_leaf_states_rates_from_psi_stem` zeroed transpiration wherever + `assim_max_ < 0`.** Transpiration there is the hydraulic supply and does not + depend on photosynthesis, so the branch made two leaves at the same operating + point disagree about whether water was moving — purely because one had + respiration switched on. The shut-down *state* it reached for is unchanged: + the `ci` root-find has no root in `[gamma*, ca]` and takes its + compensation-point fallback, which is what the branch set by hand. That + fallback did not exist when the branch was written. +4. **`optimise_psi_stem_ProfitMax` now scans a grid before refining.** + `brent_fmin` steps in from the bounds, so it returns neither an endpoint nor + the global maximum of a multi-modal objective — and at Tair 50 °C with the + thermal cost on it returned **1.643 MPa where the objective is maximised at + full closure** (−1.5314 at `psi_soil` against −1.5459 at the interior local + max), reporting an open stoma where the model says shut. The scan reuses + `prepare_profitmax`'s own grid, so this costs no extra model evaluations. + Found by comparing against Sicangco et al.'s Figure 4: their model closed at + 48–56 °C and ours did not, and the whole difference was the search. + +⚠️ **`optimise_psi_stem_TF` and `optimise_psi_stem_Sperry` still have fault 4**, +and are documented rather than fixed because neither has a scan to reuse. The +general form is now hazard 11 in the developer guide: *a bracketing optimiser +answers "where is the interior maximum", which is not the same question as "where +is the maximum".* + +## The `dgc_dT` term in the energy-balance derivative was re-derived, not sign-flipped + +`dprofit_energy_balance_term` carried a named `const double dgc_dT = 0.0` and a +note promising PLAN 13.1 would make it "a one-line change instead of a +re-derivation". It would not have been: the damping factor it multiplied puts the +new term under `A_T` where the derivation puts it under `A_prime`. The two agree +**only at `dgc_dT = 0`**, which is why nothing caught it. The derivation is +written out at the function, and at `dgc_dT = 0` it reduces to exactly the old +expression — so the prescribed-VPD behaviour is unchanged, which is the other +half of why the golden file does not move. + +## A comment block may open only one `\verbatim` run, and only CI can see it + +Doxygen 1.9.8 — what the runner installs — handles the first `\verbatim` in a +`///` block and drops the second one's OPEN, then reports `unexpected command +endverbatim` at a line in the *filtered* stream that lands in unrelated code: 120 +lines past the responsible comment, in the case that found it. Doxygen 1.17 +renders the same input in silence, so a local run is not an oracle. `docs.yml` +now counts opens per block and names the offending block. + +Established with a probe header of eight isolated constructs in one render: one +run, banner rules, `|` in prose, `|` inside a run, and emoji are all clean; only +two-runs-in-one-block errors. ⚠️ **Two earlier explanations were asserted before +being measured and are both wrong** — a setext heading swallowing the block, and +a bare `\|` opening a Markdown table. A `\|`-escaping change to +`tools/doxygen_filter.awk` had been committed on the second and written into the +filter's header as fact; **that commit is reverted**, and the filter's header now +records the measurement instead. + +⚠️ The `/*! \file */` block is exempt from the count, measured rather than +assumed: `closed_form.hpp`'s file block has carried two runs across a long green +master, so counting it would fail the build on code Doxygen renders happily. + +C++ suite **486 → 534 checks**, 0 failures. + ## A trait gradient at a collar potential the caller supplies `leaf_gradient()` and `leaf_gradient_batch()` take `psi`, and evaluate there diff --git a/PLAN.md b/PLAN.md index 1bcaabe..0fd0252 100644 --- a/PLAN.md +++ b/PLAN.md @@ -29,7 +29,7 @@ the user-visible history. | [#74](https://github.com/traitecoevo/phylloptim/issues/74) | 11f, 11g | The `stem_b` shortcut is undone by a rebuild once per observation | 11f's 24.5× is **2.4×** through `leaf_gradient_batch()`. Unmotivated until something frees a vulnerability curve; filed so the wrong figure is not quoted meanwhile | | [#6](https://github.com/traitecoevo/phylloptim/issues/6) | 12 | Real-data calibration, then inversion | **The one with the most downstream** — it specifies what is left of #4, and doing it found #38, #40, #41 and #52. ⚠️ The *synthetic* vignette is done ([#53](https://github.com/traitecoevo/phylloptim/pull/53)) and **AD lost**: read 12a before repeating the comparison | | [#52](https://github.com/traitecoevo/phylloptim/issues/52) | 12 | `leaf_gradient()` rebuilds a `Leaf` per call, with no way to pass one in | From #53. A third of a per-observation gradient's cost, and why item 12's "reuse one object" advice cannot be followed as written | -| [#7](https://github.com/traitecoevo/phylloptim/issues/7) | 13 | Energy balance, in priority order | Leaf-to-air VPD is the cheap win and is **not wired in**; free convection is not worth it | +| [#7](https://github.com/traitecoevo/phylloptim/issues/7) | 13 | Energy balance, in priority order | ⚠️ **Item 1, leaf-to-air VPD, is DONE.** What is left is temperature-dependent longwave (#28); free convection is not worth it | | [#28](https://github.com/traitecoevo/phylloptim/issues/28) | 13 | Temperature-dependent outgoing longwave in the Penman-Monteith Rn | from plant #581 / #567 review | | [#31](https://github.com/traitecoevo/phylloptim/issues/31) | 31 | `profit_psi_stem_TF` returns a plausible number below `psi_upstream` | Found writing the vignette. Profit is **discontinuous by 1.58** at the boundary | | [#34](https://github.com/traitecoevo/phylloptim/issues/34) | 6d | Delete plant's `Leaf` bindings | Unblocked (plant #591 merged, and #33 is done). The hazard-7 payoff, and the only stage that can break plant | @@ -86,12 +86,25 @@ The package carries three formulations and only one is usable: | formulation | state | |---|---| | TF24 hydraulic gain-risk | production; the whole solve is built on it | -| Sperry et al. (2017) cost | present but **hardwired to `psi_soil_[0]`**, and nothing routes to it | +| Sperry et al. (2017) ProfitMax | **runnable**: `optimise_psi_stem_ProfitMax()` on the single-potential path, with both terms normalised as the paper defines them | | Medlyn et al. (2011) USO | present but **bypasses the hydraulic solve altogether** | -So you cannot run the same drivers through two of them and compare, which is the -obvious thing to want. The goal is each as a first-class member, alongside Prentice -et al. (2014) least-cost and Cowan-Farquhar. +⚠️ **The Sperry row used to read "hardwired to `psi_soil_[0]`, and nothing routes to +it", and both halves were misleading.** Being evaluated at `psi_soil` is what Eqn 5 +of Sperry (2017) *specifies* — the cost is defined against `k(psi_soil)` — so that +was never the defect. The defect was that `profit_psi_stem_Sperry` uses the +unnormalised form `A - lambda*cost` with `lambda_` a prescribed input that nothing +sets, and `lambda_` is not a constant: the equivalent value is +`|A|max/(k_soil - kcrit)`, which moves with every driver. `optimise_psi_stem_ProfitMax` +computes it and reports it in `lambda_`, so the two forms can be checked against +each other rather than assumed equivalent. + +`optimise_psi_stem_TF()` is the same off-path solver for the TF24 cost, so **the same +drivers can now be run through both formulations on the same footing** — which this +section used to say was impossible. What is still missing is dispatch on the +PRODUCTION path (`find_root_collar_psi`, the collar first-order condition, the root +network); the goal is each as a first-class member there, alongside Prentice et al. +(2014) least-cost and Cowan-Farquhar. **What should be pluggable is λ, not the cost function.** All of these maximise a profit, so all satisfy `dA/dE = λ` and differ **only** in λ(state). Given λ, each diff --git a/R/RcppExports.R b/R/RcppExports.R index 5071c98..e51df28 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -77,6 +77,10 @@ Leaf__electron_transport <- function(obj_) { .Call('_phylloptim_Leaf__electron_transport', PACKAGE = 'phylloptim', obj_) } +Leaf__leaf_temp_from_E <- function(obj_, E) { + .Call('_phylloptim_Leaf__leaf_temp_from_E', PACKAGE = 'phylloptim', obj_, E) +} + Leaf__set_leaf_states_rates_from_psi_stem <- function(obj_, psi_stem, psi_upstream) { invisible(.Call('_phylloptim_Leaf__set_leaf_states_rates_from_psi_stem', PACKAGE = 'phylloptim', obj_, psi_stem, psi_upstream)) } @@ -137,6 +141,26 @@ Leaf__optimise_psi_stem_Sperry <- function(obj_) { invisible(.Call('_phylloptim_Leaf__optimise_psi_stem_Sperry', PACKAGE = 'phylloptim', obj_)) } +Leaf__thermal_cost_at <- function(obj_, leaf_temp) { + .Call('_phylloptim_Leaf__thermal_cost_at', PACKAGE = 'phylloptim', obj_, leaf_temp) +} + +Leaf__prepare_profitmax <- function(obj_) { + invisible(.Call('_phylloptim_Leaf__prepare_profitmax', PACKAGE = 'phylloptim', obj_)) +} + +Leaf__profit_psi_stem_ProfitMax <- function(obj_, psi_stem, psi_upstream) { + .Call('_phylloptim_Leaf__profit_psi_stem_ProfitMax', PACKAGE = 'phylloptim', obj_, psi_stem, psi_upstream) +} + +Leaf__profitmax_curve <- function(obj_, n) { + .Call('_phylloptim_Leaf__profitmax_curve', PACKAGE = 'phylloptim', obj_, n) +} + +Leaf__optimise_psi_stem_ProfitMax <- function(obj_) { + invisible(.Call('_phylloptim_Leaf__optimise_psi_stem_ProfitMax', PACKAGE = 'phylloptim', obj_)) +} + Leaf__optimise_psi_stem_TF <- function(obj_) { invisible(.Call('_phylloptim_Leaf__optimise_psi_stem_TF', PACKAGE = 'phylloptim', obj_)) } @@ -285,6 +309,30 @@ Leaf__lambda___set <- function(obj_, value) { invisible(.Call('_phylloptim_Leaf__lambda___set', PACKAGE = 'phylloptim', obj_, value)) } +Leaf__carbon_gain___get <- function(obj_) { + .Call('_phylloptim_Leaf__carbon_gain___get', PACKAGE = 'phylloptim', obj_) +} + +Leaf__carbon_gain___set <- function(obj_, value) { + invisible(.Call('_phylloptim_Leaf__carbon_gain___set', PACKAGE = 'phylloptim', obj_, value)) +} + +Leaf__hydraulic_cost_norm___get <- function(obj_) { + .Call('_phylloptim_Leaf__hydraulic_cost_norm___get', PACKAGE = 'phylloptim', obj_) +} + +Leaf__hydraulic_cost_norm___set <- function(obj_, value) { + invisible(.Call('_phylloptim_Leaf__hydraulic_cost_norm___set', PACKAGE = 'phylloptim', obj_, value)) +} + +Leaf__thermal_cost___get <- function(obj_) { + .Call('_phylloptim_Leaf__thermal_cost___get', PACKAGE = 'phylloptim', obj_) +} + +Leaf__thermal_cost___set <- function(obj_, value) { + invisible(.Call('_phylloptim_Leaf__thermal_cost___set', PACKAGE = 'phylloptim', obj_, value)) +} + Leaf__lambda_analytical___get <- function(obj_) { .Call('_phylloptim_Leaf__lambda_analytical___get', PACKAGE = 'phylloptim', obj_) } @@ -701,6 +749,46 @@ Leaf__use_energy_balance___set <- function(obj_, value) { invisible(.Call('_phylloptim_Leaf__use_energy_balance___set', PACKAGE = 'phylloptim', obj_, value)) } +Leaf__vpd_leaf___get <- function(obj_) { + .Call('_phylloptim_Leaf__vpd_leaf___get', PACKAGE = 'phylloptim', obj_) +} + +Leaf__vpd_leaf___set <- function(obj_, value) { + invisible(.Call('_phylloptim_Leaf__vpd_leaf___set', PACKAGE = 'phylloptim', obj_, value)) +} + +Leaf__use_thermal_cost___get <- function(obj_) { + .Call('_phylloptim_Leaf__use_thermal_cost___get', PACKAGE = 'phylloptim', obj_) +} + +Leaf__use_thermal_cost___set <- function(obj_, value) { + invisible(.Call('_phylloptim_Leaf__use_thermal_cost___set', PACKAGE = 'phylloptim', obj_, value)) +} + +Leaf__T50___get <- function(obj_) { + .Call('_phylloptim_Leaf__T50___get', PACKAGE = 'phylloptim', obj_) +} + +Leaf__T50___set <- function(obj_, value) { + invisible(.Call('_phylloptim_Leaf__T50___set', PACKAGE = 'phylloptim', obj_, value)) +} + +Leaf__Tcrit___get <- function(obj_) { + .Call('_phylloptim_Leaf__Tcrit___get', PACKAGE = 'phylloptim', obj_) +} + +Leaf__Tcrit___set <- function(obj_, value) { + invisible(.Call('_phylloptim_Leaf__Tcrit___set', PACKAGE = 'phylloptim', obj_, value)) +} + +Leaf__profitmax_scan_n___get <- function(obj_) { + .Call('_phylloptim_Leaf__profitmax_scan_n___get', PACKAGE = 'phylloptim', obj_) +} + +Leaf__profitmax_scan_n___set <- function(obj_, value) { + invisible(.Call('_phylloptim_Leaf__profitmax_scan_n___set', PACKAGE = 'phylloptim', obj_, value)) +} + Leaf__d___get <- function(obj_) { .Call('_phylloptim_Leaf__d___get', PACKAGE = 'phylloptim', obj_) } diff --git a/R/RcppR6.R b/R/RcppR6.R index 979afe9..c208457 100644 --- a/R/RcppR6.R +++ b/R/RcppR6.R @@ -1,6 +1,6 @@ ## Generated by RcppR6: do not edit by hand ## Version: 0.2.4 -## Hash: aef2740595352538d5bff284a5829060 +## Hash: d9d602ca256a4fa03646aa2a65245648 ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -53,6 +53,14 @@ NULL ##' `$opt_psi_stem_`, `$profit_`, `$assim_colimited_`, `$transpiration_`, ##' `$lambda`, `$g1_eff`. ##' +##' A second optimality formulation is available on the single-soil-layer +##' supply path: `$optimise_psi_stem_ProfitMax()` maximises Sperry (2017)'s +##' normalised profit, and `$profitmax_curve()` returns the whole cost, gain +##' and profit curve in one call. `$optimise_psi_stem_TF()` runs the TF24 cost +##' through the same solver, so the two can be compared at identical drivers. +##' These write `$carbon_gain_`, `$hydraulic_cost_norm_` and `$thermal_cost_`, +##' which are unitless and are NOT the same quantities as `$hydraulic_cost_`. +##' ##' All water potentials are **positive magnitudes in MPa**. There is one ##' representation throughout and it is asserted, not documented: a negative ##' `psi_crit`, `stem_b`, `root_b` or `root_psi_crit` is an error here, and a @@ -141,6 +149,9 @@ NULL electron_transport = function() { Leaf__electron_transport(self) }, + leaf_temp_from_E = function(E) { + Leaf__leaf_temp_from_E(self, E) + }, set_leaf_states_rates_from_psi_stem = function(psi_stem, psi_upstream) { Leaf__set_leaf_states_rates_from_psi_stem(self, psi_stem, psi_upstream) }, @@ -186,6 +197,21 @@ NULL optimise_psi_stem_Sperry = function() { Leaf__optimise_psi_stem_Sperry(self) }, + thermal_cost_at = function(leaf_temp) { + Leaf__thermal_cost_at(self, leaf_temp) + }, + prepare_profitmax = function() { + Leaf__prepare_profitmax(self) + }, + profit_psi_stem_ProfitMax = function(psi_stem, psi_upstream) { + Leaf__profit_psi_stem_ProfitMax(self, psi_stem, psi_upstream) + }, + profitmax_curve = function(n) { + Leaf__profitmax_curve(self, n) + }, + optimise_psi_stem_ProfitMax = function() { + Leaf__optimise_psi_stem_ProfitMax(self) + }, optimise_psi_stem_TF = function() { Leaf__optimise_psi_stem_TF(self) }, @@ -358,6 +384,27 @@ NULL Leaf__lambda___set(self, value) } }, + carbon_gain_ = function(value) { + if (missing(value)) { + Leaf__carbon_gain___get(self) + } else { + Leaf__carbon_gain___set(self, value) + } + }, + hydraulic_cost_norm_ = function(value) { + if (missing(value)) { + Leaf__hydraulic_cost_norm___get(self) + } else { + Leaf__hydraulic_cost_norm___set(self, value) + } + }, + thermal_cost_ = function(value) { + if (missing(value)) { + Leaf__thermal_cost___get(self) + } else { + Leaf__thermal_cost___set(self, value) + } + }, lambda_analytical_ = function(value) { if (missing(value)) { Leaf__lambda_analytical___get(self) @@ -722,6 +769,41 @@ NULL Leaf__use_energy_balance___set(self, value) } }, + vpd_leaf_ = function(value) { + if (missing(value)) { + Leaf__vpd_leaf___get(self) + } else { + Leaf__vpd_leaf___set(self, value) + } + }, + use_thermal_cost_ = function(value) { + if (missing(value)) { + Leaf__use_thermal_cost___get(self) + } else { + Leaf__use_thermal_cost___set(self, value) + } + }, + T50_ = function(value) { + if (missing(value)) { + Leaf__T50___get(self) + } else { + Leaf__T50___set(self, value) + } + }, + Tcrit_ = function(value) { + if (missing(value)) { + Leaf__Tcrit___get(self) + } else { + Leaf__Tcrit___set(self, value) + } + }, + profitmax_scan_n_ = function(value) { + if (missing(value)) { + Leaf__profitmax_scan_n___get(self) + } else { + Leaf__profitmax_scan_n___set(self, value) + } + }, d_ = function(value) { if (missing(value)) { Leaf__d___get(self) diff --git a/README.md b/README.md index 613c3ea..b7c1fc1 100644 --- a/README.md +++ b/README.md @@ -25,10 +25,13 @@ so it can be tested, profiled, extended and embedded on its own. **A home for several stomatal models, not just ours.** The package carries our hydraulic gain-risk formulation, and it already contains two alternatives — the Sperry et al. (2017) cost formulation and the Medlyn et al. (2011) optimal -stomatal model — inherited from plant. Today those are second-class: -`optimise_psi_stem_Sperry` is hardwired to a single soil layer, and the Medlyn path -bypasses the hydraulic solve altogether, so neither can be swapped in as a -like-for-like alternative. The goal is to make each one a **first-class member**, +stomatal model — inherited from plant. The Sperry side is now runnable as the paper +defines it: `optimise_psi_stem_ProfitMax()` maximises `CG - HC` with both terms +normalised, on the single-potential supply path, and `optimise_psi_stem_TF()` runs +the TF24 cost through the same solver so the two can be compared at identical +drivers. Both are still off the production collar solve, and the Medlyn path +bypasses the hydraulic solve altogether. The goal is to make each one a +**first-class member** there too, selectable and runnable against identical drivers, alongside Prentice et al. (2014) least-cost and Cowan-Farquhar. diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index c84c925..d3d5fcb 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -80,6 +80,14 @@ Leaf: `$opt_psi_stem_`, `$profit_`, `$assim_colimited_`, `$transpiration_`, `$lambda`, `$g1_eff`. + A second optimality formulation is available on the single-soil-layer + supply path: `$optimise_psi_stem_ProfitMax()` maximises Sperry (2017)'s + normalised profit, and `$profitmax_curve()` returns the whole cost, gain + and profit curve in one call. `$optimise_psi_stem_TF()` runs the TF24 cost + through the same solver, so the two can be compared at identical drivers. + These write `$carbon_gain_`, `$hydraulic_cost_norm_` and `$thermal_cost_`, + which are unitless and are NOT the same quantities as `$hydraulic_cost_`. + All water potentials are **positive magnitudes in MPa**. There is one representation throughout and it is asserted, not documented: a negative `psi_crit`, `stem_b`, `root_b` or `root_psi_crit` is an error here, and a @@ -171,6 +179,14 @@ Leaf: profit_: {type: double, access: field} psi_stem: {type: double, access: field} lambda_: {type: double, access: field} + # --- Sperry (2017) ProfitMax, all unitless ------------------------------- + # Written ONLY by the ProfitMax path. `hydraulic_cost_` is deliberately NOT + # reused: it holds carbon units on the TF24 path and conductance units on the + # Sperry-cost one, and a third meaning behind the same name is how a reader + # quotes the wrong number. + carbon_gain_: {type: double, access: field} + hydraulic_cost_norm_: {type: double, access: field} + thermal_cost_: {type: double, access: field} lambda_analytical_: {type: double, access: field} hydraulic_cost_: {type: double, access: field} opt_psi_stem_: {type: double, access: field} @@ -266,6 +282,20 @@ Leaf: # --- energy balance (#523) ------------------------------------------------ use_energy_balance_: {type: bool, access: field} + # The deficit Fick's law divides by (PLAN 13.1). READ-ONLY in practice: it is + # derived by set_leaf_vpd from `atm_vpd` and the leaf temperature, and every + # solve overwrites it. Bound because a caller checking why gs moved needs to + # see it, and because `vpd_leaf_ == atm_vpd_` is the statement that the + # prescribed-temperature path is unchanged. + vpd_leaf_: {type: double, access: field} + # --- Sicangco et al. (2026) instantaneous thermal cost, default OFF ------- + # ⚠️ NOT plant's ATLS thermal damage, which is a LASTING ratchet. This is a + # cost paid at the temperature the leaf is at now, with no memory. + use_thermal_cost_: {type: bool, access: field} + T50_: {type: double, access: field} + Tcrit_: {type: double, access: field} + # Points used to scan the supply stream for |A|max. Sicangco use 500-600. + profitmax_scan_n_: {type: int, access: field} d_: {type: double, access: field} wind_speed_: {type: double, access: field} Tair_: {type: double, access: field} @@ -416,6 +446,14 @@ Leaf: electron_transport: return_type: double args: [] + # The operating-point leaf temperature at a given transpiration, on the + # energy-balance path. Bound so a caller aligning an EXTERNAL energy balance to + # this one can iterate on the model's own function rather than transcribing + # `Tair + (Rn - lambda*E)*ra/rho_cp` into R -- where it would go stale the + # first time this one changed, silently. + leaf_temp_from_E: + return_type: double + args: [E: double] set_leaf_states_rates_from_psi_stem: return_type: void args: [psi_stem: double, psi_upstream: double] @@ -474,6 +512,39 @@ Leaf: optimise_psi_stem_Sperry: return_type: void args: [] + + # --- Sperry (2017) ProfitMax with BOTH terms normalised ------------------- + # + # `optimise_psi_stem_Sperry` maximises `A - lambda*cost` with lambda a + # PRESCRIBED input; these maximise `CG - (HC + TC)` with CG and HC normalised + # as Sperry defines them. The two share an argmax exactly when + # lambda = |A|max/(k_soil - kcrit), and `optimise_psi_stem_ProfitMax` reports + # that lambda in `$lambda_` so the equivalence can be checked rather than + # believed. Use these unless you have a lambda of your own: |A|max is a + # property of the drivers, so the equivalent lambda has to be recomputed at + # every observation and doing that from R costs more than the model does. + # + # Single soil layer only, like the other two -- which is not a restriction for + # this model, whose supply is the integral of one vulnerability curve with no + # root network in it. + thermal_cost_at: + return_type: double + args: [leaf_temp: double] + prepare_profitmax: + return_type: void + args: [] + profit_psi_stem_ProfitMax: + return_type: double + args: [psi_stem: double, psi_upstream: double] + # n rows of (psi, CG, HC, TC, profit), flattened COLUMN-MAJOR, in one crossing + # of the R boundary. This is Sicangco et al.'s Figures 2, 3 and S4; building it + # row by row from R would pay ~1.8 us of call overhead against a ~3 us model. + profitmax_curve: + return_type: "std::vector" + args: [n: int] + optimise_psi_stem_ProfitMax: + return_type: void + args: [] optimise_psi_stem_TF: return_type: void args: [] diff --git a/inst/include/phylloptim/constants.hpp b/inst/include/phylloptim/constants.hpp index f8355d6..36f9235 100644 --- a/inst/include/phylloptim/constants.hpp +++ b/inst/include/phylloptim/constants.hpp @@ -127,6 +127,15 @@ inline constexpr double aerodynamic_resistance_coef = 200.0; // profit and are simply rejected. inline constexpr double leaf_temp_min = -40.0; inline constexpr double leaf_temp_max = 70.0; +// Floor on the leaf-to-air vapour pressure deficit (kPa) used by Fick's law +// (PLAN 13.1). The same argument as the temperature clamp above, one step +// downstream: a leaf cooled below the dew point has a NEGATIVE deficit, and +// dividing a positive transpiration by it reports a negative stomatal +// conductance. Flooring keeps such points finite and unattractive rather than +// sign-inverted. 0.01 kPa is far below any operating deficit -- the driver +// defaults are 1.5-2 kPa -- so it binds only where the model is already outside +// what a one-way diffusion equation describes. +inline constexpr double vpd_leaf_min = 0.01; } // namespace phylloptim diff --git a/inst/include/phylloptim/leaf_model.hpp b/inst/include/phylloptim/leaf_model.hpp index edd4959..3cfbe8d 100644 --- a/inst/include/phylloptim/leaf_model.hpp +++ b/inst/include/phylloptim/leaf_model.hpp @@ -344,8 +344,42 @@ class Leaf { // set_physiology falls back to the fixed ra. Only read on the PM path. double d_ = 0.05; // characteristic leaf dimension, m double wind_speed_ = 2.0; // above-canopy wind speed U0, m s^-1 + + // --- Sicangco et al. (2026) thermal cost, default OFF ----------------------- + // An instantaneous cost of PSII damage, sigmoid in LEAF temperature: + // + // TC(Tleaf) = 1 / (1 + exp(-r*(Tleaf - T50))), r = 2/(T50 - Tcrit) + // + // and, in the same paper's ProfitMaxTC, Jmax is scaled by (1 - TC) so electron + // transport is destroyed at the same rate the cost is incurred. Both halves are + // behind this one gate. + // + // ⚠️ IT IS INSTANTANEOUS, and that is a modelling position rather than an + // omission. plant's ATLS layer models thermal damage as a LASTING ratchet -- + // damage accumulated at one time step is still there at the next. This is the + // other thing: a cost paid at the temperature the leaf is at right now, with no + // memory. They are not interchangeable and neither is a port of the other. + // + // ⚠️ AND THE PAPER AND ITS CODE DISAGREE ON r. The published Eqn 8 area reads + // "r = (T50 - Tcrit)/2"; gsthermal's F0_func and TJmax_updated both compute + // `r = 2/(T50 - Tcrit)`. The code is what produced the figures, so the code is + // what is implemented here. + bool use_thermal_cost_ = false; + double T50_ = 50.4; // deg C; Sicangco Table 2, heatwave treatment + double Tcrit_ = 46.5; // deg C; Sicangco Table 2, heatwave treatment + + // Points used to scan the transpiration supply stream for |A|max, which is what + // Sperry's carbon gain is normalised by. Sicangco's Ps_to_Pcrit defaults to 500 + // and their instantaneous simulations use 600. A member rather than a + // constructor argument on purpose: the constructor's arity is pinned by plant's + // generated RcppR6 glue and by the CI consumer program. + int profitmax_scan_n_ = 500; double PPFD_; double atm_vpd_; + // The vapour pressure deficit the DIFFUSION equations use, kPa. Off the + // energy-balance path it is exactly `atm_vpd_`; on it, the leaf-to-air deficit + // at the operating-point leaf temperature (PLAN 13.1, #7). See set_leaf_vpd. + double vpd_leaf_; double atm_o2_kpa_; // --- Temperature-response parameters --------------------------------------- @@ -413,6 +447,28 @@ class Leaf { double opt_root_psi_; double assim_max_; + // --- Sperry (2017) ProfitMax outputs, all unitless ------------------------- + // Written ONLY by the ProfitMax path (optimise_psi_stem_ProfitMax and + // profit_psi_stem_ProfitMax). Separate members rather than reusing + // `hydraulic_cost_`, which is in carbon units on the TF24 path and in + // conductance units on the Sperry-cost one: three meanings behind one name is + // hazard 8 waiting to happen. + // Normalisers, seeded by prepare_profitmax() and read by the two functions + // above. Cached rather than recomputed per candidate because they are constant + // over one solve and each costs an exp+pow. + double profitmax_A_max_; // |A|max over the supply stream, umol m^-2 s^-1 + double profitmax_k_soil_; // k(psi_soil), kg m^-2 s^-1 MPa^-1 + double profitmax_k_span_; // k(psi_soil) - k_crit, the HC denominator + // The scan prepare_profitmax() already runs, kept so the objective can be + // rebuilt on it without evaluating the model a second time. See + // optimise_psi_stem_ProfitMax for why a grid is needed at all. + std::vector profitmax_scan_psi_; + std::vector profitmax_scan_A_; + std::vector profitmax_scan_Tleaf_; + double carbon_gain_; // CG = A/|A|max over the supply stream + double hydraulic_cost_norm_; // HC = [k(psi_soil)-k(psi)]/[k(psi_soil)-k_crit] + double thermal_cost_; // TC, zero unless use_thermal_cost_ + double opt_psi_stem_; double opt_ci_; @@ -482,7 +538,7 @@ class Leaf { // * the cost is 17 double comparisons per set_physiology() call, i.e. per // driver set, NOT per inner solve iteration. Measured: within run-to-run // noise on bench_solve. - static constexpr int photo_temp_key_size = 19; + static constexpr int photo_temp_key_size = 22; std::array photo_temp_cache_key_{}; std::array photo_temp_key() const; std::vector f_r; @@ -831,9 +887,9 @@ class Leaf { // gate-off path; the derivation and the two sign checks are at the definition. double dprofit_energy_balance_term(double ci, double gc, double g_ci, double inv_atm, double gc_const, - double dgc_dpsistem, double dgc_dpsi, - double dpsistem_dpsi, double dT_dE, - double Tleaf); + double A_prime, double dgc_dpsistem, + double dgc_dpsi, double dpsistem_dpsi, + double dT_dE, double Tleaf); // The profit-maximising collar potential within [bound_a, bound_b], by a // safeguarded root-find on dprofit == 0 (PLAN 11a). Returns a bound when the // optimum is pinned to it, which is the case on 42 of the 240 feasible @@ -881,6 +937,19 @@ class Leaf { // kept); provided and unit-tested for the leaf-to-air VPD in the full cut. double saturation_vapour_pressure(double temp) const; double saturation_vapour_pressure_slope(double temp) const; + + // Seat `vpd_leaf_` at a leaf temperature. THE ONE PLACE the leaf-to-air deficit + // is defined, called wherever a leaf temperature is established: set_physiology + // (at Tair, so off-path callers get atm_vpd_ back exactly), + // set_leaf_states_rates_from_psi_stem, and dprofit_at_collar_psi. + // + // ⚠️ NOT folded into update_temperature_dependent_params, even though that is + // the other function taking a leaf temperature. The two finite-difference + // blocks call that one at T +/- h and restore the photosynthesis members by + // assignment; `vpd_leaf_` must NOT move with them, because dA/dTleaf there is + // taken at fixed ci and assim_colimited() reads no VPD at all. The VPD route + // into the derivative is carried separately, by dgc_dT. + void set_leaf_vpd(double leaf_temp); // Explicit leaf energy balance: Tleaf = Tair + (Rn - lambda*E) * ra / (rho*cp). // E is the hydraulically-pinned transpiration (kg H2O m^-2 s^-1); no PM // inversion and no A->E feedback, so this is a single algebraic forward pass. @@ -1090,9 +1159,34 @@ class Leaf { double profit_psi_stem_Sperry(double psi_stem, double psi_upstream); double profit_psi_stem_TF(double psi_stem, double psi_upstream); + // The instantaneous thermal cost at a leaf temperature, in [0,1]. Zero when the + // gate is off, so callers need not branch. + double thermal_cost_at(double leaf_temp) const; + + // Sperry (2017) ProfitMax, with BOTH terms normalised as the paper defines them + // -- see optimise_psi_stem_ProfitMax for what that buys over the lambda form. + // Seeds |A|max and the conductance span; the two below read what it seeds. + void prepare_profitmax(); + double profit_psi_stem_ProfitMax(double psi_stem, double psi_upstream); + + // The whole cost/gain/profit curve over [psi_soil, psi_crit] in ONE crossing of + // the R boundary: n rows of (psi, CG, HC, TC, profit), flattened column-major. + // This is Sicangco et al.'s Figures 2, 3 and S4, and building it row by row from + // R would pay ~1.8 us of call overhead against a ~3 us model evaluation. + std::vector profitmax_curve(int n); + // optimiser functions void optimise_psi_stem_Sperry(); void optimise_psi_stem_TF(); + void optimise_psi_stem_ProfitMax(); + + // Clear the outputs a single-layer optimiser does NOT write. Hazard 8: these + // three describe a ROOT-COLLAR solve, and optimise_psi_stem_* never runs one, + // so leaving the last find_root_collar_psi()'s values standing reports an + // operating point whose parts came from two different solves. Measured before + // this existed: E = 9.216e-5 from the Sperry solve sitting beside E_up = + // 2.626e-5 from a collar solve several calls earlier. + void clear_collar_solve_state(); // --- WHICH KIND of operating point the collar solve found ------------------- // @@ -1387,6 +1481,12 @@ inline void Leaf::setup_clean_leaf() { transpiration_= util::na_value; // kg m^-2 s^-1 profit_= util::na_value; // umol C m^-2 s^-1 lambda_= util::na_value; // umol C m^-2 s^-1 kg^-1 m^2 s^1 + carbon_gain_= util::na_value; + hydraulic_cost_norm_= util::na_value; + thermal_cost_= util::na_value; + profitmax_A_max_= util::na_value; + profitmax_k_soil_= util::na_value; + profitmax_k_span_= util::na_value; lambda_analytical_= util::na_value; // umol C m^-2 s^-1 kg^-1 m^2 s^1 hydraulic_cost_= util::na_value; // umol C m^-2 s^-1 electron_transport_= util::na_value; //electron transport rate umol m^-2 s^-1 @@ -1408,7 +1508,8 @@ inline void Leaf::setup_clean_leaf() { Rn_= util::na_value; // W m^-2 ra_= util::na_value; // s m^-1 PPFD_= util::na_value; //umol m^-2 s^-1 - atm_vpd_= util::na_value; //kPa + atm_vpd_= util::na_value; //kPa + vpd_leaf_= util::na_value; //kPa atm_o2_kpa_= util::na_value; // kPa atm_kpa_= util::na_value; // kPa umol_per_mol_to_Pa_ = util::na_value; @@ -1533,6 +1634,11 @@ inline void Leaf::set_physiology(const RootNetwork& root_network, double PPFD, c ? aerodynamic_resistance_coef * std::sqrt(d_ / wind_speed_) : aerodynamic_resistance_fixed; + // The diffusion deficit at the baseline temperature. On the PM path the solve + // re-seats this per candidate psi; here it makes assim_max_ and any bare + // stom_cond_CO2() call well defined before one has run. + set_leaf_vpd(leaf_temp_); + // Temperature/O2-dependent block. Off the PM path this is recomputed only // when (leaf_temp_, atm_o2_kpa_) changes from the previous call (see // photo_temp_cache_ in the header); same inputs -> bit-identical outputs, so @@ -2305,6 +2411,7 @@ inline double Leaf::dprofit_at_collar_psi(double opt_root_psi, bool* feasible) { if (use_energy_balance_) { Tleaf_here = leaf_temp_from_E(transpiration(psi_stem, psi), &dT_dE); update_temperature_dependent_params(Tleaf_here); + set_leaf_vpd(Tleaf_here); } const double ci = psi_stem_to_ci(psi_stem, psi); @@ -2389,15 +2496,15 @@ inline double Leaf::dprofit_at_collar_psi(double opt_root_psi, bool* feasible) { vcmax_ = vc0; jmax_ = jm0; gamma_ = ga0; ko_ = ko0; kc_ = kc0; R_d_ = rd0; km_ = km0; electron_transport_ = J0; const double Rd_T = (Rd_up - Rd_dn) / (2.0 * h); - // gc = gc_const * E, so dE/dpsi comes from the same spline derivatives the - // main branch uses; recomputed here because the main branch's locals are - // below this early return. - const double gc_c = atm_kpa_ * kg_to_mol_h2o / atm_vpd_ / H2O_CO2_stom_diff_ratio_; - const double dgc_ps = gc_c * leaf_specific_conductance_max_ * - stem_curve_integral_deriv(psi_stem); - const double dgc_p = gc_c * leaf_specific_conductance_max_ * - (-stem_curve_integral_deriv(psi)); - const double dE_dpsi = (dgc_ps * dpsistem_dpsi + dgc_p) / gc_c; + // dE/dpsi from the same spline derivatives the main branch uses; + // recomputed here because the main branch's locals are below this early + // return. It used to be written as the two gc partials divided back by + // their shared constant, which cancelled a coefficient that is no longer a + // single number now that the deficit moves with Tleaf. Stated directly. + const double dE_dpsi = + leaf_specific_conductance_max_ * + (stem_curve_integral_deriv(psi_stem) * dpsistem_dpsi - + stem_curve_integral_deriv(psi)); dprofit += -Rd_T * dT_dE * dE_dpsi; } return std::isfinite(dprofit) ? dprofit : 0.0; @@ -2422,7 +2529,7 @@ inline double Leaf::dprofit_at_collar_psi(double opt_root_psi, bool* feasible) { // (transp_from_psi(psi_stem) - transp_from_psi(psi)), so the partials use the // analytic spline derivative. const double gc_const = - atm_kpa_ * kg_to_mol_h2o / atm_vpd_ / H2O_CO2_stom_diff_ratio_; + atm_kpa_ * kg_to_mol_h2o / vpd_leaf_ / H2O_CO2_stom_diff_ratio_; const double gc = gc_const * transpiration(psi_stem, psi); const double dgc_dpsistem = gc_const * leaf_specific_conductance_max_ * stem_curve_integral_deriv(psi_stem); @@ -2455,7 +2562,7 @@ inline double Leaf::dprofit_at_collar_psi(double opt_root_psi, bool* feasible) { return base; } return base + dprofit_energy_balance_term(ci, gc, g_ci, inv_atm, gc_const, - dgc_dpsistem, dgc_dpsi, + A_prime, dgc_dpsistem, dgc_dpsi, dpsistem_dpsi, dT_dE, Tleaf_here); } @@ -2499,14 +2606,31 @@ inline double Leaf::dprofit_at_collar_psi(double opt_root_psi, bool* feasible) { // cooling back toward the optimum. That is the decoupling mechanism, and its // sign is the sharpest available test that this term is right. // -// `dgc_dT` is a named zero rather than an omission: stom_cond_CO2 currently -// divides by the prescribed AIR vpd, so gc does not depend on Tleaf. Wiring -// leaf-to-air VPD (PLAN 13.1) makes it non-zero and this becomes a one-line -// change instead of a re-derivation. +// ⚠️ `dgc_dT` IS NOW NON-ZERO, and the placeholder algebra that stood here while +// it was zero did NOT generalise. This block used to carry a named +// `const double dgc_dT = 0.0` with a note promising that PLAN 13.1 would make it +// "a one-line change instead of a re-derivation". The note was wrong: the +// damping factor it multiplied, `(gc*inv_atm - dgc_dT*(ca-ci)*inv_atm)/g_ci`, +// puts the new term under `A_T` where the derivation puts it under `A_prime`. +// The two agree only at dgc_dT = 0, which is why nothing caught it. Derived +// again, in full, below. +// +// g(ci; psi_stem, psi, T) = A(ci,T)*umol_to_mol - gc(T)*(ca-ci)*inv_atm = 0, so +// +// g_T = A_T*umol_to_mol - dgc_dT*(ca-ci)*inv_atm +// dci/dT = -g_T/g_ci +// dA/dT (total) = A_T + A_prime*dci/dT +// = [A_T*gc*inv_atm + A_prime*dgc_dT*(ca-ci)*inv_atm] / g_ci +// +// using g_ci = A_prime*umol_to_mol + gc*inv_atm to collect. At dgc_dT = 0 this is +// `A_T * (gc*inv_atm)/g_ci`, i.e. exactly what the old expression returned, so +// the prescribed-VPD behaviour is unchanged. The two checks in the paragraph +// above still hold for the first term; the second is new and has its own sign +// argument at the assignment. inline double Leaf::dprofit_energy_balance_term( double ci, double gc, double g_ci, double inv_atm, double gc_const, - double dgc_dpsistem, double dgc_dpsi, double dpsistem_dpsi, double dT_dE, - double Tleaf) { + double A_prime, double dgc_dpsistem, double dgc_dpsi, double dpsistem_dpsi, + double dT_dE, double Tleaf) { if (!use_energy_balance_ || dT_dE == 0.0) { return 0.0; } @@ -2558,10 +2682,19 @@ inline double Leaf::dprofit_energy_balance_term( const double dE_dpsi = (dgc_dpsistem * dpsistem_dpsi + dgc_dpsi) / gc_const; const double tau = dT_dE * dE_dpsi; - const double dgc_dT = 0.0; // see the note above; non-zero once PLAN 13.1 lands - const double damping = (gc * inv_atm - dgc_dT * (ca_ - ci) * inv_atm) / g_ci; + // dgc/dTleaf through the deficit alone. gc = K*E/D(T) with E hydraulically + // pinned, so dgc/dT = -gc * D'(T)/D, and D'(T) = esat'(Tleaf) since e_air is + // fixed. NEGATIVE: a hotter leaf has a larger deficit, so the same water flux + // implies a SMALLER conductance. Zero where the deficit is on its floor, for + // the same reason leaf_temp_from_E returns dT_dE = 0 on its clamp -- the + // reported value no longer responds to the input. + const double dgc_dT = + (vpd_leaf_ > vpd_leaf_min) + ? -gc * saturation_vapour_pressure_slope(Tleaf) / vpd_leaf_ + : 0.0; - return A_T * tau * damping; + return tau * + (A_T * gc * inv_atm + A_prime * dgc_dT * (ca_ - ci) * inv_atm) / g_ci; } inline double Leaf::arrh_curve(double Ea, double ref_value, double leaf_temp) const { @@ -2607,13 +2740,23 @@ inline std::array Leaf::photo_temp_key() cons gamma_25_, gamma_ha_, kc_25_, kc_ha_, ko_25_, ko_ha_, - R_d_25, rd_q10_intercept_, rd_q10_slope_}; + R_d_25, rd_q10_intercept_, rd_q10_slope_, + // The thermal cost reaches jmax_ above, so it belongs in the key. A + // bool widens to 0.0/1.0, which compares by == exactly like the rest. + double(use_thermal_cost_), T50_, Tcrit_}; } inline void Leaf::update_temperature_dependent_params(double leaf_temp) { vcmax_ = peak_arrh_curve(vcmax_ha_, vcmax_25, leaf_temp, vcmax_H_d_, vcmax_d_S_); jmax_ = peak_arrh_curve(jmax_ha_, jmax_25, leaf_temp, jmax_H_d_, jmax_d_S_); + // Sicangco et al. (2026) Eqn 12: irreversible PSII damage takes electron + // transport down with it, so the peaked Arrhenius Jmax is scaled by (1 - TC). + // Branched rather than multiplied by a gate-off zero, so the gate-off + // arithmetic is untouched. + if (use_thermal_cost_) { + jmax_ *= 1.0 - thermal_cost_at(leaf_temp); + } gamma_ = arrh_curve(gamma_ha_, gamma_25_, leaf_temp); ko_ = arrh_curve(ko_ha_, ko_25_, leaf_temp); kc_ = arrh_curve(kc_ha_, kc_25_, leaf_temp); @@ -2644,6 +2787,39 @@ inline double Leaf::saturation_vapour_pressure_slope(double temp) const { return 4098.0 * saturation_vapour_pressure(temp) / ((temp + 237.3) * (temp + 237.3)); } +// The deficit that drives diffusion out of the leaf (PLAN 13.1, issue #7). +// +// Fick's law needs the deficit between the SUB-STOMATAL cavity, saturated at the +// LEAF temperature, and the air. The driver `atm_vpd_` is the deficit at AIR +// temperature, so the two agree only when the leaf is at air temperature: +// +// e_air = esat(Tair) - atm_vpd +// vpd_leaf = esat(Tleaf) - e_air = atm_vpd + [esat(Tleaf) - esat(Tair)] +// +// ⚠️ WRITTEN AS `atm_vpd_ + (esat(T) - esat(Tair_))` FOR AN ARITHMETIC REASON. +// Off the energy-balance path the leaf IS at air temperature, so the bracket is +// a value minus itself: exactly +0.0, and `x + 0.0 == x` for every finite x. The +// prescribed-temperature path is therefore bit-identical, which is what keeps +// the golden file untouched. Computing it as `esat(T) - esat(Tair) + atm_vpd` in +// a different association would not guarantee that. +// +// ⚠️ AND IT CAN GO NON-POSITIVE. A leaf transpiring hard enough to sit well below +// air temperature has esat(Tleaf) < e_air, i.e. condensation rather than +// evaporation, which Fick's law in this direction does not describe -- and a +// deficit passing through zero is an infinite conductance. Clamped to +// `vpd_leaf_min` rather than allowed to invert the flux: the model would +// otherwise report a NEGATIVE stomatal conductance for a positive transpiration. +inline void Leaf::set_leaf_vpd(double leaf_temp) { + if (!use_energy_balance_) { + vpd_leaf_ = atm_vpd_; + return; + } + const double d = + atm_vpd_ + (saturation_vapour_pressure(leaf_temp) - + saturation_vapour_pressure(Tair_)); + vpd_leaf_ = std::max(d, vpd_leaf_min); +} + // Explicit leaf energy balance (#523): Tleaf = Tair + (Rn - lambda*E)*ra/(rho*cp). // E is the hydraulically-pinned transpiration (kg H2O m^-2 s^-1), so lambda*E is // the latent heat flux (W m^-2) and (Rn - lambda*E) the sensible heat flux H. @@ -2896,7 +3072,7 @@ inline double Leaf::transpiration_to_psi_stem(double transpiration_, double psi_ // returns stomatal conductance to CO2, mol C m^-2 LA s^-1 inline double Leaf:: stom_cond_CO2(double psi_stem, double psi_upstream) { double transpiration_ = transpiration(psi_stem, psi_upstream); - return atm_kpa_ * transpiration_ * kg_to_mol_h2o / atm_vpd_ / H2O_CO2_stom_diff_ratio_; + return atm_kpa_ * transpiration_ * kg_to_mol_h2o / vpd_leaf_ / H2O_CO2_stom_diff_ratio_; } @@ -3082,24 +3258,50 @@ inline double Leaf::psi_stem_to_ci(double psi_stem, double psi_upstream) { // given psi_stem, find assimilation, transpiration and stomal conductance to c02 inline void Leaf::set_leaf_states_rates_from_psi_stem(double psi_stem, double psi_upstream) { - // ⚠️ `Tleaf_` IS WRITTEN ON ALL THREE PATHS BELOW, which is hazard 8 applied to - // the one output whose value is a driver on one path and a solved quantity on - // the other. The two shut-down branches transpire nothing, and on the PM path a - // leaf that transpires nothing is the HOTTEST one -- so leaving `Tleaf_` to the - // third branch would report the previous candidate's temperature for exactly - // the operating points where the answer is most extreme. + // ⚠️ AN `assim_max_ < 0` EARLY EXIT USED TO SIT HERE, AND IT ZEROED THE WATER. + // It read: if gross assimilation at ci = ca cannot cover respiration, put ci at + // the compensation point and set transpiration and conductance to zero. The + // first half is right and is now the ci solver's own job; the second half is + // wrong, and wrong in a way that matters. + // + // Transpiration on this path is the HYDRAULIC SUPPLY at the candidate potential. + // It does not ask whether there is carbon to be had -- water moves down a + // potential gradient whatever photosynthesis is doing. Zeroing it made two leaves + // at the SAME operating point disagree about whether water was moving purely + // because one of them had respiration switched on, which is how it was found. + // + // Deleting it changes nothing about the shut-down STATE it was reaching for. + // Where `assim_max_ < 0` the ci root-find has no supply==demand root anywhere in + // [gamma*, ca] -- assimilation is monotone in ci and negative at both ends -- so + // psi_stem_to_ci takes its compensation-point fallback and returns ci = gamma* + // with gross assimilation zero and net = -R_d, which is exactly what this branch + // set by hand. That fallback did not exist when the branch was written; it was + // added with the Q10 respiration in #41, which is also what made this regime + // reachable at ordinary light. + // + // ⚠️ AND IT IS UNREACHABLE FROM THE COLLAR SOLVE, which is why plant is + // unaffected: prepare_collar_solve has its OWN `assim_max_ < 0` exit and returns + // false before any candidate potential is evaluated. The golden grid's minimum + // assim_max_ is 3.71, so the file never reaches this branch either and is + // bit-identical across the change -- which is a statement about the grid, not + // evidence that the change is inert. test_transpiration_survives_negative_assim + // is the test that can see it. + // + // ⚠️ `Tleaf_` IS WRITTEN ON BOTH PATHS BELOW, which is hazard 8 applied to the + // one output whose value is a driver on one path and a solved quantity on the + // other. The remaining shut-down branch transpires nothing, and on the PM path + // a leaf that transpires nothing is the HOTTEST one -- so leaving `Tleaf_` to + // the transpiring branch would report the previous candidate's temperature for + // exactly the operating points where the answer is most extreme. (#106 wrote + // this on all THREE branches it found; deleting the `assim_max_ < 0` exit + // leaves two.) if (psi_upstream >= psi_stem){ ci_ = gamma_*umol_per_mol_to_Pa_; transpiration_ = 0; stom_cond_CO2_ = 0; Tleaf_ = use_energy_balance_ ? leaf_temp_from_E(0.0) : leaf_temp_; } else{ - if(assim_max_ < 0){ - ci_ = gamma_*umol_per_mol_to_Pa_; - transpiration_ = 0; - stom_cond_CO2_ = 0; - Tleaf_ = use_energy_balance_ ? leaf_temp_from_E(0.0) : leaf_temp_; - } else{ + { // Transpiration is the hydraulic supply, independent of ci; compute it // first so the PM path can derive the operating-point leaf temperature. // Off the PM path this is a memoised no-op reorder (psi_stem_to_ci -> @@ -3118,9 +3320,15 @@ inline void Leaf::set_leaf_states_rates_from_psi_stem(double psi_stem, double ps Tleaf_ = use_energy_balance_ ? leaf_temp_from_E(transpiration_) : leaf_temp_; if (use_energy_balance_) { update_temperature_dependent_params(Tleaf_); + // ...and the deficit Fick's law divides by, which moves with it. Both + // read the STORED `Tleaf_` rather than re-deriving it, so the temperature + // the Farquhar parameters were computed at, the one the deficit was + // computed at, and the one reported to the caller cannot be three + // different numbers -- and it costs no extra `leaf_temp_from_E`. + set_leaf_vpd(Tleaf_); } ci_ = psi_stem_to_ci(psi_stem, psi_upstream); - stom_cond_CO2_ = atm_kpa_ * transpiration_ * kg_to_mol_h2o / atm_vpd_ / H2O_CO2_stom_diff_ratio_; + stom_cond_CO2_ = atm_kpa_ * transpiration_ * kg_to_mol_h2o / vpd_leaf_ / H2O_CO2_stom_diff_ratio_; } } assim_colimited_ = assim_colimited(ci_); @@ -3188,7 +3396,7 @@ inline double Leaf::g1_eff() const { if (!std::isfinite(chi) || chi >= 1.0) { return util::na_value; } - return chi * std::sqrt(atm_vpd_) / (1.0 - chi); + return chi * std::sqrt(vpd_leaf_) / (1.0 - chi); } // Pure: no write to hydraulic_cost_, so the AD pass cannot scribble model state @@ -3205,6 +3413,24 @@ inline double Leaf::hydraulic_cost_TF(double psi_stem) { return hydraulic_cost_; } +// Sicangco et al. (2026) Eqns 8-9: the fraction of maximum PSII damage sustained +// at this leaf temperature, in [0,1]. `r` follows gsthermal, not the printed Eqn +// 8 -- see the note on use_thermal_cost_. +// +// Exactly 0.0 with the gate off, so every caller can add it unconditionally. +inline double Leaf::thermal_cost_at(double leaf_temp) const { + if (!use_thermal_cost_) { + return 0.0; + } + const double span = T50_ - Tcrit_; + if (!(span > 0.0)) { + util::stop("thermal cost needs T50 > Tcrit; got T50 = " + + util::format_double(T50_) + ", Tcrit = " + + util::format_double(Tcrit_)); + } + return 1.0 / (1.0 + std::exp(-(2.0 / span) * (leaf_temp - T50_))); +} + // Profit functions inline double Leaf::profit_psi_stem_Sperry(double psi_stem, double psi_upstream) { @@ -3230,19 +3456,53 @@ double benefit_ = assim_colimited_; //optimisation functions +// Everything a single-layer optimiser leaves untouched, cleared rather than +// inherited. See the declaration for what it cost when this did not exist. +inline void Leaf::clear_collar_solve_state() { + operating_point_kind_ = OperatingPointKind::Unsolved; + opt_root_psi_ = util::na_value; + E_up_ = util::na_value; + soil_consumption_.assign(soil_consumption_.size(), util::na_value); +} + +// ⚠️ THE TWO OPTIMISERS BELOW USE A BARE BRENT SEARCH AND SHARE A HAZARD. +// Brent is a LOCAL optimiser that steps in from the bounds, so it cannot return +// an endpoint and it cannot see past a local maximum. `optimise_psi_stem_TF` is +// therefore wrong wherever the TF24 profit is maximised at full closure, and +// `optimise_psi_stem_Sperry` wherever the caller's lambda makes it so. +// `optimise_psi_stem_ProfitMax` scans a grid first for exactly this reason -- see +// the note there, with the measured case that motivated it. These two are left +// alone because neither has a scan to reuse, so fixing them costs 500 model +// evaluations per solve rather than nothing; the collar solve is unaffected, +// since maximise_profit_over_collar handles a pinned optimum explicitly. // need docs on Golden Section Search. inline void Leaf::optimise_psi_stem_Sperry() { - // These two optimise psi_stem directly and produce no root-collar operating - // point, so they have no kind to report -- clear it rather than leave the - // last collar solve's classification standing over their outputs. - operating_point_kind_ = OperatingPointKind::Unsolved; + clear_collar_solve_state(); if (!supply_is_single_layer()) { util::stop("psi soil must have only one value to use non-root-based profit optimisation methods"); } + // ⚠️ lambda_ IS AN INPUT HERE AND HAS NO DEFAULT. setup_clean_leaf sets it to + // the NA sentinel and set_physiology does not touch it, so a caller who drives + // the leaf and calls this without setting it is optimising a NaN objective -- + // and brent_fmin does not report that. It returns a plausible-looking potential + // (2.551266 MPa at the package defaults) with profit_ = NaN beside it, and the + // potential is a property of the bracket rather than of the model. Refuse. + // + // Note the asymmetry this guards: set_traits() clears lambda_ (via + // setup_clean_leaf) and set_drivers() does not, so whether it survives depends + // on the order of two calls that look interchangeable. + if (!std::isfinite(lambda_)) { + util::stop("optimise_psi_stem_Sperry needs lambda_ set: it is a PRESCRIBED " + "marginal water cost, cleared to NA by the constructor and by " + "set_traits, and never set by set_physiology. For Sperry (2017)'s " + "own normalised profit use optimise_psi_stem_ProfitMax, which " + "computes its own."); + } + opt_psi_stem_ = supply_psi_soil_scalar(); @@ -3265,10 +3525,258 @@ inline void Leaf::optimise_psi_stem_Sperry() { } +// =========================================================================== +// Sperry et al. (2017) ProfitMax, as Sicangco et al. (2026) implement it +// --------------------------------------------------------------------------- +// WHY THIS EXISTS ALONGSIDE optimise_psi_stem_Sperry, WHICH IS THE SAME MODEL. +// Sperry maximises `Profit = CG - HC` with both terms normalised, where this +// package's older entry point maximises `A - lambda*(k(psi_soil)-k(psi))`. +// Multiplying Sperry's objective by |A|max shows the two are the same function up +// to a positive scale factor, so they share an argmax EXACTLY when lambda takes +// the value below: +// +// CG = A(psi)/|A|max +// HC = [k(psi_soil)-k(psi)] / [k(psi_soil)-kcrit] +// lambda* = |A|max / [k(psi_soil) - kcrit] +// +// Checked numerically on a 4001-point grid at gross assimilation, at net, and at +// a leaf temperature of 48 C where CG is negative throughout: same grid point +// every time. So the lambda form is not wrong -- it is unusable, because lambda* +// is not a constant. |A|max moves with every driver and k(psi_soil) with the +// soil, so a caller has to rescan the supply stream and recompute lambda at every +// observation, from R, at ~1.8 us of call overhead per crossing against a ~3 us +// model. This does it in C++ and reports the profit in Sperry's own units, so a +// figure from this package can be laid over one from the paper. +// +// THE NORMALISATION IS NOT COSMETIC, which is the other reason to have it. HC +// runs 0 to 1 across the operating range HOWEVER WIDE the vulnerability curve is, +// because the denominator rescales with the curve -- so a ProfitMax cost cannot +// become numerically negligible the way the TF24 cost does as stem_b widens +// (measured at 1.3e-5 at stem_b = 200). Both terms are scale-free, which also +// means the plant's willingness to spend water does not depend on how much carbon +// is at stake in absolute terms. That is exactly what cost_scale_TF24 decides on +// the TF24 path, and it is the parameter a calibration cannot identify without +// leaf water potential. +// +// ⚠️ kcrit IS 5% OF kmax, and that is the same convention this package already +// has rather than a new one. Sperry, Sabot et al. (2020) and Sicangco all set +// kcrit = 0.05*kmax, i.e. psi_crit is P95; at this package's defaults +// f(psi_crit) = 0.0500 exactly, so nothing has to be reconciled. +inline void Leaf::prepare_profitmax() { + if (!supply_is_single_layer()) { + util::stop("psi soil must have only one value to use non-root-based profit " + "optimisation methods"); + } + const double psi_soil = supply_psi_soil_scalar(); + profitmax_k_soil_ = + leaf_specific_conductance_max_ * proportion_of_conductivity(psi_soil); + profitmax_k_span_ = + profitmax_k_soil_ - + leaf_specific_conductance_max_ * proportion_of_conductivity(psi_crit); + + // |A|max over the transpiration supply stream (Sperry Eqn 4; Sicangco use + // max|Anet| because CG_net can be negative). A scan rather than "A at psi_crit" + // on purpose: that shortcut is only valid for GROSS assimilation, where A is + // monotone in psi, and the net-assimilation arms of this paper are precisely + // where it stops being. + // + // ⚠️ THE psi_soil ENDPOINT IS SKIPPED, matching gsthermal's `A[E == 0] <- NA`. + // There E is exactly zero and this model shuts down and reports A = -R_d, a + // number that describes a leaf with closed stomata rather than a point on the + // supply stream. Including it would set |A|max from respiration whenever + // assimilation is small, which is the whole high-temperature regime. + const int n = profitmax_scan_n_; + if (n < 3) { + util::stop("profitmax_scan_n_ must be at least 3"); + } + double a_max = 0.0; + const double step = (psi_crit - psi_soil) / double(n - 1); + profitmax_scan_psi_.assign(static_cast(n), util::na_value); + profitmax_scan_A_.assign(static_cast(n), util::na_value); + profitmax_scan_Tleaf_.assign(static_cast(n), util::na_value); + for (int i = 0; i < n; ++i) { + const double p = psi_soil + step * double(i); + set_leaf_states_rates_from_psi_stem(p, psi_soil); + const std::size_t k = static_cast(i); + profitmax_scan_psi_[k] = p; + profitmax_scan_A_[k] = assim_colimited_; + profitmax_scan_Tleaf_[k] = + use_energy_balance_ ? leaf_temp_from_E(transpiration_) : leaf_temp_; + // i == 0 is the psi_soil endpoint, where E is exactly zero and this model + // reports A = -R_d: a closed stoma rather than a point on the supply stream. + // It is RECORDED (the profit at full closure is a legitimate candidate) but + // excluded from |A|max, which is what gsthermal's `A[E == 0] <- NA` does. + if (i > 0 && std::isfinite(assim_colimited_)) { + a_max = std::max(a_max, std::abs(assim_colimited_)); + } + } + profitmax_A_max_ = a_max; +} + +// CG - (HC + TC) at one candidate potential, with the normalisers prepare_profitmax +// seeded. Writes carbon_gain_, hydraulic_cost_norm_ and thermal_cost_ so a caller +// plotting the paper's Figure 2 can read the three components off the object. +inline double Leaf::profit_psi_stem_ProfitMax(double psi_stem, + double psi_upstream) { + set_leaf_states_rates_from_psi_stem(psi_stem, psi_upstream); + + carbon_gain_ = (profitmax_A_max_ > 0.0) + ? assim_colimited_ / profitmax_A_max_ + : 0.0; + hydraulic_cost_norm_ = + (profitmax_k_soil_ - + leaf_specific_conductance_max_ * proportion_of_conductivity(psi_stem)) / + profitmax_k_span_; + // At the operating-point leaf temperature, which on the energy-balance path is + // the one set_leaf_states_rates_from_psi_stem just solved for. Off it the leaf + // is at the prescribed temperature and TC is constant across the stream -- + // which is Sicangco's ProfitMaxTC run with the energy balance disabled, and is + // why their thermal cost does nothing without it. + thermal_cost_ = thermal_cost_at(use_energy_balance_ + ? leaf_temp_from_E(transpiration_) + : leaf_temp_); + + return carbon_gain_ - (hydraulic_cost_norm_ + thermal_cost_); +} + +inline std::vector Leaf::profitmax_curve(int n) { + if (n < 2) { + util::stop("profitmax_curve needs at least 2 points"); + } + prepare_profitmax(); + const double psi_soil = supply_psi_soil_scalar(); + const double step = (psi_crit - psi_soil) / double(n - 1); + std::vector out(static_cast(5 * n)); + for (int i = 0; i < n; ++i) { + const double p = psi_soil + step * double(i); + const double profit = profit_psi_stem_ProfitMax(p, psi_soil); + const std::size_t k = static_cast(i); + const std::size_t N = static_cast(n); + out[k] = p; + out[N + k] = carbon_gain_; + out[2 * N + k] = hydraulic_cost_norm_; + out[3 * N + k] = thermal_cost_; + out[4 * N + k] = profit; + } + return out; +} + +inline void Leaf::optimise_psi_stem_ProfitMax() { + clear_collar_solve_state(); + + const double psi_soil = supply_psi_soil_scalar(); // also checks single-layer + opt_psi_stem_ = psi_soil; + + if ((PPFD_ < 1.5e-8) | (psi_soil > psi_crit)) { + profit_ = 0; + transpiration_ = 0; + stom_cond_CO2_ = 0; + carbon_gain_ = 0; + hydraulic_cost_norm_ = 0; + thermal_cost_ = 0; + lambda_ = util::na_value; + return; + } + + prepare_profitmax(); + if (!(profitmax_k_span_ > 0.0) || !(profitmax_A_max_ > 0.0)) { + // No usable normalisation: either the soil is already at the critical + // potential (no conductance to spend) or nothing on the stream assimilates. + profit_ = 0; + transpiration_ = 0; + stom_cond_CO2_ = 0; + carbon_gain_ = 0; + hydraulic_cost_norm_ = 0; + thermal_cost_ = 0; + lambda_ = util::na_value; + return; + } + + // Reported so the equivalence above is inspectable rather than asserted: this + // is the lambda that makes optimise_psi_stem_Sperry find the same point. + lambda_ = profitmax_A_max_ / profitmax_k_span_; + + // ⚠️ GRID FIRST, THEN REFINE, AND A BARE BRENT SEARCH IS WRONG HERE. + // + // This used to be `brent_fmin` over [psi_soil, psi_crit] alone. Brent is a + // LOCAL optimiser that steps in from the bounds, so it cannot return an + // endpoint and it cannot see past a local maximum -- and this objective has + // both of those, in exactly the regime the model is interesting in. Measured at + // Tair 50 C with the thermal cost on: the profit runs -1.5314 at psi_soil, + // -1.5510 at 1.19, -1.5459 at 1.88, then falls away, so the GLOBAL maximum is + // the closed-stomata endpoint and there is a local one near 1.9. Brent returned + // 1.643. The model was reporting a leaf with open stomata where the objective + // says it should be shut. + // + // It is not a hypothetical: full closure at high temperature is what Sicangco + // et al. (2026) report for their CGnet arms -- "for sufficiently high + // temperatures CGnet is negative for all possible values of Psi_leaf and the + // optimum shifts toward stomatal closure" -- and their own implementation finds + // it because it takes `which.max` over a 500-point grid rather than searching. + // + // So: evaluate the objective on the scan prepare_profitmax() has ALREADY run + // (no extra model evaluations -- A and Tleaf are stored, and HC and TC are + // analytic in psi and Tleaf), take the grid argmax, and refine with Brent only + // when that argmax is interior. An endpoint argmax is returned as the endpoint, + // which is the answer rather than a failure to search. + const std::size_t n = profitmax_scan_psi_.size(); + const double inv_A = 1.0 / profitmax_A_max_; + const double inv_k = 1.0 / profitmax_k_span_; + auto grid_profit = [&](std::size_t i) { + const double A = profitmax_scan_A_[i]; + if (!std::isfinite(A)) { + return -std::numeric_limits::infinity(); + } + const double hc = (profitmax_k_soil_ - leaf_specific_conductance_max_ * + proportion_of_conductivity( + profitmax_scan_psi_[i])) * + inv_k; + return A * inv_A - (hc + thermal_cost_at(profitmax_scan_Tleaf_[i])); + }; + + std::size_t best = 0; + double best_profit = grid_profit(0); + for (std::size_t i = 1; i < n; ++i) { + const double p = grid_profit(i); + if (p > best_profit) { + best_profit = p; + best = i; + } + } + + if (best == 0 || best + 1 == n) { + // Pinned to a bound. Re-evaluate through the real objective so every reported + // field describes the returned point rather than the grid's reconstruction. + opt_psi_stem_ = profitmax_scan_psi_[best]; + profit_ = profit_psi_stem_ProfitMax(opt_psi_stem_, psi_soil); + return; + } + + double neg_profit_opt = 0.0; + opt_psi_stem_ = util::brent_fmin( + [&](double psi_stem) { return -profit_psi_stem_ProfitMax(psi_stem, psi_soil); }, + profitmax_scan_psi_[best - 1], profitmax_scan_psi_[best + 1], GSS_tol_abs, + &neg_profit_opt); + profit_ = -neg_profit_opt; + + // ⚠️ Refining inside one grid cell can come out WORSE than the grid point when + // the cell is narrow relative to GSS_tol_abs, because Brent terminates on + // bracket width. Keep whichever is better; the grid point is always a feasible + // candidate. + if (best_profit > profit_) { + opt_psi_stem_ = profitmax_scan_psi_[best]; + } + + // brent_fmin's last evaluation is not necessarily at the returned argmax, so + // re-evaluate to leave every reported field describing ONE operating point. + // Hazard 8, in the form where the fields are individually plausible. + profit_ = profit_psi_stem_ProfitMax(opt_psi_stem_, psi_soil); +} + inline void Leaf::optimise_psi_stem_TF() { - // See optimise_psi_stem_Sperry: no root-collar operating point, so no kind. - operating_point_kind_ = OperatingPointKind::Unsolved; + // See optimise_psi_stem_Sperry: no root-collar operating point. + clear_collar_solve_state(); if (!supply_is_single_layer()) { util::stop("psi soil must have only one value to use non-root-based profit optimisation methods"); diff --git a/man/Leaf.Rd b/man/Leaf.Rd index 360a3ad..68c8f2a 100644 --- a/man/Leaf.Rd +++ b/man/Leaf.Rd @@ -73,6 +73,14 @@ Set the drivers with `$set_physiology()`, solve with `$opt_psi_stem_`, `$profit_`, `$assim_colimited_`, `$transpiration_`, `$lambda`, `$g1_eff`. +A second optimality formulation is available on the single-soil-layer +supply path: `$optimise_psi_stem_ProfitMax()` maximises Sperry (2017)'s +normalised profit, and `$profitmax_curve()` returns the whole cost, gain +and profit curve in one call. `$optimise_psi_stem_TF()` runs the TF24 cost +through the same solver, so the two can be compared at identical drivers. +These write `$carbon_gain_`, `$hydraulic_cost_norm_` and `$thermal_cost_`, +which are unitless and are NOT the same quantities as `$hydraulic_cost_`. + All water potentials are **positive magnitudes in MPa**. There is one representation throughout and it is asserted, not documented: a negative `psi_crit`, `stem_b`, `root_b` or `root_psi_crit` is an error here, and a diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 01c159e..a080330 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -281,6 +281,18 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// Leaf__leaf_temp_from_E +double Leaf__leaf_temp_from_E(phylloptim::RcppR6::RcppR6 obj_, double E); +RcppExport SEXP _phylloptim_Leaf__leaf_temp_from_E(SEXP obj_SEXP, SEXP ESEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type E(ESEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__leaf_temp_from_E(obj_, E)); + return rcpp_result_gen; +END_RCPP +} // Leaf__set_leaf_states_rates_from_psi_stem void Leaf__set_leaf_states_rates_from_psi_stem(phylloptim::RcppR6::RcppR6 obj_, double psi_stem, double psi_upstream); RcppExport SEXP _phylloptim_Leaf__set_leaf_states_rates_from_psi_stem(SEXP obj_SEXP, SEXP psi_stemSEXP, SEXP psi_upstreamSEXP) { @@ -464,6 +476,63 @@ BEGIN_RCPP return R_NilValue; END_RCPP } +// Leaf__thermal_cost_at +double Leaf__thermal_cost_at(phylloptim::RcppR6::RcppR6 obj_, double leaf_temp); +RcppExport SEXP _phylloptim_Leaf__thermal_cost_at(SEXP obj_SEXP, SEXP leaf_tempSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type leaf_temp(leaf_tempSEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__thermal_cost_at(obj_, leaf_temp)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__prepare_profitmax +void Leaf__prepare_profitmax(phylloptim::RcppR6::RcppR6 obj_); +RcppExport SEXP _phylloptim_Leaf__prepare_profitmax(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Leaf__prepare_profitmax(obj_); + return R_NilValue; +END_RCPP +} +// Leaf__profit_psi_stem_ProfitMax +double Leaf__profit_psi_stem_ProfitMax(phylloptim::RcppR6::RcppR6 obj_, double psi_stem, double psi_upstream); +RcppExport SEXP _phylloptim_Leaf__profit_psi_stem_ProfitMax(SEXP obj_SEXP, SEXP psi_stemSEXP, SEXP psi_upstreamSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type psi_stem(psi_stemSEXP); + Rcpp::traits::input_parameter< double >::type psi_upstream(psi_upstreamSEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__profit_psi_stem_ProfitMax(obj_, psi_stem, psi_upstream)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__profitmax_curve +std::vector Leaf__profitmax_curve(phylloptim::RcppR6::RcppR6 obj_, int n); +RcppExport SEXP _phylloptim_Leaf__profitmax_curve(SEXP obj_SEXP, SEXP nSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< int >::type n(nSEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__profitmax_curve(obj_, n)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__optimise_psi_stem_ProfitMax +void Leaf__optimise_psi_stem_ProfitMax(phylloptim::RcppR6::RcppR6 obj_); +RcppExport SEXP _phylloptim_Leaf__optimise_psi_stem_ProfitMax(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Leaf__optimise_psi_stem_ProfitMax(obj_); + return R_NilValue; +END_RCPP +} // Leaf__optimise_psi_stem_TF void Leaf__optimise_psi_stem_TF(phylloptim::RcppR6::RcppR6 obj_); RcppExport SEXP _phylloptim_Leaf__optimise_psi_stem_TF(SEXP obj_SEXP) { @@ -868,6 +937,72 @@ BEGIN_RCPP return R_NilValue; END_RCPP } +// Leaf__carbon_gain___get +double Leaf__carbon_gain___get(phylloptim::RcppR6::RcppR6 obj_); +RcppExport SEXP _phylloptim_Leaf__carbon_gain___get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__carbon_gain___get(obj_)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__carbon_gain___set +void Leaf__carbon_gain___set(phylloptim::RcppR6::RcppR6 obj_, double value); +RcppExport SEXP _phylloptim_Leaf__carbon_gain___set(SEXP obj_SEXP, SEXP valueSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type value(valueSEXP); + Leaf__carbon_gain___set(obj_, value); + return R_NilValue; +END_RCPP +} +// Leaf__hydraulic_cost_norm___get +double Leaf__hydraulic_cost_norm___get(phylloptim::RcppR6::RcppR6 obj_); +RcppExport SEXP _phylloptim_Leaf__hydraulic_cost_norm___get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__hydraulic_cost_norm___get(obj_)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__hydraulic_cost_norm___set +void Leaf__hydraulic_cost_norm___set(phylloptim::RcppR6::RcppR6 obj_, double value); +RcppExport SEXP _phylloptim_Leaf__hydraulic_cost_norm___set(SEXP obj_SEXP, SEXP valueSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type value(valueSEXP); + Leaf__hydraulic_cost_norm___set(obj_, value); + return R_NilValue; +END_RCPP +} +// Leaf__thermal_cost___get +double Leaf__thermal_cost___get(phylloptim::RcppR6::RcppR6 obj_); +RcppExport SEXP _phylloptim_Leaf__thermal_cost___get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__thermal_cost___get(obj_)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__thermal_cost___set +void Leaf__thermal_cost___set(phylloptim::RcppR6::RcppR6 obj_, double value); +RcppExport SEXP _phylloptim_Leaf__thermal_cost___set(SEXP obj_SEXP, SEXP valueSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type value(valueSEXP); + Leaf__thermal_cost___set(obj_, value); + return R_NilValue; +END_RCPP +} // Leaf__lambda_analytical___get double Leaf__lambda_analytical___get(phylloptim::RcppR6::RcppR6 obj_); RcppExport SEXP _phylloptim_Leaf__lambda_analytical___get(SEXP obj_SEXP) { @@ -2012,6 +2147,116 @@ BEGIN_RCPP return R_NilValue; END_RCPP } +// Leaf__vpd_leaf___get +double Leaf__vpd_leaf___get(phylloptim::RcppR6::RcppR6 obj_); +RcppExport SEXP _phylloptim_Leaf__vpd_leaf___get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__vpd_leaf___get(obj_)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__vpd_leaf___set +void Leaf__vpd_leaf___set(phylloptim::RcppR6::RcppR6 obj_, double value); +RcppExport SEXP _phylloptim_Leaf__vpd_leaf___set(SEXP obj_SEXP, SEXP valueSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type value(valueSEXP); + Leaf__vpd_leaf___set(obj_, value); + return R_NilValue; +END_RCPP +} +// Leaf__use_thermal_cost___get +bool Leaf__use_thermal_cost___get(phylloptim::RcppR6::RcppR6 obj_); +RcppExport SEXP _phylloptim_Leaf__use_thermal_cost___get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__use_thermal_cost___get(obj_)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__use_thermal_cost___set +void Leaf__use_thermal_cost___set(phylloptim::RcppR6::RcppR6 obj_, bool value); +RcppExport SEXP _phylloptim_Leaf__use_thermal_cost___set(SEXP obj_SEXP, SEXP valueSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< bool >::type value(valueSEXP); + Leaf__use_thermal_cost___set(obj_, value); + return R_NilValue; +END_RCPP +} +// Leaf__T50___get +double Leaf__T50___get(phylloptim::RcppR6::RcppR6 obj_); +RcppExport SEXP _phylloptim_Leaf__T50___get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__T50___get(obj_)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__T50___set +void Leaf__T50___set(phylloptim::RcppR6::RcppR6 obj_, double value); +RcppExport SEXP _phylloptim_Leaf__T50___set(SEXP obj_SEXP, SEXP valueSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type value(valueSEXP); + Leaf__T50___set(obj_, value); + return R_NilValue; +END_RCPP +} +// Leaf__Tcrit___get +double Leaf__Tcrit___get(phylloptim::RcppR6::RcppR6 obj_); +RcppExport SEXP _phylloptim_Leaf__Tcrit___get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__Tcrit___get(obj_)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__Tcrit___set +void Leaf__Tcrit___set(phylloptim::RcppR6::RcppR6 obj_, double value); +RcppExport SEXP _phylloptim_Leaf__Tcrit___set(SEXP obj_SEXP, SEXP valueSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< double >::type value(valueSEXP); + Leaf__Tcrit___set(obj_, value); + return R_NilValue; +END_RCPP +} +// Leaf__profitmax_scan_n___get +int Leaf__profitmax_scan_n___get(phylloptim::RcppR6::RcppR6 obj_); +RcppExport SEXP _phylloptim_Leaf__profitmax_scan_n___get(SEXP obj_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__profitmax_scan_n___get(obj_)); + return rcpp_result_gen; +END_RCPP +} +// Leaf__profitmax_scan_n___set +void Leaf__profitmax_scan_n___set(phylloptim::RcppR6::RcppR6 obj_, int value); +RcppExport SEXP _phylloptim_Leaf__profitmax_scan_n___set(SEXP obj_SEXP, SEXP valueSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< phylloptim::RcppR6::RcppR6 >::type obj_(obj_SEXP); + Rcpp::traits::input_parameter< int >::type value(valueSEXP); + Leaf__profitmax_scan_n___set(obj_, value); + return R_NilValue; +END_RCPP +} // Leaf__d___get double Leaf__d___get(phylloptim::RcppR6::RcppR6 obj_); RcppExport SEXP _phylloptim_Leaf__d___get(SEXP obj_SEXP) { @@ -2479,6 +2724,7 @@ static const R_CallMethodDef CallEntries[] = { {"_phylloptim_Leaf__assim_colimited", (DL_FUNC) &_phylloptim_Leaf__assim_colimited, 2}, {"_phylloptim_Leaf__assim_minus_stom_cond_CO2", (DL_FUNC) &_phylloptim_Leaf__assim_minus_stom_cond_CO2, 4}, {"_phylloptim_Leaf__electron_transport", (DL_FUNC) &_phylloptim_Leaf__electron_transport, 1}, + {"_phylloptim_Leaf__leaf_temp_from_E", (DL_FUNC) &_phylloptim_Leaf__leaf_temp_from_E, 2}, {"_phylloptim_Leaf__set_leaf_states_rates_from_psi_stem", (DL_FUNC) &_phylloptim_Leaf__set_leaf_states_rates_from_psi_stem, 3}, {"_phylloptim_Leaf__E_from_Soil_to_Root_Collar", (DL_FUNC) &_phylloptim_Leaf__E_from_Soil_to_Root_Collar, 3}, {"_phylloptim_Leaf__find_root_collar_psi", (DL_FUNC) &_phylloptim_Leaf__find_root_collar_psi, 1}, @@ -2494,6 +2740,11 @@ static const R_CallMethodDef CallEntries[] = { {"_phylloptim_Leaf__profit_psi_stem_TF", (DL_FUNC) &_phylloptim_Leaf__profit_psi_stem_TF, 3}, {"_phylloptim_Leaf__lambda_TF24", (DL_FUNC) &_phylloptim_Leaf__lambda_TF24, 2}, {"_phylloptim_Leaf__optimise_psi_stem_Sperry", (DL_FUNC) &_phylloptim_Leaf__optimise_psi_stem_Sperry, 1}, + {"_phylloptim_Leaf__thermal_cost_at", (DL_FUNC) &_phylloptim_Leaf__thermal_cost_at, 2}, + {"_phylloptim_Leaf__prepare_profitmax", (DL_FUNC) &_phylloptim_Leaf__prepare_profitmax, 1}, + {"_phylloptim_Leaf__profit_psi_stem_ProfitMax", (DL_FUNC) &_phylloptim_Leaf__profit_psi_stem_ProfitMax, 3}, + {"_phylloptim_Leaf__profitmax_curve", (DL_FUNC) &_phylloptim_Leaf__profitmax_curve, 2}, + {"_phylloptim_Leaf__optimise_psi_stem_ProfitMax", (DL_FUNC) &_phylloptim_Leaf__optimise_psi_stem_ProfitMax, 1}, {"_phylloptim_Leaf__optimise_psi_stem_TF", (DL_FUNC) &_phylloptim_Leaf__optimise_psi_stem_TF, 1}, {"_phylloptim_Leaf__medlyn_model_gs", (DL_FUNC) &_phylloptim_Leaf__medlyn_model_gs, 2}, {"_phylloptim_Leaf__solve_medlyn_ci_numerical", (DL_FUNC) &_phylloptim_Leaf__solve_medlyn_ci_numerical, 1}, @@ -2531,6 +2782,12 @@ static const R_CallMethodDef CallEntries[] = { {"_phylloptim_Leaf__psi_stem__set", (DL_FUNC) &_phylloptim_Leaf__psi_stem__set, 2}, {"_phylloptim_Leaf__lambda___get", (DL_FUNC) &_phylloptim_Leaf__lambda___get, 1}, {"_phylloptim_Leaf__lambda___set", (DL_FUNC) &_phylloptim_Leaf__lambda___set, 2}, + {"_phylloptim_Leaf__carbon_gain___get", (DL_FUNC) &_phylloptim_Leaf__carbon_gain___get, 1}, + {"_phylloptim_Leaf__carbon_gain___set", (DL_FUNC) &_phylloptim_Leaf__carbon_gain___set, 2}, + {"_phylloptim_Leaf__hydraulic_cost_norm___get", (DL_FUNC) &_phylloptim_Leaf__hydraulic_cost_norm___get, 1}, + {"_phylloptim_Leaf__hydraulic_cost_norm___set", (DL_FUNC) &_phylloptim_Leaf__hydraulic_cost_norm___set, 2}, + {"_phylloptim_Leaf__thermal_cost___get", (DL_FUNC) &_phylloptim_Leaf__thermal_cost___get, 1}, + {"_phylloptim_Leaf__thermal_cost___set", (DL_FUNC) &_phylloptim_Leaf__thermal_cost___set, 2}, {"_phylloptim_Leaf__lambda_analytical___get", (DL_FUNC) &_phylloptim_Leaf__lambda_analytical___get, 1}, {"_phylloptim_Leaf__lambda_analytical___set", (DL_FUNC) &_phylloptim_Leaf__lambda_analytical___set, 2}, {"_phylloptim_Leaf__hydraulic_cost___get", (DL_FUNC) &_phylloptim_Leaf__hydraulic_cost___get, 1}, @@ -2635,6 +2892,16 @@ static const R_CallMethodDef CallEntries[] = { {"_phylloptim_Leaf__f_r__set", (DL_FUNC) &_phylloptim_Leaf__f_r__set, 2}, {"_phylloptim_Leaf__use_energy_balance___get", (DL_FUNC) &_phylloptim_Leaf__use_energy_balance___get, 1}, {"_phylloptim_Leaf__use_energy_balance___set", (DL_FUNC) &_phylloptim_Leaf__use_energy_balance___set, 2}, + {"_phylloptim_Leaf__vpd_leaf___get", (DL_FUNC) &_phylloptim_Leaf__vpd_leaf___get, 1}, + {"_phylloptim_Leaf__vpd_leaf___set", (DL_FUNC) &_phylloptim_Leaf__vpd_leaf___set, 2}, + {"_phylloptim_Leaf__use_thermal_cost___get", (DL_FUNC) &_phylloptim_Leaf__use_thermal_cost___get, 1}, + {"_phylloptim_Leaf__use_thermal_cost___set", (DL_FUNC) &_phylloptim_Leaf__use_thermal_cost___set, 2}, + {"_phylloptim_Leaf__T50___get", (DL_FUNC) &_phylloptim_Leaf__T50___get, 1}, + {"_phylloptim_Leaf__T50___set", (DL_FUNC) &_phylloptim_Leaf__T50___set, 2}, + {"_phylloptim_Leaf__Tcrit___get", (DL_FUNC) &_phylloptim_Leaf__Tcrit___get, 1}, + {"_phylloptim_Leaf__Tcrit___set", (DL_FUNC) &_phylloptim_Leaf__Tcrit___set, 2}, + {"_phylloptim_Leaf__profitmax_scan_n___get", (DL_FUNC) &_phylloptim_Leaf__profitmax_scan_n___get, 1}, + {"_phylloptim_Leaf__profitmax_scan_n___set", (DL_FUNC) &_phylloptim_Leaf__profitmax_scan_n___set, 2}, {"_phylloptim_Leaf__d___get", (DL_FUNC) &_phylloptim_Leaf__d___get, 1}, {"_phylloptim_Leaf__d___set", (DL_FUNC) &_phylloptim_Leaf__d___set, 2}, {"_phylloptim_Leaf__wind_speed___get", (DL_FUNC) &_phylloptim_Leaf__wind_speed___get, 1}, diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index 6c365ec..a1168d4 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -80,6 +80,10 @@ double Leaf__electron_transport(phylloptim::RcppR6::RcppR6 obj return obj_->electron_transport(); } // [[Rcpp::export]] +double Leaf__leaf_temp_from_E(phylloptim::RcppR6::RcppR6 obj_, double E) { + return obj_->leaf_temp_from_E(E); +} +// [[Rcpp::export]] void Leaf__set_leaf_states_rates_from_psi_stem(phylloptim::RcppR6::RcppR6 obj_, double psi_stem, double psi_upstream) { obj_->set_leaf_states_rates_from_psi_stem(psi_stem, psi_upstream); } @@ -140,6 +144,26 @@ void Leaf__optimise_psi_stem_Sperry(phylloptim::RcppR6::RcppR6 obj_->optimise_psi_stem_Sperry(); } // [[Rcpp::export]] +double Leaf__thermal_cost_at(phylloptim::RcppR6::RcppR6 obj_, double leaf_temp) { + return obj_->thermal_cost_at(leaf_temp); +} +// [[Rcpp::export]] +void Leaf__prepare_profitmax(phylloptim::RcppR6::RcppR6 obj_) { + obj_->prepare_profitmax(); +} +// [[Rcpp::export]] +double Leaf__profit_psi_stem_ProfitMax(phylloptim::RcppR6::RcppR6 obj_, double psi_stem, double psi_upstream) { + return obj_->profit_psi_stem_ProfitMax(psi_stem, psi_upstream); +} +// [[Rcpp::export]] +std::vector Leaf__profitmax_curve(phylloptim::RcppR6::RcppR6 obj_, int n) { + return obj_->profitmax_curve(n); +} +// [[Rcpp::export]] +void Leaf__optimise_psi_stem_ProfitMax(phylloptim::RcppR6::RcppR6 obj_) { + obj_->optimise_psi_stem_ProfitMax(); +} +// [[Rcpp::export]] void Leaf__optimise_psi_stem_TF(phylloptim::RcppR6::RcppR6 obj_) { obj_->optimise_psi_stem_TF(); } @@ -308,6 +332,33 @@ void Leaf__lambda___set(phylloptim::RcppR6::RcppR6 obj_, doubl obj_->lambda_ = value; } +// [[Rcpp::export]] +double Leaf__carbon_gain___get(phylloptim::RcppR6::RcppR6 obj_) { + return obj_->carbon_gain_; +} +// [[Rcpp::export]] +void Leaf__carbon_gain___set(phylloptim::RcppR6::RcppR6 obj_, double value) { + obj_->carbon_gain_ = value; +} + +// [[Rcpp::export]] +double Leaf__hydraulic_cost_norm___get(phylloptim::RcppR6::RcppR6 obj_) { + return obj_->hydraulic_cost_norm_; +} +// [[Rcpp::export]] +void Leaf__hydraulic_cost_norm___set(phylloptim::RcppR6::RcppR6 obj_, double value) { + obj_->hydraulic_cost_norm_ = value; +} + +// [[Rcpp::export]] +double Leaf__thermal_cost___get(phylloptim::RcppR6::RcppR6 obj_) { + return obj_->thermal_cost_; +} +// [[Rcpp::export]] +void Leaf__thermal_cost___set(phylloptim::RcppR6::RcppR6 obj_, double value) { + obj_->thermal_cost_ = value; +} + // [[Rcpp::export]] double Leaf__lambda_analytical___get(phylloptim::RcppR6::RcppR6 obj_) { return obj_->lambda_analytical_; @@ -776,6 +827,51 @@ void Leaf__use_energy_balance___set(phylloptim::RcppR6::RcppR6 obj_->use_energy_balance_ = value; } +// [[Rcpp::export]] +double Leaf__vpd_leaf___get(phylloptim::RcppR6::RcppR6 obj_) { + return obj_->vpd_leaf_; +} +// [[Rcpp::export]] +void Leaf__vpd_leaf___set(phylloptim::RcppR6::RcppR6 obj_, double value) { + obj_->vpd_leaf_ = value; +} + +// [[Rcpp::export]] +bool Leaf__use_thermal_cost___get(phylloptim::RcppR6::RcppR6 obj_) { + return obj_->use_thermal_cost_; +} +// [[Rcpp::export]] +void Leaf__use_thermal_cost___set(phylloptim::RcppR6::RcppR6 obj_, bool value) { + obj_->use_thermal_cost_ = value; +} + +// [[Rcpp::export]] +double Leaf__T50___get(phylloptim::RcppR6::RcppR6 obj_) { + return obj_->T50_; +} +// [[Rcpp::export]] +void Leaf__T50___set(phylloptim::RcppR6::RcppR6 obj_, double value) { + obj_->T50_ = value; +} + +// [[Rcpp::export]] +double Leaf__Tcrit___get(phylloptim::RcppR6::RcppR6 obj_) { + return obj_->Tcrit_; +} +// [[Rcpp::export]] +void Leaf__Tcrit___set(phylloptim::RcppR6::RcppR6 obj_, double value) { + obj_->Tcrit_ = value; +} + +// [[Rcpp::export]] +int Leaf__profitmax_scan_n___get(phylloptim::RcppR6::RcppR6 obj_) { + return obj_->profitmax_scan_n_; +} +// [[Rcpp::export]] +void Leaf__profitmax_scan_n___set(phylloptim::RcppR6::RcppR6 obj_, int value) { + obj_->profitmax_scan_n_ = value; +} + // [[Rcpp::export]] double Leaf__d___get(phylloptim::RcppR6::RcppR6 obj_) { return obj_->d_; diff --git a/tests/cpp/test_leaf.cpp b/tests/cpp/test_leaf.cpp index 51441d3..5847a63 100644 --- a/tests/cpp/test_leaf.cpp +++ b/tests/cpp/test_leaf.cpp @@ -1674,13 +1674,36 @@ void test_energy_balance_gate_off_is_inert() { // the same drivers. Without the contrast the test could pass on the Arrhenius // optimum alone, which has nothing to do with energy balance. // -// ⚠️ What this can and cannot show. stom_cond_CO2 still divides by the -// PRESCRIBED air VPD however hot the leaf gets (PLAN 13.1), so only the -// A(Tleaf) channel is exercised here; the leaf-to-air VPD channel is absent. -// The sweep therefore holds VPD fixed, which is the protocol of the paper's -// Figure S2 rather than its Figure 2. +// ⚠️ THE SIGNATURE IS IN TRANSPIRATION, NOT IN CONDUCTANCE, AND THAT CHANGED HERE. +// This test used to assert that CONDUCTANCE rises where assimilation falls. It +// did, and the reason was a bug: stom_cond_CO2 divided by the prescribed AIR +// deficit however hot the leaf got, so gs was simply E rescaled by a constant and +// inherited E's shape exactly. With the leaf-to-air deficit wired in (PLAN 13.1) +// the two come apart, because +// +// gs = P*E / (1.6 * D_leaf), D_leaf = atm_vpd + esat(Tleaf) - esat(Tair) +// +// and over the decoupled window D_leaf grows FASTER than E does. Measured on this +// sweep, between the thermal optimum and E's own peak: E x1.034, D_leaf x1.089, +// so gs x0.950 -- falling throughout. The gs counts went from a positive number +// to gate-on 0, gate-off 1. +// +// That is not a loss of the mechanism, and the assertions below are rewritten to +// say which is which: +// +// * the ENERGY-BALANCE mechanism is intact and is what the gate buys. E keeps +// rising above the thermal optimum, and it does so over a wider window with +// the gate on (4 points) than off (1). +// * the CONDUCTANCE signature, which is what a gas-exchange dataset reports, +// does not survive the correction at these drivers. Anyone comparing this +// model to a measured dgs/dT needs to know that, and the third assertion pins +// the ratio that causes it rather than the count it produces. +// +// The sweep holds AIR VPD fixed, which is the protocol of the paper's Figure S2 +// rather than its Figure 2 -- and note that fixing the air deficit no longer +// fixes the deficit the leaf actually sees. void test_energy_balance_stomatal_decoupling() { - printf("decoupling: gs rises where A falls, and only with the gate on\n"); + printf("decoupling: E rises where A falls; gs does not, once D moves with Tleaf\n"); // ⚠️ THE FIXTURE DEFAULTS CANNOT SHOW THIS, and the reason is physical rather // than numerical. kmax = K_s*theta/h is 3.14e-5 at the defaults, and a leaf // that conductive simply cannot move enough water to cool itself: measured in @@ -1692,7 +1715,7 @@ void test_energy_balance_stomatal_decoupling() { // until a test passes -- at the defaults the honest answer is "no decoupling", // and that is what the probe found. auto sweep = [&](bool gate) { - std::vector Ta, A, gs; + std::vector Ta, A, gs, E, D; for (double t = 15.0; t <= 45.0; t += 1.0) { Drivers d; d.PPFD = 1200.0; d.leaf_temp = t; d.atm_vpd = 1.0; d.h = 1.5; @@ -1712,39 +1735,61 @@ void test_energy_balance_stomatal_decoupling() { l.find_root_collar_psi(); if (!std::isfinite(l.assim_colimited_) || !std::isfinite(l.stom_cond_CO2_)) continue; Ta.push_back(t); A.push_back(l.assim_colimited_); gs.push_back(l.stom_cond_CO2_); + E.push_back(l.transpiration_); D.push_back(l.vpd_leaf_); } - return std::make_tuple(Ta, A, gs); + return std::make_tuple(Ta, A, gs, E, D); }; - auto count_decoupled = [](const std::vector& Ta, - const std::vector& A, - const std::vector& gs) { - if (Ta.size() < 5) return 0; + // Points above the thermal optimum where the response variable RISES while + // assimilation falls. `y` is E for the mechanism and gs for the signature. + auto count_decoupled = [](const std::vector& A, + const std::vector& y) { + if (A.size() < 5) return 0; const size_t peak = std::distance(A.begin(), std::max_element(A.begin(), A.end())); int n = 0; - for (size_t i = peak + 2; i + 1 < Ta.size(); ++i) { - const double dA = (A[i + 1] - A[i - 1]); - const double dg = (gs[i + 1] - gs[i - 1]); - if (dA < 0.0 && dg > 0.0) ++n; + for (size_t i = peak + 2; i + 1 < A.size(); ++i) { + if ((A[i + 1] - A[i - 1]) < 0.0 && (y[i + 1] - y[i - 1]) > 0.0) ++n; } return n; }; const auto on = sweep(true); const auto off = sweep(false); - const int n_on = count_decoupled(std::get<0>(on), std::get<1>(on), std::get<2>(on)); - const int n_off = count_decoupled(std::get<0>(off), std::get<1>(off), std::get<2>(off)); - printf(" decoupled points: gate on %d, gate off %d\n", n_on, n_off); - const auto& A_on = std::get<1>(on); + const auto& gs_on = std::get<2>(on); + const auto& E_on = std::get<3>(on); + const auto& D_on = std::get<4>(on); + + const int nE_on = count_decoupled(A_on, E_on); + const int nE_off = count_decoupled(std::get<1>(off), std::get<3>(off)); + const int ng_on = count_decoupled(A_on, gs_on); + const int ng_off = count_decoupled(std::get<1>(off), std::get<2>(off)); + printf(" decoupled points -- in E: gate on %d, gate off %d;" + " in gs: gate on %d, gate off %d\n", nE_on, nE_off, ng_on, ng_off); + ok(A_on.size() > 5, "the decoupling sweep produced a curve"); const size_t peak = std::distance(A_on.begin(), std::max_element(A_on.begin(), A_on.end())); ok(peak > 0 && peak + 1 < A_on.size(), "assimilation has an interior thermal optimum"); - ok(n_on > 0, "with the gate ON, conductance rises where assimilation falls"); - ok(n_on > n_off, - "and it does so more than with the gate off -- the energy balance is why"); + ok(nE_on > 0, "with the gate ON, transpiration rises where assimilation falls"); + ok(nE_on > nE_off, + "and it does so over a wider window than with the gate off -- " + "the energy balance is why"); + + // WHY the conductance signature does not follow, stated as the ratio rather + // than as the count it happens to produce here. Over the window from the + // thermal optimum to E's own peak, the deficit outgrows the flux, and + // gs = P*E/(1.6*D) therefore falls even though E is rising. + const size_t Epeak = std::distance(E_on.begin(), + std::max_element(E_on.begin(), E_on.end())); + ok(Epeak > peak, "transpiration peaks ABOVE the assimilation optimum"); + const double E_growth = E_on[Epeak] / E_on[peak]; + const double D_growth = D_on[Epeak] / D_on[peak]; + printf(" optimum -> E peak: E x%.4f, D_leaf x%.4f, gs x%.4f\n", + E_growth, D_growth, gs_on[Epeak] / gs_on[peak]); + ok(D_growth > E_growth, + "the leaf-to-air deficit outgrows the flux, so gs falls where E rises"); } // The closed-form fast path (leaf/closed_form.hpp). Two things matter: that it is @@ -3028,6 +3073,345 @@ void test_out_of_domain_under_rescale() { "the rescale is named so the two domains are not confused"); } +// =========================================================================== +// Leaf-to-air VPD (PLAN 13.1, #7) +// --------------------------------------------------------------------------- +void test_leaf_to_air_vpd() { + printf("leaf-to-air VPD: the deficit Fick's law divides by\n"); + Drivers d; + d.atm_vpd = 1.5; + + // Gate off, the leaf IS at air temperature, so the two deficits must agree -- + // and EXACTLY, not nearly. That is what keeps the prescribed-temperature path, + // and therefore the golden file, untouched by this change. + phylloptim::Leaf off = make_pm_leaf(d, {2.0}, {1.0}, false); + ok(off.vpd_leaf_ == off.atm_vpd_, + "gate off: vpd_leaf_ is bit-identical to the driver"); + off.find_root_collar_psi(); + ok(off.vpd_leaf_ == off.atm_vpd_, + "gate off: and still is after a solve"); + + // Gate on, with the leaf running hot: the deficit it sees is larger than the + // air's, so the same water flux implies a SMALLER conductance. + phylloptim::Leaf on = make_pm_leaf(d, {2.0}, {1.0}, true); + on.Rn_ = 400.0; + on.ra_ = 12.0; + on.find_root_collar_psi(); + const double Tleaf = on.leaf_temp_from_E(on.transpiration_); + ok(Tleaf > on.Tair_, "the gate-on leaf is hotter than the air here"); + ok(on.vpd_leaf_ > on.atm_vpd_, + "gate on: a hotter leaf sees a larger deficit than the air's"); + // And the value is the definition, not an approximation of it. + near(on.vpd_leaf_, + on.atm_vpd_ + on.saturation_vapour_pressure(Tleaf) - + on.saturation_vapour_pressure(on.Tair_), + 1e-12, "vpd_leaf_ = atm_vpd + esat(Tleaf) - esat(Tair)"); + printf(" Tair %.1f C, Tleaf %.2f C: D_air %.3f kPa -> D_leaf %.3f kPa (x%.2f)\n", + on.Tair_, Tleaf, on.atm_vpd_, on.vpd_leaf_, on.vpd_leaf_ / on.atm_vpd_); + + // The consequence, stated as an identity rather than a direction: gs is E + // rescaled by the deficit, so getting the deficit wrong scales gs by the ratio. + near(on.stom_cond_CO2_, + on.atm_kpa_ * on.transpiration_ * phylloptim::kg_to_mol_h2o / + on.vpd_leaf_ / phylloptim::H2O_CO2_stom_diff_ratio, + 1e-12, "gs is the transpiration divided by the deficit it sees"); +} + +// Build a leaf on the single-potential path, which is what the ProfitMax entry +// points require and what Sicangco's model is (one soil potential, no root +// resistance network). +phylloptim::Leaf make_single_leaf(const Drivers &d, double psi_soil, + bool gate = false) { + phylloptim::Leaf l; + l.setup_transpiration(100); + l.setup_root_vulnerability(100); + l.use_energy_balance_ = gate; + l.set_supply_single(0.0); + phylloptim::RootNetwork rn; + rn.r_R_V_sum = std::vector{1.0e3}; + rn.r_R_H_min = std::vector{0.0}; + rn.r_R_V = std::vector{1.0e3}; + rn.c_r_V = std::vector{0.0}; + rn.c_r_H = std::vector{0.0}; + l.set_physiology(rn, d.PPFD, {psi_soil}, {1.0}, d.K_s * d.theta / d.h, + d.atm_vpd, d.ca, d.leaf_temp, d.atm_o2_kpa, d.atm_kpa); + return l; +} + +// =========================================================================== +// Sperry (2017) ProfitMax +// --------------------------------------------------------------------------- +// THE LOAD-BEARING CLAIM is that the normalised objective and this package's +// older `A - lambda*cost` form are the same function up to a positive scale, so +// they share an argmax when lambda = |A|max/(k_soil - kcrit). Everything the +// Sicangco et al. (2026) replication does rests on it. +void test_profitmax_matches_the_lambda_form() { + printf("ProfitMax: the normalised objective and the lambda form agree\n"); + Drivers d; + d.PPFD = 1500.0; + + for (double psi_soil : {0.5, 2.0}) { + for (double t : {25.0, 40.0}) { + d.leaf_temp = t; + phylloptim::Leaf l = make_single_leaf(d, psi_soil); + + l.optimise_psi_stem_ProfitMax(); + const double psi_pm = l.opt_psi_stem_; + const double lambda_star = l.lambda_; + ok(std::isfinite(lambda_star) && lambda_star > 0.0, + "ProfitMax reports the equivalent lambda"); + + l.lambda_ = lambda_star; + l.optimise_psi_stem_Sperry(); + + // Both are Brent searches terminating on bracket width GSS_tol_abs, so + // agreement to that scale is the most that can be asked of them. + near(psi_pm, l.opt_psi_stem_, 5.0e-3, + "the two objectives find the same collar potential"); + } + } +} + +void test_profitmax_normalisation() { + printf("ProfitMax: what the normalisation does and does not remove\n"); + Drivers d; + d.PPFD = 1500.0; + + phylloptim::Leaf l = make_single_leaf(d, 0.5); + const std::vector curve = l.profitmax_curve(101); + const std::size_t n = 101; + ok(curve.size() == 5 * n, "profitmax_curve returns five columns"); + + // HC runs from 0 at the soil potential to 1 at psi_crit, monotonically. That is + // the definition and it is what makes the cost unable to vanish. + near(curve[2 * n + 0], 0.0, 1e-12, "HC is zero at the soil potential"); + near(curve[2 * n + (n - 1)], 1.0, 1e-9, "HC is one at psi_crit"); + bool hc_monotone = true; + for (std::size_t i = 1; i < n; ++i) { + if (!(curve[2 * n + i] >= curve[2 * n + i - 1])) hc_monotone = false; + } + ok(hc_monotone, "HC increases monotonically along the supply stream"); + + // ⚠️ AND HC DOES NOT DEPEND ON kmax AT ALL. Both the numerator and the + // denominator carry one factor of leaf_specific_conductance_max_, so it cancels + // exactly. That is worth pinning because Sicangco's ProfitMaxkmax(T) arm gives + // kmax a temperature response and the paper describes it as changing the cost: + // it does not. It changes the SUPPLY, and reaches the cost only through which + // potentials the leaf can reach and how hot it gets there. + Drivers d2 = d; + d2.K_s = d.K_s * 3.0; + phylloptim::Leaf l2 = make_single_leaf(d2, 0.5); + const std::vector curve2 = l2.profitmax_curve(101); + double hc_worst = 0.0; + for (std::size_t i = 0; i < n; ++i) { + hc_worst = std::max(hc_worst, std::abs(curve2[2 * n + i] - curve[2 * n + i])); + } + ok(hc_worst < 1e-12, "HC is invariant to a 3x change in kmax"); + printf(" HC under kmax x3: worst difference %.3e\n", hc_worst); + + // The optimiser lands on the curve's own maximum. + l.optimise_psi_stem_ProfitMax(); + double best = -1e300; + std::size_t at = 0; + for (std::size_t i = 0; i < n; ++i) { + if (curve[4 * n + i] > best) { best = curve[4 * n + i]; at = i; } + } + ok(std::abs(l.opt_psi_stem_ - curve[at]) < 2.0 * (curve[1] - curve[0]), + "the optimiser lands within a grid step of the curve's maximum"); + ok(l.profit_ >= best - 1e-9, "and at no lower profit than the grid's best"); +} + +void test_profitmax_thermal_cost() { + printf("ProfitMax: the thermal cost, and that it is inert when off\n"); + Drivers d; + d.PPFD = 1500.0; + d.leaf_temp = 48.0; + + phylloptim::Leaf off = make_single_leaf(d, 0.5); + off.optimise_psi_stem_ProfitMax(); + ok(off.thermal_cost_ == 0.0, "gate off: TC is exactly zero"); + + phylloptim::Leaf on = make_single_leaf(d, 0.5); + on.use_thermal_cost_ = true; + on.T50_ = 50.4; + on.Tcrit_ = 46.5; + // The gate reaches jmax_ through the temperature block, so the drivers have to + // be re-supplied for the cache key to notice. Doing it the way a caller would. + on = make_single_leaf(d, 0.5); + on.use_thermal_cost_ = true; + on.T50_ = 50.4; + on.Tcrit_ = 46.5; + phylloptim::RootNetwork rn; + rn.r_R_V_sum = std::vector{1.0e3}; + rn.r_R_H_min = std::vector{0.0}; + rn.r_R_V = std::vector{1.0e3}; + rn.c_r_V = std::vector{0.0}; + rn.c_r_H = std::vector{0.0}; + on.set_physiology(rn, d.PPFD, {0.5}, {1.0}, d.K_s * d.theta / d.h, d.atm_vpd, + d.ca, d.leaf_temp, d.atm_o2_kpa, d.atm_kpa); + on.optimise_psi_stem_ProfitMax(); + ok(on.thermal_cost_ > 0.0 && on.thermal_cost_ < 1.0, + "gate on at 48 C: TC is in (0,1)"); + ok(on.jmax_ < off.jmax_, "and Jmax is scaled down by (1 - TC)"); + near(on.jmax_, off.jmax_ * (1.0 - on.thermal_cost_at(d.leaf_temp)), 1e-12, + "by exactly that factor"); + printf(" Tleaf %.1f C, Tcrit %.1f, T50 %.1f: TC %.4f, Jmax %.3f -> %.3f\n", + d.leaf_temp, on.Tcrit_, on.T50_, on.thermal_cost_, off.jmax_, on.jmax_); + + // ⚠️ THE COST AT Tcrit IS A FIXED 11.9%, AND THAT CONTRADICTS THE PAPER'S PROSE. + // Sicangco et al. write that "(1 - TC) equals one for temperatures below Tcrit, + // [so] Eqns 10 and 11 yield the same result under such conditions". The equation + // they cite does not do that: with r = 2/(T50 - Tcrit), the argument at Tcrit is + // exactly -2 whatever the two thresholds are, so + // + // TC(Tcrit) = 1/(1 + e^2) = 0.1192... + // + // independent of parameterisation. A leaf sitting AT its critical temperature + // has already lost 11.9% of Jmax in this model, and a leaf 4 K below it still + // pays 3%. Pinned as an identity because it is a property of the functional + // form rather than of the values in Table 2. + near(on.thermal_cost_at(on.Tcrit_), 1.0 / (1.0 + std::exp(2.0)), 1e-14, + "TC at Tcrit is 1/(1+e^2), not zero"); + { + phylloptim::Leaf wide = make_single_leaf(d, 0.5); + wide.use_thermal_cost_ = true; + wide.T50_ = 55.0; + wide.Tcrit_ = 43.0; + near(wide.thermal_cost_at(wide.Tcrit_), on.thermal_cost_at(on.Tcrit_), 1e-14, + "and is the same 11.9% for a threshold pair three times as wide"); + } + near(on.thermal_cost_at(on.T50_), 0.5, 1e-12, "and is exactly 0.5 at T50"); +} + +void test_single_layer_optimisers_clear_collar_state() { + printf("single-layer optimisers do not inherit a collar solve's outputs\n"); + Drivers d; + phylloptim::Leaf l = make_single_leaf(d, 0.5); + l.find_root_collar_psi(); + ok(std::isfinite(l.opt_root_psi_) && std::isfinite(l.E_up_), + "the collar solve wrote a collar operating point"); + + l.optimise_psi_stem_ProfitMax(); + ok(!std::isfinite(l.opt_root_psi_), "ProfitMax clears opt_root_psi_"); + ok(!std::isfinite(l.E_up_), "ProfitMax clears E_up_"); + bool consumption_cleared = true; + for (double c : l.soil_consumption_) { + if (std::isfinite(c)) consumption_cleared = false; + } + ok(consumption_cleared, "ProfitMax clears soil_consumption_"); + ok(std::isfinite(l.transpiration_) && l.transpiration_ > 0.0, + "while still writing its own transpiration"); +} + +void test_sperry_refuses_an_unset_lambda() { + printf("optimise_psi_stem_Sperry refuses an unset lambda\n"); + Drivers d; + phylloptim::Leaf l = make_single_leaf(d, 0.5); + ok(!std::isfinite(l.lambda_), "lambda_ starts unset"); + bool threw = false; + try { + l.optimise_psi_stem_Sperry(); + } catch (const std::exception &e) { + threw = true; + ok(std::string(e.what()).find("lambda_") != std::string::npos, + "and the message names lambda_"); + } + ok(threw, "rather than searching a NaN objective"); +} + +// ⚠️ THE TEST THE GOLDEN FILE CANNOT BE. `set_leaf_states_rates_from_psi_stem` +// used to zero transpiration wherever `assim_max_ < 0`, and the golden grid's +// minimum assim_max_ is 3.71, so it never reached the branch. The collar solve +// cannot reach it either -- prepare_collar_solve exits first -- so the ONLY way +// to see this is to call the forward evaluation directly in that regime, which is +// what the single-layer optimisers and the ProfitMax curve do. +void test_transpiration_survives_negative_assim() { + printf("water moves whether or not there is carbon to be had\n"); + Drivers d; + d.PPFD = 1500.0; + d.leaf_temp = 50.0; // hot enough that A(ci = ca) cannot cover R_d + + phylloptim::Leaf l = make_single_leaf(d, 0.5); + ok(l.assim_max_ < 0.0, "the regime is reached: assim_max_ is negative"); + + const double psi = 3.0; + l.set_leaf_states_rates_from_psi_stem(psi, 0.5); + const double E = l.transpiration_; + ok(E > 0.0, "transpiration follows the hydraulic supply, not the carbon"); + near(E, l.transpiration(psi, 0.5), 1e-12, + "and equals the supply function exactly"); + ok(l.stom_cond_CO2_ > 0.0, "so the conductance is positive too"); + + // The carbon state is the one the branch used to set by hand, and it now comes + // from the ci solver's own compensation-point fallback. + near(l.ci_, l.gamma_ * l.umol_per_mol_to_Pa_, 1e-9, + "ci sits at the compensation point"); + near(l.assim_colimited_, -l.R_d_, 1e-9, + "and net assimilation is exactly -R_d"); + printf(" Tleaf %.0f C: assim_max_ %.3f, E %.3e kg m-2 s-1, A %.3f\n", + d.leaf_temp, l.assim_max_, E, l.assim_colimited_); + + // ⚠️ AND TWO LEAVES AT THE SAME OPERATING POINT NOW AGREE ABOUT THE WATER. This + // is how the old behaviour was found: an arm optimised with respiration off and + // scored with it on reported a potential that moves water beside a transpiration + // of exactly zero. + Drivers dg = d; + phylloptim::Leaf gross = make_single_leaf(dg, 0.5); + gross.set_traits(96, 2.680147, 3.898245, 5.870283, 2.680147, 3.898245, + 5.870283, 1.5, 157.44, 0.30, 0.7, 0.99, 7.5, /*R_d_25=*/0.0); + phylloptim::RootNetwork rn; + rn.r_R_V_sum = std::vector{1.0e3}; + rn.r_R_H_min = std::vector{0.0}; + rn.r_R_V = std::vector{1.0e3}; + rn.c_r_V = std::vector{0.0}; + rn.c_r_H = std::vector{0.0}; + gross.set_physiology(rn, dg.PPFD, {0.5}, {1.0}, dg.K_s * dg.theta / dg.h, + dg.atm_vpd, dg.ca, dg.leaf_temp, dg.atm_o2_kpa, + dg.atm_kpa); + gross.set_leaf_states_rates_from_psi_stem(psi, 0.5); + ok(gross.assim_max_ > 0.0, "with R_d_25 = 0 the same drivers are NOT shut down"); + near(gross.transpiration_, E, 1e-12, + "and both leaves report the same transpiration at the same potential"); +} + +// ⚠️ THE OBJECTIVE IS NOT UNIMODAL AND ITS MAXIMUM CAN BE AN ENDPOINT. This is +// the test that a bare Brent search fails: at a leaf hot enough that net +// assimilation is negative everywhere, the profit is highest at FULL CLOSURE and +// there is a local maximum out in the interior. Brent steps in from the bounds +// and cannot return an endpoint, so it used to report the local one -- an open +// stoma where the model says the leaf should be shut. +void test_profitmax_finds_a_closed_optimum() { + printf("ProfitMax finds a boundary optimum, which Brent alone cannot\n"); + Drivers d; + d.PPFD = 1500.0; + d.leaf_temp = 50.0; + + phylloptim::Leaf l = make_single_leaf(d, 0.5); + l.use_thermal_cost_ = true; + l.optimise_psi_stem_ProfitMax(); + + // Reconstruct the objective on a coarse grid and find its global maximum + // independently of the solver. + const std::vector curve = l.profitmax_curve(201); + const std::size_t n = 201; + std::size_t best = 0; + for (std::size_t i = 1; i < n; ++i) { + if (curve[4 * n + i] > curve[4 * n + best]) best = i; + } + printf(" grid argmax at psi = %.4f (index %zu of %zu), solver returned %.4f\n", + curve[best], best, n, l.opt_psi_stem_); + + // profitmax_curve re-prepares, so re-solve before reading the operating point. + l.optimise_psi_stem_ProfitMax(); + const double step = curve[1] - curve[0]; + ok(std::abs(l.opt_psi_stem_ - curve[best]) < 3.0 * step, + "the solver lands on the objective's GLOBAL maximum, not a local one"); + + // And the profit it reports is at least the grid's best. + ok(l.profit_ >= curve[4 * n + best] - 1e-9, + "at no lower profit than the grid's best"); +} + void benchmark() { printf("\ntiming\n"); Drivers d; @@ -3105,6 +3489,14 @@ int main() { test_infeasible_is_a_distinct_failure(); test_out_of_domain_names_the_spline(); test_out_of_domain_under_rescale(); + test_leaf_to_air_vpd(); + test_profitmax_matches_the_lambda_form(); + test_profitmax_normalisation(); + test_profitmax_thermal_cost(); + test_single_layer_optimisers_clear_collar_state(); + test_sperry_refuses_an_unset_lambda(); + test_transpiration_survives_negative_assim(); + test_profitmax_finds_a_closed_optimum(); benchmark(); printf("\n%d checks, %d failures\n", checks, failures); diff --git a/tools/doxygen_filter.awk b/tools/doxygen_filter.awk index 1b5abd0..ea63f17 100644 --- a/tools/doxygen_filter.awk +++ b/tools/doxygen_filter.awk @@ -35,6 +35,25 @@ # one unreadable line. A run of lines indented relative to the surrounding # prose is reproduced exactly instead. # +# ⚠️ ONE RUN PER COMMENT BLOCK. Doxygen 1.9 -- which is what CI installs -- +# handles the FIRST `\verbatim` in a block and then drops the second one's +# OPEN, reporting `unexpected command endverbatim` at its close. The line it +# names is in the filtered stream and lands in unrelated code, so it reads +# like a parser bug somewhere else entirely: the report that found this was +# 120 lines past the comment responsible. Doxygen 1.17 renders the same input +# in silence, so a local `doxygen` run says nothing. +# +# Established by a probe header of eight isolated constructs rendered in CI: +# one run, banner rules, `|` in prose, `|` inside the run, and an emoji are +# all clean; every block with two runs errors, and only at the SECOND close. +# So write one indented display per comment block -- put the equations +# together rather than one on each side of a paragraph. `docs.yml` asserts +# it, because nothing else can. +# +# The `/*! \file */` block rule 2 emits is EXEMPT, and that is measured +# rather than assumed: closed_form.hpp's file block has carried two runs +# across a long green master. So the assertion counts `///` blocks only. +# # A run only OPENS after a blank line and only on a line that is not a list # item. Both conditions are load-bearing. Without the first, the hanging # indent under a bullet (` * POSITIVE magnitudes -- ...` followed by