--- title: "Leaf gas exchange from R" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Leaf gas exchange from R} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4, dpi = 120, fig.align = "center") ``` ```{r setup, message = FALSE} library(phylloptim) library(dplyr) library(tidyr) library(ggplot2) theme_set(theme_bw(base_size = 11) + theme(panel.grid.minor = element_blank(), strip.background = element_rect(fill = "grey95", colour = NA))) ``` ## What this model does differently Most leaf gas-exchange models obtain stomatal conductance from an empirical function of assimilation, VPD and CO₂. Ball-Berry, Medlyn and their descendants all work that way. This model has no such function. Instead it builds an explicit hydraulic path from the soil, through the roots and the root collar, up the stem to the leaf. Then it asks which leaf water potential maximises profit, where profit is assimilation minus a hydraulic cost that rises as the xylem loses conductivity. Conductance is whatever falls out of the answer. Stomatal closure under drought is therefore not fitted here. It is derived, from two Weibull vulnerability curves, the soil potential, and the plant's own conductance. ## What is on this page | Section | What it covers | |---|---| | *One solve* to *Light and VPD together* | driving the model and reading the result | | *Traits and numerical settings are separate* | which arguments are physiology and which are numerics | | *Water potentials are positive magnitudes* | the one convention that will break a ported script | | *Two supply paths, two routes* | the root network `plant` uses, and the collar route | | *Trait gradients* | exact derivatives, and when they are worth using | | *A faster method* | the closed form: what it costs in accuracy, and when it refuses | `vignette("the-models")` is the reference for what is being solved. `vignette("fitting")` works through a calibration. ## One solve ```{r} leaf_solve(psi_soil = 2.0, PPFD = 900) |> as_tibble() |> glimpse() ``` Drivers in, operating point out. `psi_soil` is a positive magnitude in MPa, which a later section covers. The defaults are `plant`'s model, `TF24_floor`, which at its default price is TF24, maximised over the leaf potential with the potential upstream of the stem held at `psi_soil` (the *leaf route*). `collar` is `NaN` because that route never solves for a collar potential; the root network `plant` uses, and the *collar route* that solves for one, come later. Most of those columns are the usual gas-exchange quantities. Two are less usual. `lambda` is the marginal cost of water, dA/dE, at the operating point. It is how much carbon the leaf gives up per unit of water it declines to transpire, and it is the quantity the optimisation equalises. `g1_eff` re-expresses the solved conductance as the `g1` of the Medlyn model, which is a convenient common scale if that is the model you are used to. It is computed from the solve rather than supplied, so it varies with conditions. ## A drought response `leaf_solve()` is vectorised over its drivers, so a response curve is one call. ```{r} drought <- leaf_solve(psi_soil = seq(0.2, 5.5, length.out = 40), PPFD = 900) |> as_tibble() drought |> select(psi_soil, A, gc) |> pivot_longer(c(A, gc), names_to = "variable") |> mutate(variable = recode(variable, A = "assimilation~(mu*mol~m^-2~s^-1)", gc = "stomatal~conductance~(mol~m^-2~s^-1)")) |> ggplot(aes(psi_soil, value)) + geom_line(linewidth = 0.8, colour = "steelblue4") + facet_wrap(~variable, scales = "free_y", labeller = label_parsed) + labs(x = "soil water potential (MPa)", y = NULL) ``` Both fall as the soil dries, because each unit of water costs more as the stem nears the end of its vulnerability curve. Neither is held up by a minimum-conductance floor: `g_min` is zero by default. ## The same sweep on a familiar scale `g1_eff` across the same drought, over the range where the stomata are open: ```{r, fig.height = 3.4} drought |> filter(gc > 0, is.finite(g1_eff)) |> ggplot(aes(psi_soil, g1_eff)) + geom_line(linewidth = 0.8, colour = "firebrick") + labs(x = "soil water potential (MPa)", y = expression(effective~g[1]~(kPa^0.5))) ``` It declines as the soil dries, because it is computed from the solve rather than being a parameter of it. ## Light and VPD together Because `leaf_solve()` recycles its drivers elementwise, a factorial design is an `expand_grid()` and a single call. There is no loop, and nothing to join back afterwards, because the drivers come back as columns of the answer. ```{r} expand_grid(PPFD = seq(50, 1800, length.out = 30), atm_vpd = c(0.5, 1, 2, 3)) |> with(leaf_solve(psi_soil = 1.0, PPFD = PPFD, atm_vpd = atm_vpd)) |> ggplot(aes(PPFD, A, colour = factor(atm_vpd))) + geom_line(linewidth = 0.8) + scale_colour_viridis_d(end = 0.9, direction = -1) + labs(x = expression(PPFD~(mu*mol~m^-2~s^-1)), y = expression(A~(mu*mol~m^-2~s^-1)), colour = "VPD (kPa)") ``` The curves saturate at different heights because the hydraulic cost of the transpiration needed to support a given assimilation rises with VPD. ## Traits and numerical settings are separate The C++ constructor takes fifteen positional arguments, four of which are tolerances sitting among the physiology. There are also four traits it does not take at all, which have to be assigned afterwards. `leaf_traits()` and `leaf_control()` split the first fifteen and supply the missing four, so a calibration loop varying traits never has to know either fact. ```{r} vulnerable <- leaf_traits(stem_P50 = 2.2, root_P50 = 2.2) bind_rows( default = leaf_solve(psi_soil = 3.0, PPFD = 900), vulnerable = leaf_solve(psi_soil = 3.0, PPFD = 900, traits = vulnerable), .id = "traits" ) |> select(traits, psi_leaf, A, gc, g1_eff) ``` `stem_P50`/`stem_c` and `root_P50`/`root_c` are two different vulnerability curves. The stem pair drives the hydraulic cost; the root pair drives uptake. They are named explicitly rather than left as an unmarked `b`/`c` pair, because using one where the other belongs is easy to do and produces a plausible wrong answer rather than an error. Tolerances go the other way, through `leaf_control()`: ```{r} tibble( setting = c("default", "tight"), psi_leaf = c( leaf_solve(psi_soil = 2.0, PPFD = 900)$psi_leaf, leaf_solve(psi_soil = 2.0, PPFD = 900, control = leaf_control(GSS_tol_abs = 1e-5))$psi_leaf ) ) ``` Those two are identical, because the solve root-finds its own first-order condition rather than searching profit, to about 1e-12. `GSS_tol_abs` does not reach the leaf route at all. It is read in two places, both on the collar route: the threshold below which the feasible collar interval is treated as a single point, and the tolerance of a fallback search used when neither bracket endpoint has a usable gradient. This is not only a precision story. It is what makes the trait gradients below possible. An argmax resolved only to 1e-3 moves in discontinuous jumps as a trait varies, so its derivative comes back as noise, or worse, as something smooth, plausible and sign-inverted. ## Water potentials are positive magnitudes There is one representation throughout: every ψ is a positive number of MPa, including `psi_soil`, `psi_crit` and the collar potential. Where an equation needs one potential to oppose another, the minus sign is in the equation. This is asserted rather than documented, which matters if you are porting a script that passed signed potentials: ```{r, error = TRUE} leaf_solve(psi_soil = -2.0, PPFD = 900) ``` That is an error rather than a quiet `abs()`. A model accepting both conventions would let such a script return a plausible wrong number instead. ## Two supply paths, two routes Everything above used the single-potential supply, which collapses the whole soil-to-collar journey into one series resistance. The other path is the multi-layer root network that `plant` uses. Both are driven through the same `root_network` argument, because both take their resistances as a per-call driver. On the leaf route that resistance never reaches the solve: the potential upstream of the stem is pinned at `psi_soil`. The **collar route** maximises over the root-collar potential instead, so water has to be pulled through the roots, and the resistance costs the leaf. Every cost curve takes either route; the collar route is the only one on the multi-layer supply. ```{r} leaf_solve(psi_soil = 1.5, PPFD = 900, route = "collar", root_network = series_resistance(1e3)) |> as_tibble() |> select(psi_soil, psi_leaf, collar, A, gc, g1_eff) ``` The resistance is per unit leaf area, because nothing in this model is allowed to scale with plant size. On the collar route, raising it makes water harder to get, and the leaf closes down: ```{r, fig.height = 3.4} tibble(resistance = 10^seq(2.5, 4.5, length.out = 25)) |> rowwise() |> mutate(A = leaf_solve(psi_soil = 1.5, PPFD = 900, route = "collar", root_network = series_resistance(resistance))$A) |> ungroup() |> ggplot(aes(resistance, A)) + geom_line(linewidth = 0.8, colour = "steelblue4") + scale_x_log10() + labs(x = expression(soil-to-collar~resistance~(MPa~s~mol^-1~m^2)), y = expression(A~(mu*mol~m^-2~s^-1))) ``` The single potential is what makes comparison against other optimality models meaningful. Medlyn, Prentice least-cost and Cowan-Farquhar are all formulated against a single soil potential, so comparing them against a multi-layer root network would compare two things at once. The supply is chosen when the leaf is built rather than afterwards. There is no settable `leaf$supply_kind`: ```{r, error = TRUE} l <- leaf_model() l$supply_kind <- "multilayer" ``` Flipping the tag alone would leave the other path's state configured and silently ignored. Both `leaf_supply_*()` functions therefore reconfigure the object completely, so it can never be in a state where the tag and the supply disagree. ## Multiple soil layers On the multi-layer supply, pass a vector for a layered profile, or a list of vectors for several profiles. Layers default to 1 m thick with the root carbon split evenly. ```{r} leaf_solve(psi_soil = list(c(1.0, 1.5, 2.0), c(3.0, 3.5, 4.0)), PPFD = 900, supply = leaf_supply_multilayer(), route = "collar") |> as_tibble() |> select(psi_soil, layers, A, gc, uptake) ``` `psi_soil` in the result is the wettest layer, with `layers` alongside, so the column stays numeric and the frame stays rectangular. ## The stateful interface `leaf_solve()` is a wrapper over an R6 object you can drive yourself. This is what `plant` does, and it is what you want when you care about the intermediate state rather than only the operating point. It is also the cheaper route per row: driving one object by hand, as below, measures at about half of `leaf_solve()`'s cost per row. Both are dominated by the same thing, which is that a call across the R boundary costs about half as much as the whole C++ solve. If you need an inner loop faster than that, what has to go is the number of R calls per row rather than which R function makes them. The measured figures are in the README's [Performance from R](https://github.com/traitecoevo/phylloptim#performance-from-r). ```{r} l <- leaf_model() set_drivers(l, psi_soil = 2.0, PPFD = 900) l$optimise() tibble(profit = l$profit_, cost = l$hydraulic_cost_, A = l$assim_colimited_) ``` The object exposes the whole model, not just the operating point. The stem vulnerability curve is defined for any potential: ```{r, fig.height = 3.4} tibble(psi = seq(0, 6, length.out = 200)) |> mutate(k = vapply(psi, l$proportion_of_conductivity, numeric(1))) |> ggplot(aes(psi, k)) + geom_line(linewidth = 0.8, colour = "darkgreen") + labs(x = expression(psi~(MPa)), y = "proportion of conductivity", title = "stem vulnerability") ``` The profit function is not defined everywhere. Water flows from the soil to the leaf, so `psi_leaf` must be at least the upstream potential. Below it the leaf would be wetter than the soil, so `transpiration()`, `stom_cond_CO2()` and `psi_leaf_to_ci()` stop with an error rather than return a negative flux. `profit_psi_leaf_TF24()` reads that pair as zero flux, as the solve's evaluator does, so it returns a number, but not one on the curve the solve maximises. `psi_leaf` equal to the upstream potential is the shut stomata's state, with no stomatal flux. Plot the domain the solve actually searches: ```{r, fig.height = 3.4} psi_upstream <- 2.0 tibble(psi = seq(psi_upstream, 6, length.out = 200)) |> mutate(profit = vapply(psi, l$profit_psi_leaf_TF24, numeric(1), psi_upstream = psi_upstream)) |> ggplot(aes(psi, profit)) + geom_line(linewidth = 0.8, colour = "darkgreen") + geom_vline(xintercept = l$opt_psi_leaf_, linetype = 2, colour = "grey40") + labs(x = expression(psi[leaf]~(MPa)), y = "profit", title = "the function being maximised") ``` The dashed line is where the solve landed, and the maximum it sits on is flat. That flatness has a practical consequence: the profit *value* is much better determined than its *location*. An error `dp` in the value displaces the location by `sqrt(dp/k)` for curvature `k`, so locating the peak by comparing heights costs precision that solving `dprofit = 0` keeps. It is why the solve does not search this curve. `vignette("the-models")` Appendix C has the measurements. ## Trait gradients, and where they stop working Because the operating point is now located precisely, it can be *differentiated*. `leaf_gradient()` returns the derivative of each solved output with respect to each trait: ```{r} g <- leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = c("vcmax_25", "jmax_25", "stem_P50", "TF24_cost_scale")) round(g$gradient, 5) ``` ### Choosing which gradients to compute `pars` selects them, and by default you get all of them. Omit it and `leaf_gradient()` differentiates every parameter it can: ```{r} dim(leaf_gradient(psi_soil = 2.0, PPFD = 900)$gradient) rownames(leaf_gradient(psi_soil = 2.0, PPFD = 900)$gradient) ``` That is `r length(leaf_traits())` traits from `leaf_traits()`, plus `leaf_specific_conductance_max`, which is not a trait but is fitted like one, and `TF24_floor_lambda_o`, the default model's price of water: `r length(leaf_traits()) + 2L` in all. The collar route over the single-potential path adds `resistance`, its series resistance, which the leaf route never reads; `model = "TF24"` drops the price, and `model = "CF77"` swaps in `CF77_lambda_`, that curve's prescribed marginal value of water. These counts are computed rather than written out, for the reason the hidden chunk below records. `R_d_25`, dark respiration at 25 °C, is one of the traits and behaves like any other. So `dY/dvcmax_25` is a partial derivative at fixed respiration. A fit that wants respiration to follow Vcmax moves both and adds the two columns. Selecting fewer parameters costs no less. The work is one solve plus one pass that carries the derivative in every parameter at once, so asking for one parameter costs what asking for all of them does. ```{r} #| echo: false bench <- function(pars) { t0 <- Sys.time() for (i in 1:20) leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = pars) as.numeric(Sys.time() - t0, units = "secs") / 20 * 1000 } tibble::tibble( # Counted, not written down: this label said 16 for a while when the answer was # 14, and the prose above it went on saying 14 after the answer became 16. That # is exactly the way a number beside a chunk rots -- the chunk runs, the prose # next to it is never asserted -- so the prose uses inline R now too. request = c(sprintf("`pars = NULL` (all %d)", nrow(leaf_gradient(psi_soil = 2.0, PPFD = 900)$gradient)), "4 parameters", "1 parameter"), `ms per call` = round(c(bench(NULL), bench(c("vcmax_25", "jmax_25", "stem_P50", "TF24_cost_scale")), bench("vcmax_25")), 2) ) |> knitr::kable() ``` So name the parameters you are fitting for the shape of the result, not for speed. The five outputs are not selectable. Every call returns the derivative of `A`, `gc`, `psi_leaf`, `collar` and `profit` with respect to each requested parameter, because they are read off the same solve and cost nothing extra. `pars` chooses the rows of that matrix, never the columns. The first four are what a gas-exchange calibration observes. `profit` is the objective, and it is there for a demographic consumer: `plant` bills carbon from the leaf's `profit` rather than its assimilation. It is also the one output the envelope theorem reaches, so at an interior optimum its gradient is the direct partial alone. `?leaf_gradient` explains that, along with the measurement that makes it worth doing rather than merely true. ### A collar potential you supply If your model carries the collar potential rather than solving for it each step, as `plant` does, pass it and the gradient is taken there. This and the rest of the section are on the collar route over the multi-layer supply: ```{r} collar <- list(supply = leaf_supply_multilayer(), route = "collar") g <- leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "vcmax_25", psi = 2.7, supply = collar$supply, route = collar$route) g$status ``` `psi` is exogenous, so the answer is the partial derivative at fixed collar, plus `dpsi_dtheta` if the collar you imposed itself moves with the traits. `M` and `dY_dpsi` come back too, so a caller integrating its own sensitivity of `psi` has the coefficients. Hand back the collar the solver found, with the response it derived, and the two paths agree to rounding: ```{r} a <- leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "vcmax_25", supply = collar$supply, route = collar$route) b <- leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "vcmax_25", psi = a$psi, dpsi_dtheta = -a$M / a$H, supply = collar$supply, route = collar$route) all.equal(a$gradient, b$gradient, tolerance = 1e-12) ``` A `psi` outside the feasible collar interval is refused: the leaf does not move it onto the interval's end, where the outputs would follow the *bound* rather than your `dpsi_dtheta`. `status` reports `"refused"`, `g$psi` is the one you gave, and everything else is `NA`. A feasible `psi` whose leaf potential lies in the stem curve's tail, just inside the dry end, is `"no-gradient"`: it has a value, but every row would be rounding there. Naming something that cannot be differentiated is an error rather than a silent omission, and the message says what is available instead: ```{r, error = TRUE} leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "resistance", supply = collar$supply, route = collar$route) ``` These are not finite differences of the solve. The outputs are read at the profit-maximising collar potential, so a trait moves them twice over: directly, and by moving the optimum itself. Differentiating the optimality condition `dprofit/dψ = 0`, rather than the solved output, captures both terms, and one pass of automatic differentiation carries every parameter at once. For `TF24_cost_scale`, `TF24_beta2`, `stem_P50` and `stem_c` the second term is the entire answer. Change any of them and assimilation moves only because the leaf chooses a different operating point. The interesting part is where that derivation stops being valid. It assumes the optimum is interior, and in dry soil it often is not: the leaf ends up at the very edge of the range of collar potentials it can physically occupy. There `dprofit/dψ` is not zero at the answer, and the formula is not merely imprecise but wrong by orders of magnitude, without failing loudly. The solve records which kind of point it found, and at a pinned one `leaf_gradient()` differentiates the bound that holds the leaf instead: ```{r, fig.height = 3.4} sweep <- lapply(seq(1.0, 5.0, by = 0.125), function(p) { g <- leaf_gradient(psi_soil = p, PPFD = 900, atm_vpd = 3.0, pars = "stem_P50", supply = collar$supply, route = collar$route) tibble(psi_soil = p, dA_dstem_P50 = g$gradient["stem_P50", "A"], status = g$status) }) |> bind_rows() ggplot(sweep, aes(psi_soil, dA_dstem_P50)) + geom_line(colour = "grey60") + geom_point(aes(colour = status), size = 2) + labs(x = expression(psi[soil]~(MPa)), y = expression(dA/d*italic(stem_P50)), colour = NULL, title = "sensitivity to stem vulnerability, and where the optimum pins") ``` On the leaf route this sweep stays interior throughout: the potential upstream of the stem is pinned at `psi_soil`, so there is no collar interval to pin against. Read that left to right. In wet soil, a less vulnerable stem is worth about two units of assimilation per MPa of `stem_P50`. The sensitivity rises as the soil dries, because vulnerability matters more when water is scarce, peaks near 2.9 MPa, eases off, and then falls away sharply. Past about 3.9 MPa the optimum is a bound: the leaf is stuck at the wettest collar potential it can physically reach, and a stem trait it has no freedom to exploit stops mattering at all. The tail is exactly zero. The transition itself is a genuine kink in trait space, and a real non-differentiable set. It belongs to the constrained problem rather than to this implementation. A gradient-based fit that wanders across it is optimising through a corner, and no amount of care in the derivative changes that. What `leaf_gradient()` can do is tell you when you are there, which a naive finite difference cannot: `status` marks these points, from the kind of point the solve recorded. ## A faster method: the closed form Everything above reaches the operating point the same way: a root-find of the first-order condition over the feasible interval. That is the `method = "exact"` path, and it is the default. `set_model()` takes a third argument, and `"closed"` selects a different route to the same optimum. The idea is the one `vignette("the-models")` builds the whole comparison on. Every model here satisfies `dA/dE = lambda` at its optimum, and *given* `lambda` the answer collapses to the Medlyn relation `ci/ca = xi/(xi + sqrt(D))` with `xi = sqrt(Q/lambda)`. The difficulty is that `lambda` depends on the `psi` you are solving for — but where it is a power law in `psi` at the wet end, that circularity inverts explicitly. One Newton step on the supply-minus-demand residual finishes the job, and a search over the bracket is not needed at all. ```{r} cf <- leaf_model(supply = leaf_supply_singlelayer()) set_drivers(cf, psi_soil = 0.5, PPFD = 1500, atm_vpd = 1, root_network = series_resistance(1500)) cf$set_model("TF24", "leaf", "closed") cf$optimise() tibble(method = cf$model_method(), psi_leaf = cf$opt_psi_leaf_, A = cf$assim_colimited_, fell_back = cf$last_solve_fell_back_) ``` ### When it applies, and when it refuses It is refused rather than approximated wherever it does not hold, and the refusals carry their reasons: - **`TF24`, `CF77` and `TF24_floor` only, of the nine cost curves.** The three-way split is worth knowing, because only one part of it is a gap in the code. A closed form exists exactly when `h'(A)`, the derivative of the benefit link, does not depend on the solution. For `JS22`, `CMax` and `ProfitMax` it does not — so the inversion *exists* for them and is simply not written yet. For `SOX` and `JW26` it does: their link is `log`, so `h'(A) = 1/A`, and `lambda` therefore carries the assimilation the solve is for. Given `lambda -> xi -> ci -> A` the dependency closes into a fixed point, which can be iterated but not inverted. Those two **cannot** have a closed form. ⚠️ **`TF24_floor` needs `$TF24_floor_lambda_o > 0`, and its tail is worse than the other two's.** Its wet-end `lambda` is `lambda_o` — a constant, so the same `n = 0` case as `CF77` — and at `lambda_o = 0` that limit is zero, `xi = sqrt(Q/lambda)` diverges, and there is nothing to expand about; that case is `TF24`, which has its own start, and the refusal points there. Where it does apply, the leading order discards the whole `Θ̃'(ψ)/K(ψ)` term, which at the DEFAULT `TF24_cost_scale` dominates `lambda_o` at the operating point. Measured over a 150-row driver grid, the worst served error is 2.9e-01 at `TF24_cost_scale = 0.5` (against `CF77`'s own 2.4e-01) and **5.2e-01** at the default 7.5. The `ci/ca` guard does not catch that, because it reads the closed form's own `lambda` and stays plausible while `psi*` is wrong — so on this curve the fast arm is a measurement tool rather than a production path. - **The leaf route only, and one soil layer only.** The inversion is against a supply measured from `psi_soil`; it says nothing about the soil-to-collar path, and the multi-layer closed form does not exist. - **Not with the energy balance on.** Fick's law divides by the leaf-to-air deficit, which depends on leaf temperature, which depends on transpiration — which is what the inversion solves for. So `D` is not known at the moment of inversion and the form is not closed. This is a statement about the form rather than about its accuracy: measured with the refusal bypassed, the error reaches 2.3 MPa in `psi_leaf` even in the wettest bin, against 0.04 MPa on the prescribed-temperature path. ```{r, error = TRUE} cf$set_model("SOX", "leaf", "closed") cf$optimise() ``` ### What it costs, and what it buys The guard on validity tests an **output** — `ci/ca` at the answer — so it can only be applied after solving, and a rejected row falls back to the exact solve. That makes the fallback fraction `phi` part of the price: the realised speedup is `1/[phi + (1 - phi)/s]` rather than the ceiling `s`. Measured in `tests/cpp/bench_solve` at 0.16.2, with both arms in one process on a 96-point single-layer driver grid: | curve | realised speedup | `phi` | ceiling at `phi = 0` | |---|---|---|---| | `TF24` | **2.2×** | 0.27 | 4.0× | | `CF77` | **2.0×** | 0.39 | 5.3× | Two-fold, in other words — not the order of magnitude an earlier version of this code advertised. Those figures predated the current solver, which root-finds the first-order condition where it used to search the objective, and they also priced the inversion alone rather than a solve a caller can read outputs from. `phi` is reported per leaf, so you can measure it on your own drivers rather than taking the grid's: ```{r} sweep_leaf <- leaf_model(supply = leaf_supply_singlelayer()) sweep_leaf$set_model("TF24", "leaf", "closed") for (p in seq(0.25, 5, length.out = 20)) { set_drivers(sweep_leaf, psi_soil = p, PPFD = 1500, atm_vpd = 1, root_network = series_resistance(1500)) sweep_leaf$optimise() } tibble(calls = sweep_leaf$closed_form_calls_, fell_back = sweep_leaf$closed_form_fallbacks_, phi = sweep_leaf$closed_form_fallback_fraction()) ``` ### Accuracy, and one failure that is not a tolerance Three approximations sit inside the closed form: the power law for `lambda`, the single Newton step, and the Medlyn collapse itself — whose `3*Gstar` group is the electron-transport-limited `dA/dci`, so it is not exact under colimitation. The last of those is the binding one, and the error is organised by `ci/ca`: ```{r} grid <- expand.grid(psi_soil = c(0.25, 0.5, 1, 2, 3, 4, 5), PPFD = c(100, 500, 1500), atm_vpd = c(0.5, 1, 2, 3)) solve_both <- function(psi_soil, PPFD, atm_vpd) { one <- function(method) { l <- leaf_model(supply = leaf_supply_singlelayer()) set_drivers(l, psi_soil = psi_soil, PPFD = PPFD, atm_vpd = atm_vpd, root_network = series_resistance(1500)) l$set_model("TF24", "leaf", method) l$optimise() l } ex <- one("exact") cl <- one("closed") tibble(ci_ca = ex$ci_ / ex$ca_, fell_back = cl$last_solve_fell_back_, err_A = abs(cl$assim_colimited_ - ex$assim_colimited_) / max(1, abs(ex$assim_colimited_)), err_psi = abs(cl$opt_psi_leaf_ - ex$opt_psi_leaf_)) } res <- Map(solve_both, grid$psi_soil, grid$PPFD, grid$atm_vpd) |> bind_rows() |> filter(is.finite(ci_ca)) res |> mutate(bin = cut(ci_ca, c(0, 0.25, 0.4, 0.5, 0.6, 0.8, 1))) |> group_by(bin) |> summarise(n = n(), fell_back = sum(fell_back), mean_err_A = mean(err_A), worst_err_A = max(err_A), worst_err_psi = max(err_psi), .groups = "drop") ``` Read it as three regimes. Above `ci/ca` about 0.6 the closed form agrees with the exact solve to well under a percent in assimilation on average, and a few percent at worst. Between 0.5 and 0.6 the mean is a couple of percent. And below about 0.4 the guard has fired on nearly every row, so the tabulated error collapses to the solver's own floor — those rows *are* the exact solve. ⚠️ **The failure below `ci/ca` about 0.3 is topological, not a matter of tolerance.** There the exact optimum turns back toward shutdown while the inversion runs on to `psi_crit`: the fixed point the closed form is looking for has ceased to exist, so the answer is wrong rather than imprecise, and no amount of refinement recovers it. That is why there is a guard and a fallback at all rather than a tighter convergence criterion — and it is why `$last_solve_fell_back_` is worth reading if you are going to trust an individual row. ### The caveat that would stop you using this in a demographic model `plant` differentiates the operating point with respect to plant state, and that needs the argmax to vary *smoothly* with its inputs. Sweeping `kmax` — which is how height reaches the leaf — over ±30% and taking the mean absolute second difference of `opt_psi_leaf_`: | the sweep | exact | closed | |---|---|---| | never crosses the guard | 6.1e-06 | 7.2e-06 | | crosses the guard | 1.0e-05 | 1.0e-03 | The closed form's own arithmetic is barely rougher than the solve, at about 1.2×. What costs two orders of magnitude is the **fallback boundary**: the guard tests an output, so the set of drivers on which the method switches is a surface in driver space, and the argmax jumps as a smooth input carries the leaf across it. A mixture of two methods is discontinuous even where both are individually smooth. So this is usable for a large batch of independent operating points, and it is **not** usable underneath a derivative. `leaf_gradient()` reflects that: seating a model resets the method to `"exact"`, so a gradient always differentiates the exact solve whatever was seated before. ## What the model assumes Every model simplifies. The assumptions most worth keeping in mind when reading output from this one are: - **Stomatal behaviour is optimal and instantaneous.** The leaf sits at the maximum of its objective at every point, with no cost of adjustment and no lag. - **The leaf is intensive.** Every input is per unit leaf area, or is a dimensionless driver. Nothing scales with plant size, so whole-plant allometry belongs to the caller. - **The hydraulic path is two Weibull curves in steady state.** There is no capacitance and no refilling. - **The default cost curve is one hypothesis among several.** `TF24` is the default because it is what `plant` runs, not because it has been shown to be the right account of what a leaf is spending. `vignette("the-models")` has the other six and how they differ. - **The default traits are *Eucalyptus saligna*.** They are a starting point rather than a recommendation, and the same is true of both default supply networks. ## Where to go next - **I want to know what is being solved.** `vignette("the-models")` is the reference: the objective, the nine cost curves, the first-order conditions, and the two optional gates. - **I want to fit this to data.** `vignette("fitting")` works through a calibration, and covers when the exact gradient earns its complexity. - **I want to call it without R.** The model is a set of self-contained C++ headers, so the same solve is available from a C++ program or a Python extension. `vignette("cpp-interface")` covers it. - **I want to know how this differs from what I already use.** `COMPARISON.md` sets it against `plantecophys`, `tealeaves` and `bigleaf`. One loose end worth flagging if you are supplying your own hydraulics. Both supply paths take resistances through the same `root_network` argument, so a caller with measured or fitted resistances never has to go through a root carbon model. The single exception is `gravity_head`, which `leaf_supply_singlelayer()` still takes as configuration, because that path has no depth profile to derive it from.