From 05f0eaee732a0f52776638c30744c6f77d02bb80 Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Wed, 19 Aug 2026 12:21:31 +1000 Subject: [PATCH 1/9] Wire in leaf-to-air VPD, and Sperry's ProfitMax as he defines it Three changes, all on the energy-balance and single-layer paths; the golden file is bit-identical and the R gradient golden is bit-exact. 1. Fick's law now divides by the LEAF-to-air deficit, not the air's (PLAN 13.1, #7). Off the energy balance the leaf is at air temperature and `vpd_leaf_ == atm_vpd_` exactly, so nothing on the prescribed-temperature path moves. On it, the deficit is 1.5x the air's at Tair 25 and 3-4x by Tair 45, so gs and A were overstated by that factor. The collar first-order condition gains the matching `dgc_dT` term -- and the placeholder algebra that stood there while `dgc_dT` was zero did NOT generalise, so it is re-derived rather than sign-flipped. 2. `optimise_psi_stem_ProfitMax()` runs Sperry (2017) with both terms normalised, computing |A|max over the supply stream itself. It reports the lambda that makes `optimise_psi_stem_Sperry` find the same point, so the equivalence is checkable. `profitmax_curve()` returns the whole cost/gain/profit curve in one crossing. 3. An optional instantaneous thermal cost (Sicangco et al. 2026), default off: TC sigmoid in leaf temperature, Jmax scaled by (1-TC). NOT plant's ATLS ratchet, which has memory. Plus three bugs on the single-layer optimisers: a NaN `lambda_` was searched silently, and `opt_root_psi_`/`E_up_`/`soil_consumption_` survived from an earlier collar solve. The decoupling acceptance test now asserts on transpiration rather than conductance, because the deficit correction removes the conductance signature at those drivers. Measured in the test. Co-Authored-By: Claude Opus 5 (1M context) --- COMPARISON.md | 4 +- PLAN.md | 23 +- R/RcppExports.R | 84 +++++ R/RcppR6.R | 73 +++- README.md | 11 +- inst/RcppR6_classes.yml | 55 +++ inst/include/phylloptim/constants.hpp | 9 + inst/include/phylloptim/leaf_model.hpp | 456 +++++++++++++++++++++++-- src/RcppExports.cpp | 254 ++++++++++++++ src/RcppR6.cpp | 92 +++++ tests/cpp/test_leaf.cpp | 343 +++++++++++++++++-- 11 files changed, 1333 insertions(+), 71 deletions(-) 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/PLAN.md b/PLAN.md index 6cb9f4b..f6a4bb9 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 bffcb1a..0e7440c 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -133,6 +133,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_)) } @@ -229,6 +249,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_) } @@ -637,6 +681,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 f798ac6..d010780 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: abbaf8a5425a985fb81802acd9b06f64 +## Hash: 31cdcc80e0dd69f79183778f0325ac1a ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -183,6 +183,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) }, @@ -264,6 +279,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) @@ -621,6 +657,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 33947ed..05455a1 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 4af1b02..b65a47d 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -133,6 +133,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} @@ -222,6 +230,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 +438,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 d99448c..3ce2dab 100644 --- a/inst/include/phylloptim/constants.hpp +++ b/inst/include/phylloptim/constants.hpp @@ -106,6 +106,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 8c7eae0..d916eb0 100644 --- a/inst/include/phylloptim/leaf_model.hpp +++ b/inst/include/phylloptim/leaf_model.hpp @@ -305,8 +305,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 --------------------------------------- @@ -374,6 +408,22 @@ 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 + 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_; @@ -434,7 +484,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; @@ -739,9 +789,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 @@ -789,6 +839,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. @@ -982,9 +1045,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 ------------------- // @@ -1258,6 +1346,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 @@ -1278,7 +1372,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; @@ -1403,6 +1498,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 @@ -2148,6 +2248,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); @@ -2232,15 +2333,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; @@ -2265,7 +2366,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); @@ -2298,7 +2399,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); } @@ -2342,14 +2443,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; } @@ -2401,10 +2519,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 { @@ -2450,13 +2577,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); @@ -2487,6 +2624,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. @@ -2731,7 +2901,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; } @@ -2938,10 +3108,13 @@ inline void Leaf::set_leaf_states_rates_from_psi_stem(double psi_stem, double ps // is a single forward pass: recompute the Farquhar temperature params at // this candidate's Tleaf before solving for ci. Defeats the photo_temp // cache by design -- Tleaf varies per operating point. - update_temperature_dependent_params(leaf_temp_from_E(transpiration_)); + const double Tleaf_candidate = leaf_temp_from_E(transpiration_); + update_temperature_dependent_params(Tleaf_candidate); + // ...and the deficit Fick's law divides by, which moves with it. + set_leaf_vpd(Tleaf_candidate); } 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_); @@ -3009,7 +3182,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 @@ -3026,6 +3199,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) { @@ -3051,19 +3242,42 @@ 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); +} // 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(); @@ -3086,10 +3300,180 @@ 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, CG = A(psi)/|A|max, HC = [k(psi_soil)-k(psi)] / [k(psi_soil)-kcrit] +// +// while 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* = |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); + for (int i = 1; i < n; ++i) { + set_leaf_states_rates_from_psi_stem(psi_soil + step * double(i), psi_soil); + if (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_; + + double neg_profit_opt = 0.0; + opt_psi_stem_ = util::brent_fmin( + [&](double psi_stem) { return -profit_psi_stem_ProfitMax(psi_stem, psi_soil); }, + psi_soil, psi_crit, GSS_tol_abs, &neg_profit_opt); + profit_ = -neg_profit_opt; + + // 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_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/src/RcppExports.cpp b/src/RcppExports.cpp index a086b7e..80c6d79 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -452,6 +452,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) { @@ -713,6 +770,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) { @@ -1835,6 +1958,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) { @@ -2271,6 +2504,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}, @@ -2295,6 +2533,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}, @@ -2397,6 +2641,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 f83f997..6547070 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -136,6 +136,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(); } @@ -239,6 +259,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_; @@ -698,6 +745,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 e0c0dd7..bbfbc52 100644 --- a/tests/cpp/test_leaf.cpp +++ b/tests/cpp/test_leaf.cpp @@ -1593,13 +1593,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 @@ -1611,7 +1634,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; @@ -1631,39 +1654,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 @@ -2556,6 +2601,252 @@ 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"); +} + void benchmark() { printf("\ntiming\n"); Drivers d; @@ -2626,6 +2917,12 @@ int main() { test_bad_input_throws(); 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(); benchmark(); printf("\n%d checks, %d failures\n", checks, failures); From 56450b486aff0aa5ffbff8d205d1551f5cf72203 Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Wed, 19 Aug 2026 12:32:02 +1000 Subject: [PATCH 2/9] Water moves whether or not there is carbon to be had `set_leaf_states_rates_from_psi_stem` zeroed transpiration and conductance wherever `assim_max_ < 0`. Transpiration on that path is the hydraulic supply at the candidate potential 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: where assim_max_ < 0 the ci root-find has no root in [gamma*, ca] and takes its compensation-point fallback, giving ci = gamma* and net A = -R_d, which is what the branch set by hand. That fallback did not exist when the branch was written. Unreachable from the collar solve -- prepare_collar_solve has its own `assim_max_ < 0` exit and returns before any candidate is evaluated -- so plant is unaffected. The golden grid's minimum assim_max_ is 3.71, so the file is bit-identical and says nothing about this; test_transpiration_survives_negative_assim is the test that can. Co-Authored-By: Claude Opus 5 (1M context) --- R/RcppExports.R | 4 ++ R/RcppR6.R | 5 ++- inst/RcppR6_classes.yml | 8 ++++ inst/include/phylloptim/leaf_model.hpp | 34 +++++++++++++--- src/RcppExports.cpp | 13 ++++++ src/RcppR6.cpp | 4 ++ tests/cpp/test_leaf.cpp | 56 ++++++++++++++++++++++++++ 7 files changed, 118 insertions(+), 6 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 0e7440c..0c8386f 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)) } diff --git a/R/RcppR6.R b/R/RcppR6.R index d010780..80d3f2b 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: 31cdcc80e0dd69f79183778f0325ac1a +## Hash: 8c0706020c281dd3208a35ae9f2967de ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -141,6 +141,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) }, diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index b65a47d..2370d29 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -383,6 +383,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] diff --git a/inst/include/phylloptim/leaf_model.hpp b/inst/include/phylloptim/leaf_model.hpp index d916eb0..ff21fba 100644 --- a/inst/include/phylloptim/leaf_model.hpp +++ b/inst/include/phylloptim/leaf_model.hpp @@ -3087,16 +3087,40 @@ 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) { + // ⚠️ 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. if (psi_upstream >= psi_stem){ ci_ = gamma_*umol_per_mol_to_Pa_; transpiration_ = 0; stom_cond_CO2_ = 0; } else{ - if(assim_max_ < 0){ - ci_ = gamma_*umol_per_mol_to_Pa_; - transpiration_ = 0; - stom_cond_CO2_ = 0; - } 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 -> diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 80c6d79..6ccd2d2 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) { @@ -2490,6 +2502,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}, diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index 6547070..30fb70a 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); } diff --git a/tests/cpp/test_leaf.cpp b/tests/cpp/test_leaf.cpp index bbfbc52..9f1dd2f 100644 --- a/tests/cpp/test_leaf.cpp +++ b/tests/cpp/test_leaf.cpp @@ -2847,6 +2847,61 @@ void test_sperry_refuses_an_unset_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"); +} + void benchmark() { printf("\ntiming\n"); Drivers d; @@ -2923,6 +2978,7 @@ int main() { test_profitmax_thermal_cost(); test_single_layer_optimisers_clear_collar_state(); test_sperry_refuses_an_unset_lambda(); + test_transpiration_survives_negative_assim(); benchmark(); printf("\n%d checks, %d failures\n", checks, failures); From 38ed0baa01eedc1f03252617113ef42712ea5d9e Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Wed, 19 Aug 2026 12:50:48 +1000 Subject: [PATCH 3/9] Document the ProfitMax surface in the Leaf roxygen Regenerating man/ also picks up a stale gradient_par_names.Rd, which said fifteen names where its own roxygen source has said sixteen since R_d_25 was added. Co-Authored-By: Claude Opus 5 (1M context) --- R/RcppR6.R | 10 +++++++++- inst/RcppR6_classes.yml | 8 ++++++++ man/Leaf.Rd | 8 ++++++++ man/gradient_par_names.Rd | 4 ++-- 4 files changed, 27 insertions(+), 3 deletions(-) diff --git a/R/RcppR6.R b/R/RcppR6.R index 80d3f2b..936efc9 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: 8c0706020c281dd3208a35ae9f2967de +## Hash: 7a170c78f1f607f16c6e5c8e564b4d1d ##' @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 diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index 2370d29..06f35b3 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 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/man/gradient_par_names.Rd b/man/gradient_par_names.Rd index 058d15d..7a1d697 100644 --- a/man/gradient_par_names.Rd +++ b/man/gradient_par_names.Rd @@ -7,10 +7,10 @@ gradient_par_names() } \value{ -A character vector of fifteen names. +A character vector of sixteen names. } \description{ -The thirteen [leaf_traits()] in `set_traits()`'s argument order, then +The fourteen [leaf_traits()] in `set_traits()`'s argument order, then `leaf_specific_conductance_max` and `resistance`. Exported so that R's own copy of this order can be compared against it in a test: R passes integer positions into this enumeration, so appending to it is safe and reordering it From cf0587b938b8e6cdd57038eccc6b3eec82b96028 Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Wed, 19 Aug 2026 13:13:02 +1000 Subject: [PATCH 4/9] ProfitMax: scan the grid, then refine brent_fmin is a local optimiser that steps in from the bounds, so it could return neither an endpoint nor the global maximum of a multi-modal objective -- and this objective is both, in the regime the model is interesting in. Measured at Tair 50 C with the thermal cost on: the profit is highest at full closure (-1.5314 at psi_soil) with a local maximum near 1.9, and the solver returned 1.643. It reported an open stoma where the objective says shut. The scan prepare_profitmax() already runs now stores A and Tleaf per point, so the objective is rebuilt on it for no extra model evaluations. An endpoint argmax is returned as the endpoint. optimise_psi_stem_TF and _Sperry share the hazard and are documented rather than changed: neither has a scan to reuse. The collar solve is unaffected -- it handles a pinned optimum explicitly. Co-Authored-By: Claude Opus 5 (1M context) --- inst/include/phylloptim/leaf_model.hpp | 105 +++++++++++++++++++++++-- tests/cpp/test_leaf.cpp | 39 +++++++++ 2 files changed, 139 insertions(+), 5 deletions(-) diff --git a/inst/include/phylloptim/leaf_model.hpp b/inst/include/phylloptim/leaf_model.hpp index ff21fba..852caf3 100644 --- a/inst/include/phylloptim/leaf_model.hpp +++ b/inst/include/phylloptim/leaf_model.hpp @@ -420,6 +420,12 @@ class Leaf { 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_ @@ -3275,6 +3281,17 @@ inline void Leaf::clear_collar_solve_state() { 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() { @@ -3390,9 +3407,22 @@ inline void Leaf::prepare_profitmax() { } double a_max = 0.0; const double step = (psi_crit - psi_soil) / double(n - 1); - for (int i = 1; i < n; ++i) { - set_leaf_states_rates_from_psi_stem(psi_soil + step * double(i), psi_soil); - if (std::isfinite(assim_colimited_)) { + 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_)); } } @@ -3482,16 +3512,81 @@ inline void Leaf::optimise_psi_stem_ProfitMax() { // 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); }, - psi_soil, psi_crit, GSS_tol_abs, &neg_profit_opt); + 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_psi_stem_ProfitMax(opt_psi_stem_, psi_soil); + profit_ = profit_psi_stem_ProfitMax(opt_psi_stem_, psi_soil); } inline void Leaf::optimise_psi_stem_TF() { diff --git a/tests/cpp/test_leaf.cpp b/tests/cpp/test_leaf.cpp index 9f1dd2f..8746610 100644 --- a/tests/cpp/test_leaf.cpp +++ b/tests/cpp/test_leaf.cpp @@ -2902,6 +2902,44 @@ void test_transpiration_survives_negative_assim() { "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; @@ -2979,6 +3017,7 @@ int main() { 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); From 702ee9cc86f8fcf5e92761af8dc8973fdd0903de Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Wed, 19 Aug 2026 13:13:29 +1000 Subject: [PATCH 5/9] Record the boundary-optimum hazard in the guide Co-Authored-By: Claude Opus 5 (1M context) --- .claude/CLAUDE.md | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index 37f3843..c08d434 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -735,6 +735,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 `tests/validate/` holds the harnesses. Read the header of From 478200b54236f8f3a12f824fbbdff162c72fc331 Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Wed, 19 Aug 2026 15:41:46 +1000 Subject: [PATCH 6/9] Escape `|` in the doxygen filter MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Doxygen reads a `|` in a paragraph as a table delimiter, and `\verbatim` is not allowed inside a table cell: it drops the OPEN and then reports `unexpected command endverbatim` at a line that is neither -- here, 80 lines further on, inside unrelated code. `|` in prose is common in these headers and was harmless until a comment that had one also gained an indented block, which the filter turns into a verbatim. Found by comparing every verbatim-carrying comment block on this branch against every one on master: a pipe was the only structural feature present in mine and absent from all of theirs. ⚠️ Doxygen 1.17 renders the same input in silence; only the 1.9 that CI installs objects, so a local `doxygen` run is not evidence either way. Co-Authored-By: Claude Opus 5 (1M context) --- tools/doxygen_filter.awk | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/tools/doxygen_filter.awk b/tools/doxygen_filter.awk index 1b5abd0..6352901 100644 --- a/tools/doxygen_filter.awk +++ b/tools/doxygen_filter.awk @@ -45,11 +45,21 @@ # preformatted text and stop rendering as lists. # # ESCAPING. Text outside a verbatim run is escaped, because none of it was -# written with Doxygen in mind: `\`, `@`, `#`, `%`, `&`, `<` and `>` all mean -# something to Doxygen and here they never do. Without it, `#include ` +# written with Doxygen in mind: `\`, `@`, `#`, `%`, `&`, `<`, `>` and `|` all +# mean something to Doxygen and here they never do. Without it, `#include ` # becomes a broken link to an entity called "include" followed by a swallowed # HTML tag, and the `\int` in the roots.hpp head-loss note becomes an unknown -# command; both were observed before this was added. Verbatim runs are NOT +# command; both were observed before this was added. +# +# ⚠️ `|` WAS ADDED LAST AND COST A CI-ONLY DEBUGGING ROUND. Doxygen reads a `|` +# in a paragraph as a table delimiter, and `\verbatim` is not allowed inside a +# table cell: it drops the OPEN, and then reports `unexpected command +# endverbatim` at a line that is neither -- in the case that found this, 80 lines +# further on, inside unrelated code. Doxygen 1.17 renders the same input in +# silence, so a local `doxygen` run says nothing; only the 1.9 that CI installs +# objects. The comment that triggered it wrote the maximum of A over the supply +# stream as `|A|max`. `|` in prose is common in these headers and was harmless +# until one such comment also contained an indented block. Verbatim runs are NOT # escaped -- Doxygen reproduces them literally, so an escape would show up as a # stray backslash. If you genuinely want a Doxygen command, write a `///` # comment and rule 1 will leave it alone. @@ -89,6 +99,7 @@ function escape(s) { gsub(/&/, "\\&", s) gsub(//, "\\>", s) + gsub(/\|/, "\\|", s) return s } From b3666a7fb4acc260e7f02515078722c18ca8f6f1 Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Wed, 19 Aug 2026 17:00:11 +1000 Subject: [PATCH 7/9] Keep one verbatim display per comment block Doxygen 1.9, which is what CI installs, drops the second `\verbatim` open in a `///` block and reports `unexpected command endverbatim` at a filtered-stream line 120 lines away in unrelated code. 1.17 renders it in silence, so the docs job was the only oracle. The ProfitMax banner had two displays; they are now one. docs.yml asserts the shape and names the block, so the next one fails legibly. The `|`-escaping added earlier is reverted: it was a wrong hypothesis, and master's prose is full of unescaped pipes. --- .claude/CLAUDE.md | 9 +++++ .github/workflows/docs.yml | 47 ++++++++++++++++++++++++++ inst/include/phylloptim/leaf_model.hpp | 12 +++---- tools/doxygen_filter.awk | 36 ++++++++++++-------- 4 files changed, 84 insertions(+), 20 deletions(-) diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index c08d434..f99f769 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -236,6 +236,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. diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 62d8f71..d62a8d9 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -59,6 +59,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/inst/include/phylloptim/leaf_model.hpp b/inst/include/phylloptim/leaf_model.hpp index 852caf3..f31357b 100644 --- a/inst/include/phylloptim/leaf_model.hpp +++ b/inst/include/phylloptim/leaf_model.hpp @@ -3345,14 +3345,14 @@ 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, CG = A(psi)/|A|max, HC = [k(psi_soil)-k(psi)] / [k(psi_soil)-kcrit] -// -// while this package's older entry point maximises `A - lambda*(k(psi_soil)-k(psi))`. +// 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 +// 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 diff --git a/tools/doxygen_filter.awk b/tools/doxygen_filter.awk index 6352901..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 @@ -45,21 +64,11 @@ # preformatted text and stop rendering as lists. # # ESCAPING. Text outside a verbatim run is escaped, because none of it was -# written with Doxygen in mind: `\`, `@`, `#`, `%`, `&`, `<`, `>` and `|` all -# mean something to Doxygen and here they never do. Without it, `#include ` +# written with Doxygen in mind: `\`, `@`, `#`, `%`, `&`, `<` and `>` all mean +# something to Doxygen and here they never do. Without it, `#include ` # becomes a broken link to an entity called "include" followed by a swallowed # HTML tag, and the `\int` in the roots.hpp head-loss note becomes an unknown -# command; both were observed before this was added. -# -# ⚠️ `|` WAS ADDED LAST AND COST A CI-ONLY DEBUGGING ROUND. Doxygen reads a `|` -# in a paragraph as a table delimiter, and `\verbatim` is not allowed inside a -# table cell: it drops the OPEN, and then reports `unexpected command -# endverbatim` at a line that is neither -- in the case that found this, 80 lines -# further on, inside unrelated code. Doxygen 1.17 renders the same input in -# silence, so a local `doxygen` run says nothing; only the 1.9 that CI installs -# objects. The comment that triggered it wrote the maximum of A over the supply -# stream as `|A|max`. `|` in prose is common in these headers and was harmless -# until one such comment also contained an indented block. Verbatim runs are NOT +# command; both were observed before this was added. Verbatim runs are NOT # escaped -- Doxygen reproduces them literally, so an escape would show up as a # stray backslash. If you genuinely want a Doxygen command, write a `///` # comment and rule 1 will leave it alone. @@ -99,7 +108,6 @@ function escape(s) { gsub(/&/, "\\&", s) gsub(//, "\\>", s) - gsub(/\|/, "\\|", s) return s } From a8c65d86a9fa2b57916446c591023abbafe72c60 Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Wed, 19 Aug 2026 20:30:44 +1000 Subject: [PATCH 8/9] Add NEWS for leaf-to-air VPD and ProfitMax The branch changes results on the energy-balance path -- gs and A were overstated by the leaf-to-air deficit ratio, 3.2-4.0x over Tair 25-45 -- and had no NEWS entry, where #89 and #90 both do. Six sections: the VPD fix and its floor, ProfitMax and the unitless members it writes, the four single-layer optimiser fixes, the dgc_dT re-derivation, and the one-verbatim-run-per-block finding. Every number is from the PR comments, checked against the code rather than copied. No behaviour change. Co-Authored-By: Claude Opus 5 (1M context) --- NEWS.md | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/NEWS.md b/NEWS.md index b8b6914..b0d0aaf 100644 --- a/NEWS.md +++ b/NEWS.md @@ -162,6 +162,51 @@ 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](https://github.com/traitecoevo/phylloptim/issues/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 (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 From 517dcb82ca97265487384cf807361226355f9d0a Mon Sep 17 00:00:00 2001 From: Daniel Falster Date: Wed, 19 Aug 2026 20:38:17 +1000 Subject: [PATCH 9/9] Wrap the new NEWS section to match the file The section was written unwrapped, one line per paragraph, beside a file wrapped at ~80: longest prose line 666 characters against 85 in the neighbouring section. Wrapped to 81. Content is unchanged. The heading's issue link becomes a bare `(#7)`, matching the file's other headings, and the full link moves into the prose beside PLAN 13.1 -- a word-level diff of the whole file shows those two moves and nothing else. Co-Authored-By: Claude Opus 5 (1M context) --- NEWS.md | 129 ++++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 102 insertions(+), 27 deletions(-) diff --git a/NEWS.md b/NEWS.md index b0d0aaf..dd2b627 100644 --- a/NEWS.md +++ b/NEWS.md @@ -162,48 +162,123 @@ 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](https://github.com/traitecoevo/phylloptim/issues/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 (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. +## 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**. +`$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. +⚠️ **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. +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. +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".* +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. +`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. +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.