diff --git a/NEWS.md b/NEWS.md index 0ed7814..b8b6914 100644 --- a/NEWS.md +++ b/NEWS.md @@ -162,6 +162,72 @@ 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. +## A trait gradient at a collar potential the caller supplies + +`leaf_gradient()` and `leaf_gradient_batch()` take `psi`, and evaluate there +instead of solving for the profit-maximising collar +([#88](https://github.com/traitecoevo/phylloptim/issues/88)). Until now both +called `find_root_collar_psi()` unconditionally and read `opt_root_psi_` back, +so a model that **tracks** the optimum rather than finding it — plant's TF24f +carries the collar as an ODE state, `dpsi/dt = k * dprofit/dpsi` — could not ask +this package for a trait gradient at the point it was actually operating at. It +got a confident answer about the re-solved optimum instead, with nothing saying +so. + +The maths simplifies rather than breaks: `psi` is exogenous, so the indirect term +is whatever the caller says it is, via `dpsi_dtheta` (default zero, the partial +at fixed collar). Nothing is derived from `-M/H`, so nothing needs stationarity, +and `method` is refused — the two routes it chooses between are both about a +solved optimum. + +`M`, `H`, `dY_dpsi` and `psi` now come back in the result on both paths. A caller +whose `psi` is *dynamic* cannot supply `dpsi_dtheta` as a constant: for the +gradient-ascent law above it obeys `ds/dt = k(M + H s)`, and those are its +coefficients (traitecoevo/plant#614). + +⚠️ **`stationarity` is still computed on the prescribed path, and now means +something better.** It no longer routes anything — it measures how far the collar +you supplied sits from the optimum. It makes exactly one decision, `profit`'s: +at a stationary point the envelope theorem applies and the analytic zero is used; +away from one the *exact* `dprofit/dpsi` is used rather than a difference of it. +One rule, both paths — which is why `psi = ` with +`dpsi_dtheta = -M/H` reproduces the solving path **bit-for-bit**, asserted with +`identical()` rather than a tolerance. + +⚠️ **A clamped `psi` returns no gradient, and this is the case to understand.** +The collar actually used is `psi` clamped into the feasible interval, so it moves +with the *bound* rather than with `dpsi_dtheta` — the active-set problem arriving +through the clamp instead of through the optimiser, where the direct term alone +is plausible and wrong. `status` reports `"clamped"`, the gradient is `NA`, and +`psi` in the result is the collar that was used. Reported rather than thrown +because a tracking model reaches these points routinely: the clamp is how TF24f +pulls an out-of-range state back inside. It also fires for a `psi` within one +step of an end, where `dY/dpsi` cannot be centred. + +The solving path is unchanged and bit-identical, including `gradient_golden.tsv`. + +⚠️ **An INFEASIBLE prescribed `psi` is `"no-gradient"`, not a sentinel zero.** +`dprofit_droot_collar_psi` returns a hard `0.0` on its shut-down and +reversed-gradient exits, and a bare zero is indistinguishable from a stationary +point. The solving path got away with reading the value alone because `H` +collapses to zero with it and `usable` catches the pair; the prescribed path +never divides by `H`, so it would have adopted the sentinel *as* `dprofit/dpsi` +— silently losing profit's indirect term at exactly the dry points a tracking +model lives in. Most such points are caught as `"clamped"` first, but not the one +where the caller hands back the collar the shut-down state itself seated. + +**`Leaf$dprofit_droot_collar_psi_checked()` is new and is what makes that +possible.** The `bool* feasible` out-parameter has been there since #79 and the +C++ vignette has always said a composite ignoring it inherits the bug — but +RcppR6 has no form for a `bool*`, so the generated binding dropped it and every +R-side composite *was* that composite. It returns `{dprofit, feasible}`. + +C++ consumers get the same through `gradient::Prescribed` and the new +`psi`/`dpsi_dtheta` arguments to `gradient::batch`. `Status` gains `Prescribed` +and `Clamped` — ⚠️ **appended after `Error`, so no existing integer value +moves**, and `status_name`'s switch is exhaustive with no `default:` so the next +member added is a compiler diagnostic rather than a silent `"error"` label. + ## The gradient differentiates `profit`, which is what a demographic caller bills `leaf_gradient()` and `leaf_gradient_batch()` return a fifth column. The four diff --git a/R/RcppExports.R b/R/RcppExports.R index db96ca4..50d682b 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -105,6 +105,10 @@ Leaf__dprofit_droot_collar_psi <- function(obj_, opt_root_psi) { .Call('_phylloptim_Leaf__dprofit_droot_collar_psi', PACKAGE = 'phylloptim', obj_, opt_root_psi) } +Leaf__dprofit_droot_collar_psi_checked <- function(obj_, opt_root_psi) { + .Call('_phylloptim_Leaf__dprofit_droot_collar_psi_checked', PACKAGE = 'phylloptim', obj_, opt_root_psi) +} + Leaf__psi_stem_to_ci <- function(obj_, psi_stem, psi_upstream) { .Call('_phylloptim_Leaf__psi_stem_to_ci', PACKAGE = 'phylloptim', obj_, psi_stem, psi_upstream) } @@ -802,8 +806,8 @@ gradient_batch_check <- function(drivers) { invisible(.Call('_phylloptim_gradient_batch_check', PACKAGE = 'phylloptim', drivers)) } -gradient_batch_run <- function(obj_, drivers, theta, pars, step, stationarity_tol, method, fast_stem_curve) { - .Call('_phylloptim_gradient_batch_run', PACKAGE = 'phylloptim', obj_, drivers, theta, pars, step, stationarity_tol, method, fast_stem_curve) +gradient_batch_run <- function(obj_, drivers, theta, pars, step, stationarity_tol, method, fast_stem_curve, psi, dpsi_dtheta) { + .Call('_phylloptim_gradient_batch_run', PACKAGE = 'phylloptim', obj_, drivers, theta, pars, step, stationarity_tol, method, fast_stem_curve, psi, dpsi_dtheta) } #' Root resistances from a root carbon profile diff --git a/R/RcppR6.R b/R/RcppR6.R index f798ac6..30f4472 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: b33a1650c9cebdd0b589a1e1622ce834 ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -162,6 +162,9 @@ NULL dprofit_droot_collar_psi = function(opt_root_psi) { Leaf__dprofit_droot_collar_psi(self, opt_root_psi) }, + dprofit_droot_collar_psi_checked = function(opt_root_psi) { + Leaf__dprofit_droot_collar_psi_checked(self, opt_root_psi) + }, psi_stem_to_ci = function(psi_stem, psi_upstream) { Leaf__psi_stem_to_ci(self, psi_stem, psi_upstream) }, diff --git a/R/gradient-batch.R b/R/gradient-batch.R index 5f641cc..fe2e0ca 100644 --- a/R/gradient-batch.R +++ b/R/gradient-batch.R @@ -187,6 +187,17 @@ print.leaf_batch <- function(x, ...) { ##' [leaf_batch()]. One crossing of the R boundary per call, in place of 112 per ##' observation. ##' +##' @param psi one collar water potential per observation to evaluate at +##' (positive magnitudes, MPa) instead of solving for the profit-maximising +##' one. See [leaf_gradient()]'s section on it. Recycling a single value across +##' the batch is deliberately not offered — rows sharing an imposed collar is a +##' coincidence rather than a design, and broadcasting one would make a length +##' mismatch look intentional. +##' @param dpsi_dtheta how each prescribed `psi` responds to each parameter: +##' an `n` × `length(pars)` matrix, or one vector of `length(pars)` shared by +##' every observation. Defaults to zero. Must be finite, for the reason +##' [leaf_gradient()] records. Only meaningful with `psi`. +##' ##' @section What it computes, and where the maths is written down: ##' The same five derivatives at the same solved operating point, by the same two ##' routes, with the same active-set test deciding between them. **Read @@ -242,10 +253,15 @@ print.leaf_batch <- function(x, ...) { ##' \item{`gradient`}{an `n` × `length(pars)` × 5 array, the last dimension ##' being `A`, `gc`, `psi_stem`, `collar` and `profit`} ##' \item{`value`}{an `n` × 5 matrix of the solved outputs} -##' \item{`method`}{`"ift"` or `"fd"` per observation} -##' \item{`status`}{`"interior"`, `"pinned"`, `"no-gradient"` or `"error"`} +##' \item{`method`}{`"ift"`, `"fd"` or `"prescribed"` per observation} +##' \item{`status`}{`"interior"`, `"pinned"`, `"no-gradient"`, `"error"`, +##' or with `psi` one of `"prescribed"`, `"clamped"`, `"no-gradient"`} ##' \item{`H`, `stationarity`}{the curvature and the implied Newton step ##' `method` was decided on, per observation} +##' \item{`M`}{an `n` × `length(pars)` matrix of mixed partials +##' `d2profit/dpsi dtheta`} +##' \item{`dY_dpsi`}{an `n` × 5 matrix of `dY/dpsi` at fixed traits} +##' \item{`psi`}{the collar each row was evaluated at} ##' \item{`message`}{`""`, or why that row failed} ##' } ##' @@ -261,6 +277,8 @@ leaf_gradient_batch <- function(batch, traits = NULL, pars = NULL, theta = NULL, + psi = NULL, + dpsi_dtheta = NULL, step = 1e-6, stationarity_tol = 1e-8, method = c("auto", "ift", "fd"), @@ -310,17 +328,96 @@ leaf_gradient_batch <- function(batch, theta <- .gradient_check_theta(theta, batch, par_names) } + psi <- .gradient_check_psi_batch(psi, batch$n, method) + dpsi_dtheta <- .gradient_dpsi_dtheta_batch(dpsi_dtheta, pars, batch$n, + !is.null(psi)) + res <- gradient_batch_run(batch$leaf, batch$drivers, theta, match(pars, par_names) - 1L, step, - stationarity_tol, method, fast_stem_curve) + stationarity_tol, method, fast_stem_curve, + psi, dpsi_dtheta) dimnames(res$gradient) <- list(NULL, pars, .gradient_output_names()) dimnames(res$value) <- list(NULL, .gradient_output_names()) + dimnames(res$M) <- list(NULL, pars) + dimnames(res$dY_dpsi) <- list(NULL, .gradient_output_names()) res } # --- internals --------------------------------------------------------------- +# The batch counterparts of `.gradient_check_psi()` and +# `.gradient_dpsi_dtheta()`: one collar potential per observation, and one +# dpsi/dtheta per observation per parameter. Recycling a length-1 `psi` across a +# batch is deliberately NOT offered -- a batch whose rows share an imposed collar +# is a coincidence rather than a design, and silently broadcasting one would make +# a length mismatch look intentional. +.gradient_check_psi_batch <- function(psi, n, method) { + if (is.null(psi)) { + return(NULL) + } + if (!identical(method, "auto")) { + stop("`method` cannot be given with `psi`; see ?leaf_gradient.", + call. = FALSE) + } + if (!(is.numeric(psi) && length(psi) == n)) { + stop("`psi` must be one collar potential per observation (", n, + " for this batch); got ", length(psi), call. = FALSE) + } + if (!all(is.finite(psi) & psi > 0)) { + stop("`psi` must be finite and positive: it is a collar water potential ", + "in MPa, as a positive magnitude.", call. = FALSE) + } + as.numeric(psi) +} + +.gradient_dpsi_dtheta_batch <- function(dpsi_dtheta, pars, n, prescribed) { + if (is.null(dpsi_dtheta)) { + # NULL is zero on the C++ side, so nothing is allocated for the common case + # of a partial derivative at fixed collar. + return(NULL) + } + if (!prescribed) { + stop("`dpsi_dtheta` is the trait response of a collar potential YOU ", + "imposed, so it needs `psi`.", call. = FALSE) + } + # Finite, for the reason `.gradient_dpsi_dtheta()` records: an infinite value + # returns an all-Inf row with a finite, plausible `profit` in it. Checked here + # as well as there because the MATRIX form never reaches that function -- only + # the named-vector form does, further down. + if (!is.numeric(dpsi_dtheta) || !all(is.finite(dpsi_dtheta))) { + stop("`dpsi_dtheta` must be numeric and finite -- no NA, NaN or Inf. It ", + "is dpsi/dtheta for a collar YOU imposed; a non-finite value returns ", + "an all-infinite gradient row with a plausible `profit` in it.", + call. = FALSE) + } + if (is.null(dim(dpsi_dtheta))) { + if (length(dpsi_dtheta) != length(pars)) { + stop("`dpsi_dtheta` as a vector is one value per parameter (", + length(pars), "), shared by every observation; got ", + length(dpsi_dtheta), ". Pass an ", n, " x ", length(pars), + " matrix to vary it by observation.", call. = FALSE) + } + if (!is.null(names(dpsi_dtheta))) { + dpsi_dtheta <- .gradient_dpsi_dtheta(dpsi_dtheta, pars, TRUE) + } + return(matrix(dpsi_dtheta, nrow = n, ncol = length(pars), byrow = TRUE)) + } + if (!identical(dim(dpsi_dtheta), c(n, length(pars)))) { + stop("`dpsi_dtheta` must be ", n, " x ", length(pars), + " (observations x `pars`); got ", + paste(dim(dpsi_dtheta), collapse = " x "), call. = FALSE) + } + if (!is.null(colnames(dpsi_dtheta))) { + if (!setequal(colnames(dpsi_dtheta), pars)) { + stop("`dpsi_dtheta` has column names, so they must be exactly `pars`.", + call. = FALSE) + } + dpsi_dtheta <- dpsi_dtheta[, pars, drop = FALSE] + } + matrix(as.numeric(dpsi_dtheta), nrow = n, ncol = length(pars)) +} + # `theta` from a leaf_traits() plus the batch's own conductance and resistance. # The counterpart of `.gradient_theta()`, which builds the same vector for one # observation -- and it reads the two non-traits back out of the resolved drivers diff --git a/R/gradient.R b/R/gradient.R index 6f2a5c8..b9f5a6e 100644 --- a/R/gradient.R +++ b/R/gradient.R @@ -159,6 +159,54 @@ set_traits <- function(x, traits) { ##' `"interior"`; everywhere else `profit` goes through the same finite-difference ##' fallback as the other four. ##' +##' @section A collar potential you supply, instead of the one this solves for: +##' Pass `psi` and the outputs are evaluated at that collar rather than at the +##' argmax. This is for a caller whose model **tracks** the optimum instead of +##' finding it — plant's TF24f carries the collar potential as an ODE state, +##' `dpsi/dt = k * dprofit/dpsi`, so at finite gain its `dprofit/dpsi` is +##' deliberately non-zero and *is* the acclimation rate. +##' +##' There the derivation above simplifies rather than breaks. `psi` is exogenous, +##' so the indirect term is whatever the caller says it is: +##' +##' \deqn{dY/d\theta = \partial Y/\partial\theta|_\psi + +##' (\partial Y/\partial\psi)(d\psi/d\theta)} +##' +##' with `dpsi_dtheta` supplied — defaulting to zero, the partial at fixed +##' collar. Nothing is derived from `-M/H`, so nothing needs stationarity, and +##' `method` is refused because the two routes it chooses between are both about +##' a solved optimum. +##' +##' `M`, `H` and `dY_dpsi` come back in the result, because a caller with a +##' *dynamic* `psi` cannot supply `dpsi_dtheta` as a constant: for the +##' gradient-ascent law above it obeys `ds/dt = k(M + H s)`, and those are its +##' coefficients. Its fixed point is `-M/H`, which is what the solving path +##' returns — the two agree in the large-gain limit. +##' +##' ⚠️ **`stationarity` is still reported, and is used for exactly one thing.** +##' It no longer routes anything; it now measures how far the collar you supplied +##' sits from the optimum. The one decision it makes is `profit`'s: at a +##' stationary point the envelope theorem applies and the analytic zero is used, +##' and away from one the exact `dprofit/dpsi` is used instead of differencing +##' it. So `psi = ` with `dpsi_dtheta = -M/H` reproduces the +##' solving path **bit-for-bit**, which the tests assert. +##' +##' ⚠️ **An infeasible `psi` gets no gradient either**, reported as +##' `"no-gradient"`. `dprofit/dpsi` is a sentinel `0.0` on the shut-down and +##' reversed-gradient exits rather than a derivative, and this path never divides +##' by `H`, so it would otherwise adopt that zero as the real thing. +##' +##' ⚠️ **A clamped `psi` gets no gradient, and that is the interesting case.** +##' The collar actually used is `psi` clamped into the feasible interval, so it +##' moves with the *bound* rather than with `dpsi_dtheta` — the active-set +##' problem, arriving through the clamp instead of through the optimiser. The +##' direct term alone would be plausible and wrong. `status` reports `"clamped"` +##' and the gradient is all `NA`; `psi` in the result is the collar that was +##' used. This is reported rather than thrown because a tracking model reaches +##' these points routinely — the clamp is how TF24f pulls an out-of-range state +##' back inside — and it fires for a `psi` within one step of an end too, since +##' `dY/dpsi` cannot be centred there. +##' ##' @section The active set, which is the reason this function is careful: ##' Stationarity is the premise of the whole derivation, and it fails when the ##' optimum is pinned to an end of the feasible collar interval. There `psi*` is a @@ -301,6 +349,16 @@ set_traits <- function(x, traits) { ##' [root_network_from_carbon()] is homogeneous of degree 1 in each constant, so ##' scaling `beta_R_H` scales `r_R_H_min` by the same factor and needs no ##' rebuild — but the two solves are still two solves. +##' @param psi a collar water potential to evaluate at, in MPa as a positive +##' magnitude, instead of solving for the profit-maximising one. `NULL` (the +##' default) solves. See the section above; `method` cannot be given with it. +##' @param dpsi_dtheta how the prescribed `psi` itself responds to each +##' parameter, as one value per entry of `pars` (or one value recycled, or a +##' vector named by `pars`). Defaults to zero, which is the partial derivative +##' at fixed collar. Must be finite — an infinite value would return an +##' all-infinite gradient row with a finite, plausible-looking `profit` in it, +##' because the envelope theorem assigns that column rather than multiplying +##' through. Only meaningful with `psi`. ##' @param step relative step for the trait difference. The default `1e-06` is ##' near the middle of the five decades over which the mixed partial was ##' measured stable; it is also used, relative to the collar potential, for the @@ -331,11 +389,23 @@ set_traits <- function(x, traits) { ##' `profit` (also umol C m^-2 s^-1 per trait unit)} ##' \item{`value`}{the solved outputs the gradient is taken at} ##' \item{`method`}{`"ift"` if the implicit-function composite was used, -##' `"fd"` if the fallback was} -##' \item{`status`}{`"interior"`, `"pinned"` or `"no-gradient"`} -##' \item{`H`}{the curvature of profit in the collar potential at `psi*`} -##' \item{`stationarity`}{the implied Newton step, in MPa, that `method` was -##' decided on} +##' `"fd"` if the fallback was, `"prescribed"` if `psi` was given} +##' \item{`status`}{`"interior"`, `"pinned"` or `"no-gradient"` on the +##' solving path; `"prescribed"`, `"clamped"` or `"no-gradient"` with +##' `psi`} +##' \item{`H`}{the curvature of profit in the collar potential} +##' \item{`stationarity`}{the implied Newton step, in MPa: what `method` was +##' decided on, or with `psi` how far the collar you gave is from the +##' optimum} +##' \item{`M`}{the mixed partial `d2profit/dpsi dtheta`, one per entry of +##' `pars`, or `NA` where the fallback ran} +##' \item{`dY_dpsi`}{`dY/dpsi` at fixed traits, one per output. `profit`'s +##' entry is the analytic zero at a stationary point and the exact +##' `dprofit/dpsi` elsewhere -- see the envelope section} +##' \item{`psi`}{the collar the outputs were evaluated at: `psi*` on the +##' solving path, and on the prescribed path the value actually used, +##' which differs from the `psi` argument exactly when `status` is +##' `"clamped"`} ##' } ##' ##' @seealso [leaf_solve()] for the operating point itself, [set_traits()]. @@ -366,6 +436,8 @@ leaf_gradient <- function(psi_soil, control = leaf_control(), supply = leaf_supply_multilayer(), pars = NULL, + psi = NULL, + dpsi_dtheta = NULL, step = 1e-6, stationarity_tol = 1e-8, method = c("auto", "ift", "fd"), @@ -377,6 +449,12 @@ leaf_gradient <- function(psi_soil, if (!(is.numeric(step) && length(step) == 1L && step > 0)) { stop("`step` must be a single positive number", call. = FALSE) } + # ⚠️ THE RETURN IS THE POINT, NOT JUST THE CHECK. `psi` is compared against + # `opt_root_psi_` to detect the clamp, and a caller writing `psi = 3L` -- or the + # entirely ordinary `for (p in 2:5)`, since `2:5` is integer -- would otherwise + # be reported CLAMPED at a collar it was given exactly. The checker coerces; + # discarding what it returns was the bug. + psi <- .gradient_check_psi(psi, dpsi_dtheta, method) # --- reusing a Leaf (#52) ------------------------------------------------- # `x` is a VESSEL, not the point being differentiated: `traits` still says where @@ -441,6 +519,12 @@ leaf_gradient <- function(psi_soil, # Shared with leaf_gradient_batch(), so the two entry points cannot disagree # about which parameters exist or explain a rejection differently. .gradient_check_pars(pars, identical(supply$kind, "single")) + # ⚠️ HERE, NOT WHERE IT IS USED. `dpsi_dtheta` is only READ inside the + # composite, which a clamped or shut-down operating point never reaches -- so + # validating it there made a wrong-length vector an error at some psi and + # silently accepted at others. Argument checking cannot be conditional on what + # the model does with the arguments. + dpsi_dtheta <- .gradient_dpsi_dtheta(dpsi_dtheta, pars, !is.null(psi)) # ONE leaf for the whole gradient, re-traited rather than reconstructed. This is # the measurement that reordered PLAN 11d: a fresh Leaf costs ~155 us against @@ -473,15 +557,43 @@ leaf_gradient <- function(psi_soil, }, silent = TRUE) }, add = TRUE) } - l$find_root_collar_psi() - - psi_star <- l$opt_root_psi_ + # --- seat the operating point --------------------------------------------- + # Two ways in, and which one ran is what everything below branches on. The + # default SOLVES for the collar potential; `psi` IMPOSES one, which is what a + # caller tracking the optimum rather than finding it has (#88). + # + # ⚠️ `psi_star` keeps its name on both paths and it is no longer always the + # argmax. It is "the collar the outputs were evaluated at", which is what every + # use of it below actually means. + prescribed <- !is.null(psi) + if (prescribed) { + l$evaluate_root_collar_psi(psi) + # Exact equality, for the reason `.gradient_outputs_at` documents: the clamp + # is a min/max, so an unclamped target returns bit-identically. + clamped <- !identical(l$opt_root_psi_, psi) + psi_star <- l$opt_root_psi_ + } else { + l$find_root_collar_psi() + clamped <- FALSE + psi_star <- l$opt_root_psi_ + } value <- .gradient_outputs(l) - # Is the premise true here? See the active-set section: the test is the implied - # Newton step, which is a distance in MPa and so needs no scale of its own. + # The curvature and the residual, at whichever collar was seated. Both are + # differences of `dprofit_droot_collar_psi`, which takes psi as an ARGUMENT -- + # so neither needs the point to be an optimum, and both mean the same thing on + # the two paths. h_psi <- max(abs(psi_star), 1) * step - resid <- l$dprofit_droot_collar_psi(psi_star) + # ⚠️ WITH ITS FEASIBILITY, not bare. `dprofit_droot_collar_psi` returns a hard + # 0.0 SENTINEL on its shut-down and reversed-gradient exits, and a bare zero is + # indistinguishable from a stationary point -- the header says so, and says a + # composite that ignores the flag inherits the bug. The solving path got away + # with reading the value alone because `H` collapses to zero too and `usable` + # catches the pair; the prescribed path does NOT divide by `H`, so it would + # have adopted the sentinel as if it were dprofit/dpsi. + checked <- l$dprofit_droot_collar_psi_checked(psi_star) + resid <- checked[[1L]] + feasible <- checked[[2L]] == 1 H <- (l$dprofit_droot_collar_psi(psi_star + h_psi) - l$dprofit_droot_collar_psi(psi_star - h_psi)) / (2 * h_psi) # H == 0 with resid == 0 is the shut-down signature: dprofit returns a sentinel @@ -489,21 +601,60 @@ leaf_gradient <- function(psi_soil, # would not be a maximum. Both mean the composite has nothing to stand on. usable <- is.finite(H) && H < 0 && is.finite(resid) stationarity <- if (usable) abs(resid / H) else Inf - status <- if (!usable) "no-gradient" else - if (stationarity > stationarity_tol) "pinned" else "interior" + + if (prescribed) { + # ⚠️ THE STATIONARITY TEST DOES NOT ROUTE HERE, AND IS STILL WORTH TAKING. + # What it decides on the solving path -- composite or fallback -- is + # meaningless at a collar the caller chose: there is no argmax to be pinned + # against, and differencing the solve would answer a question about the + # optimum instead of about this point. So `method` is rejected upstream. + # + # The NUMBER keeps its meaning, though, and gains a better one: it is how far + # the point you handed over sits from the optimum, in MPa. It is reported, + # and it is used for exactly one thing -- see `envelope` below. + # + # ⚠️ `no-gradient` REACHES THIS PATH TOO, and it is not the same condition as + # the solving path's. There, `usable` also demands `H < 0` -- a MAXIMUM test, + # which is exactly what a caller-chosen collar has no business satisfying: a + # prescribed psi away from the optimum may sit where profit is convex, and + # that is fine, because nothing here divides by `H`. What genuinely disables + # the point is INFEASIBILITY: the shut-down and reversed-gradient exits, where + # `dprofit` is a sentinel rather than a derivative. Almost every such point is + # already caught as `clamped` -- the shut-down state seats a collar of its own + # choosing -- but "almost" is not a guarantee, since a caller can pass exactly + # that collar back. + status <- if (clamped) "clamped" else if (!feasible) "no-gradient" else + "prescribed" + } else { + status <- if (!usable) "no-gradient" else + if (stationarity > stationarity_tol) "pinned" else "interior" + } # `status` describes the POINT and is reported whichever route runs; `use_ift` # is the route. They differ only when the caller has forced one. - use_ift <- switch(method, + use_ift <- if (prescribed) !clamped && feasible else switch(method, auto = identical(status, "interior"), ift = TRUE, fd = FALSE) - if (use_ift && !usable) { + if (use_ift && !prescribed && !usable) { stop("leaf_gradient(): method = \"ift\" was asked for at a point with no ", "usable curvature (H = ", format(H), "), so -M/H has nothing to stand ", "on. This is a shut-down or otherwise determined operating point; use ", "method = \"auto\".", call. = FALSE) } + # ⚠️ A CLAMPED PRESCRIBED PSI GETS NO GRADIENT, RATHER THAN THE DIRECT TERM. + # It is not a failure -- the outputs at the clamped collar are perfectly good, + # and TF24f relies on the clamp to pull an out-of-range tracked state back + # inside. It is that the derivative is not the one this can compute: the collar + # actually used is `min(max(psi, a(theta)), b(theta))`, so it moves with the + # BOUND, and dY/dtheta picks up the bound's derivative rather than the caller's + # `dpsi_dtheta`. That is the active-set problem arriving through the clamp + # instead of through the optimiser, and the direct term alone would be + # plausible and wrong in exactly the documented way. + # + # Reported rather than thrown, because a fit will visit these points routinely + # and `status` is how a batch tells its caller which rows to distrust. The + # all-NA result is assembled below, with the second way of reaching it. if (use_ift) { # dY/dpsi at fixed traits. evaluate_root_collar_psi CLAMPS its target into the @@ -528,26 +679,37 @@ leaf_gradient <- function(psi_soil, hi <- .gradient_outputs_at(l, psi_star + h_psi) lo <- .gradient_outputs_at(l, psi_star - h_psi) if (is.null(hi) || is.null(lo)) { - if (identical(method, "ift")) { + if (prescribed) { + # Same reasoning as the clamp on `psi` itself, one step out: a psi that + # is inside the interval but within `h_psi` of an end cannot have dY/dpsi + # centred on it, and a one-sided difference over a shortened interval is + # the failure this detector exists for. There is no fallback to offer -- + # differencing the solve would answer about the optimum -- so the row is + # reported as clamped. + clamped <- TRUE + status <- "clamped" + use_ift <- FALSE + } else if (identical(method, "ift")) { stop("leaf_gradient(): method = \"ift\" was asked for at a point whose ", "feasible collar interval is narrower than one step, so dY/dpsi ", "cannot be centred on psi*. Use method = \"auto\".", call. = FALSE) + } else { + use_ift <- FALSE + status <- "pinned" } - use_ift <- FALSE - status <- "pinned" } else { dY_dpsi <- (hi - lo) / (2 * h_psi) # `dprofit_droot_collar_psi` is EXACT in psi -- forward AD plus the IFT at # the ci root-find -- so for profit alone the package has something better # than a difference of the same quantity, and it is already computed. The - # other four have no such route and must be differenced. + # other four have no such route and must be differenced. One rule, both + # paths, which is what keeps `psi = psi*` reproducing the solve exactly. # - # ⚠️ NOTHING CONSUMES THIS TODAY, and it is here rather than deleted for #88. - # The only reader is `.gradient_ift()` with `envelope = FALSE`, i.e. a FORCED - # method = "ift" at a pinned point -- and the 288-point grid test records that - # forcing it there throws ("narrower than one step") at all 42 pinned rows - # before this value is reached. So it is unexercised, not load-bearing: do - # not read a green suite as evidence about it. + # ⚠️ THIS IS THE READER #87 SAID DID NOT EXIST YET. Until the prescribed path + # landed, the only consumer was `.gradient_ift(envelope = FALSE)` -- a forced + # method = "ift" at a pinned point, which throws at all 42 pinned rows of the + # grid. A prescribed psi away from the optimum is not stationary, so it takes + # this branch for real, and the exactness now matters. dY_dpsi[["profit"]] <- resid } } @@ -559,28 +721,50 @@ leaf_gradient <- function(psi_soil, # near-zero dY/dpsi -- the same treatment `collar` gets, and for the same # reason: an identity is stated, not arrived at. # - # ⚠️ It is conditional on stationarity and not on `use_ift`. The identity comes - # from dprofit/dpsi == 0; at a pinned optimum psi* is a theta-dependent BOUND, - # dprofit/dpsi is not zero there, and the indirect term survives. Someone who - # forces method = "ift" at such a point already gets a confidently wrong number - # and should not get a differently wrong one for this column alone. + # ⚠️ The test is `stationarity`, NOT `status`, so the two paths agree at a psi + # that happens to BE psi*: `status` carries "prescribed" on one and "interior" + # on the other, and routing on it would break the bit-for-bit equivalence while + # looking equivalent. It is also not `use_ift`: at a pinned optimum psi* is a + # theta-dependent BOUND, dprofit/dpsi is not zero, and the indirect term + # survives -- someone forcing method = "ift" there already gets a confidently + # wrong number and should not get a differently wrong one for this column. # # The measured size of what this removes is in ?leaf_gradient. envelope <- isTRUE(usable && stationarity <= stationarity_tol) - grad <- if (use_ift) { - .gradient_ift(l, reset, theta, pars, psi_star, H, dY_dpsi, step, - fast_stem_curve, envelope = envelope) + if (prescribed && !use_ift) { + grad <- matrix(NA_real_, length(pars), length(.gradient_output_names()), + dimnames = list(pars, .gradient_output_names())) + M <- stats::setNames(rep(NA_real_, length(pars)), pars) + dY_dpsi <- stats::setNames(rep(NA_real_, length(.gradient_output_names())), + .gradient_output_names()) + } else if (use_ift) { + # ONE composite for both paths, which is what makes the equivalence above a + # real assertion rather than two implementations that happen to agree. The + # only difference is where dpsi/dtheta comes from: derived by the implicit + # function theorem when the collar was solved for, and supplied by the caller + # when it was imposed -- because then it is not this function's to know. + fit <- .gradient_ift(l, reset, theta, pars, psi_star, H, dY_dpsi, step, + fast_stem_curve, dpsi_dtheta = dpsi_dtheta, + envelope = envelope) + grad <- fit$gradient + M <- fit$M } else { - .gradient_fd(l, reset, theta, pars, step, fast_stem_curve) + grad <- .gradient_fd(l, reset, theta, pars, step, fast_stem_curve) + M <- stats::setNames(rep(NA_real_, length(pars)), pars) + dY_dpsi <- stats::setNames(rep(NA_real_, length(.gradient_output_names())), + .gradient_output_names()) } list(gradient = grad, value = value, - method = if (use_ift) "ift" else "fd", + method = if (prescribed) "prescribed" else if (use_ift) "ift" else "fd", status = status, H = H, - stationarity = stationarity) + stationarity = stationarity, + M = M, + dY_dpsi = dY_dpsi, + psi = psi_star) } # --- internals --------------------------------------------------------------- @@ -906,8 +1090,14 @@ leaf_gradient <- function(psi_soil, # neither of which re-solves the model: `dprofit` at the UNPERTURBED psi* gives # the mixed partial, and the outputs at that same psi* give the direct term. .gradient_ift <- function(l, reset, theta, pars, psi_star, H, dY_dpsi, step, - fast_stem_curve = TRUE, envelope = FALSE) { + fast_stem_curve = TRUE, dpsi_dtheta = NULL, + envelope = FALSE) { seat <- .gradient_reseat_base(reset, theta, pars, fast_stem_curve) + # The mixed partials, kept rather than consumed. `-M/H` is what this function + # needs, but `M` and `H` are also what a caller integrating its own sensitivity + # of psi needs (traitecoevo/plant#614), and they are not recoverable from the + # gradient once divided. + M <- numeric(length(pars)) out <- t(vapply(seq_along(pars), function(k) { p <- pars[[k]] seat(k) @@ -929,19 +1119,25 @@ leaf_gradient <- function(psi_soil, up <- side(1) dn <- side(-1) # M = d2profit/dpsi dtheta, with psi held FIXED at psi*. - dpsi_dtheta <- -((up[["dprofit"]] - dn[["dprofit"]]) / (2 * h)) / H + M[[k]] <<- (up[["dprofit"]] - dn[["dprofit"]]) / (2 * h) + # Where the collar was SOLVED for, psi* moves with theta and the implicit + # function theorem says how. Where it was IMPOSED, it does not move unless + # the caller says it does -- so the caller supplies dpsi/dtheta, defaulting + # to zero, and this function has no business deriving one. + d_psi <- if (is.null(dpsi_dtheta)) -M[[k]] / H else dpsi_dtheta[[k]] direct <- (up[-1L] - dn[-1L]) / (2 * h) # `collar` is not an output of the evaluation -- it IS psi*, held fixed, so # its direct term is zero by construction and the composite reduces to # dpsi*/dtheta. Setting it explicitly says so, rather than relying on the # difference of two identical numbers. - g <- direct + dY_dpsi * dpsi_dtheta - g[["collar"]] <- dpsi_dtheta + g <- direct + dY_dpsi * d_psi + g[["collar"]] <- d_psi # The envelope theorem, ASSIGNED for the same reason `collar` is: profit's # indirect term is identically zero at a stationary point, so stating that is - # better than multiplying a measured near-zero by dpsi*/dtheta. It is also - # immune to a non-finite dpsi*/dtheta, where `0 * x` would be NaN in this one - # column while the other four carried +-Inf. + # better than multiplying a measured near-zero by dpsi/dtheta. It is also + # immune to a non-finite `d_psi` -- which a CALLER supplies on the prescribed + # path -- where `0 * x` would be NaN in this one column while the other four + # carried +-Inf. if (envelope) { g[["profit"]] <- direct[["profit"]] } @@ -949,7 +1145,95 @@ leaf_gradient <- function(psi_soil, }, numeric(length(.gradient_output_names())))) reset(theta) rownames(out) <- pars - out + list(gradient = out, M = stats::setNames(M, pars)) +} + +# `psi` and `dpsi_dtheta` only mean anything together, and neither means anything +# alongside `method`. +.gradient_check_psi <- function(psi, dpsi_dtheta, method) { + if (is.null(psi)) { + if (!is.null(dpsi_dtheta)) { + stop("`dpsi_dtheta` is the trait response of a collar potential YOU ", + "imposed, so it needs `psi`. Without one the collar is solved for ", + "and its response is derived by the implicit function theorem.", + call. = FALSE) + } + return(NULL) + } + if (!(is.numeric(psi) && length(psi) == 1L && is.finite(psi) && psi > 0)) { + stop("`psi` must be a single finite positive number: it is a collar water ", + "potential in MPa, as a positive magnitude (see the sign convention ", + "in ?leaf_model).", call. = FALSE) + } + if (!identical(method, "auto")) { + stop("`method` cannot be given with `psi`. The two routes it chooses ", + "between are about a SOLVED operating point -- \"fd\" differences the ", + "solve, which would answer about the optimum rather than about the ", + "collar you gave -- so at a prescribed psi neither applies and the ", + "reported method is \"prescribed\".", call. = FALSE) + } + # ⚠️ COERCED, AND THE CALLER MUST USE THE RETURN. `psi` is compared against + # `opt_root_psi_` with `identical()` to detect the clamp, and `identical(3, 3L)` + # is FALSE -- so an integer `psi` was reported clamped at a collar it had been + # given exactly, with the "where it was pulled to" diagnostic showing no + # movement and nothing to notice. `.gradient_check_psi_batch()` already returned + # `as.numeric(psi)`, and C++'s `util::identical` is `a == b`, so R alone had it. + as.numeric(psi) +} + +# dpsi/dtheta as one value per parameter, in `pars` order. NULL on the solving +# path means "derive it"; NULL on the prescribed path means "the collar does not +# move with theta", which is zero and not the same statement. +.gradient_dpsi_dtheta <- function(dpsi_dtheta, pars, prescribed) { + if (!prescribed) { + return(NULL) + } + if (is.null(dpsi_dtheta)) { + return(stats::setNames(numeric(length(pars)), pars)) + } + # ⚠️ FINITE, AND `anyNA` WAS NOT ENOUGH -- `is.finite` covers NA and NaN, + # so this is one check rather than two, and it adds the case that matters. An + # infinite `dpsi_dtheta` does not fail loudly. The composite is + # `direct + dY_dpsi * dpsi_dtheta`, so four columns come back +-Inf -- and + # `profit`, which the envelope theorem ASSIGNS from the direct term at a + # stationary psi, comes back FINITE AND PLAUSIBLE beside them. Measured at + # `psi = psi*` with `dpsi_dtheta = Inf`: A, gc, psi_stem and collar all Inf, + # profit 0.0105. Reading `profit` alone is plant's own case (#87), so the one + # column that survives is the one most likely to be believed. The sibling + # argument is already checked this way -- `.gradient_check_psi()` demands + # `is.finite(psi)`. + if (!is.numeric(dpsi_dtheta) || !all(is.finite(dpsi_dtheta))) { + stop("`dpsi_dtheta` must be numeric and finite -- no NA, NaN or Inf. It ", + "is dpsi/dtheta for a collar YOU imposed; a non-finite value returns ", + "an all-infinite gradient row with a plausible `profit` in it.", + call. = FALSE) + } + # ⚠️ NAMES ARE CHECKED BEFORE RECYCLING, and the order matters. Recycling first + # overwrote whatever the caller wrote with `pars`, so `dpsi_dtheta = + # c(stem_b = 1)` against `pars = "vcmax_25"` was silently applied to vcmax_25 -- + # a named argument quietly meaning a different parameter, which is the exact + # failure naming it was supposed to prevent. + if (!is.null(names(dpsi_dtheta))) { + # Matched, not assumed aligned: reordering it silently would be the same + # class of bug as reordering `pars`. + if (!setequal(names(dpsi_dtheta), pars)) { + stop("`dpsi_dtheta` is named, so its names must be exactly `pars`. ", + "Missing: ", paste(setdiff(pars, names(dpsi_dtheta)), + collapse = ", "), + "; unexpected: ", paste(setdiff(names(dpsi_dtheta), pars), + collapse = ", "), call. = FALSE) + } + return(dpsi_dtheta[pars]) + } + if (length(dpsi_dtheta) == 1L) { + return(stats::setNames(rep(dpsi_dtheta, length(pars)), pars)) + } + if (length(dpsi_dtheta) != length(pars)) { + stop("`dpsi_dtheta` must be length 1 or one value per parameter (", + length(pars), " for this `pars`); got ", length(dpsi_dtheta), + call. = FALSE) + } + stats::setNames(dpsi_dtheta, pars) } # The fallback: a central difference of the whole solve. Correct at a pinned diff --git a/README.md b/README.md index e98d9f8..613c3ea 100644 --- a/README.md +++ b/README.md @@ -355,6 +355,21 @@ exception: it *is* the objective, so at an interior optimum the second route contributes nothing and its gradient is the direct partial alone. That is the envelope theorem, and the only place this package uses it. +If your model **tracks** the optimum instead of finding it, pass the collar +potential you are operating at and the derivation simplifies rather than breaks: + +```r +leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = c("vcmax_25", "stem_b"), + psi = 2.7) # evaluate here, do not solve +``` + +`psi` is exogenous, so the answer is the partial at fixed collar — plus whatever +`dpsi_dtheta` you supply, if the collar you imposed itself moves with the traits. +`M`, `H` and `dY_dpsi` come back so a caller integrating its own sensitivity of +`psi` has the coefficients. This is what plant's TF24f needs, and giving back the +collar the solver found, with the response it derived, reproduces the solving +path exactly. + That derivation assumes the optimum is interior, and at the dry end it often is not: with the optimum pinned to the edge of the feasible range the formula returns a confidently wrong number, off by up to seven orders of magnitude. So the diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index 4af1b02..969aa97 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -391,6 +391,9 @@ Leaf: dprofit_droot_collar_psi: return_type: double args: [opt_root_psi: double] + dprofit_droot_collar_psi_checked: + return_type: "std::vector" + args: [opt_root_psi: double] psi_stem_to_ci: return_type: double args: [psi_stem: double, psi_upstream: double] diff --git a/inst/include/phylloptim/gradient.hpp b/inst/include/phylloptim/gradient.hpp index 8abbdc0..3270ac2 100644 --- a/inst/include/phylloptim/gradient.hpp +++ b/inst/include/phylloptim/gradient.hpp @@ -224,15 +224,29 @@ enum class Method { Auto, Ift, Fd }; // What the OPERATING POINT is, not whether the call worked -- except for `Error`, // which is the batch's per-row failure and is the reason this is a status rather // than an exception (see `batch` below). -enum class Status { Interior, Pinned, NoGradient, Error }; - +// `Prescribed` and `Clamped` are the two the caller-supplied-psi path can end in +// (#88): the collar was imposed, and in the second case the feasible interval +// moved it, so the derivative belongs to the BOUND and is not reported. +// ⚠️ APPEND AFTER `Error`, DO NOT INSERT BEFORE IT. The names are +// source-compatible either way, but the VALUES are not: a consumer holding the +// integer would break quietly if `Error` moved from 3 to 5. `Prescribed` and +// `Clamped` are therefore last, out of narrative order, on purpose. +enum class Status { Interior, Pinned, NoGradient, Error, Prescribed, Clamped }; + +// ⚠️ NO `default:`. An exhaustive switch turns the next member added here into a +// compiler diagnostic; a `default: return "error"` would silently LABEL it +// "error" instead, which is exactly the kind of quiet mislabelling `Prescribed` +// and `Clamped` were added to avoid. inline std::string status_name(Status s) { switch (s) { case Status::Interior: return "interior"; case Status::Pinned: return "pinned"; case Status::NoGradient: return "no-gradient"; - default: return "error"; + case Status::Prescribed: return "prescribed"; + case Status::Clamped: return "clamped"; + case Status::Error: return "error"; } + return "error"; // unreachable; silences -Wreturn-type on a bad cast } struct Settings { @@ -242,6 +256,18 @@ struct Settings { bool fast_stem_curve = true; }; +// A collar potential the caller imposes, in place of the one `at` would solve +// for, plus how that collar responds to each parameter (#88). +// +// ⚠️ `dpsi_dtheta` is npars long and in `pars` order, NOT n_pars: it is one value +// per parameter ASKED FOR, because that is the vector a caller integrating its +// own sensitivity state is carrying. Null means zero -- the partial at fixed +// collar -- which is a different statement from the solving path's "derive it". +struct Prescribed { + double psi = 0.0; + const double* dpsi_dtheta = nullptr; +}; + struct Result { // The five outputs the gradient is taken at. double value[n_outputs]; @@ -252,17 +278,30 @@ struct Result { bool used_ift = false; double H = util::na_value; double stationarity = util::na_value; + // The mixed partials d2profit/dpsi dtheta, npars long, and dY/dpsi at fixed + // traits. Kept rather than consumed: `-M/H` is what the composite needs, but M + // and H are also the coefficients of a caller's own sensitivity ODE for psi + // (traitecoevo/plant#614), and they cannot be recovered once divided. + std::vector M; + double dY_dpsi[n_outputs]; + // The collar the outputs were evaluated at -- psi* when solved, and when + // prescribed the value actually USED, which differs from the requested one + // exactly when `status` is Clamped. + double psi = util::na_value; std::string message; void reset(std::size_t npars) { for (int j = 0; j < n_outputs; ++j) { value[j] = util::na_value; + dY_dpsi[j] = util::na_value; } grad.assign(npars * n_outputs, util::na_value); + M.assign(npars, util::na_value); status = Status::Error; used_ift = false; H = util::na_value; stationarity = util::na_value; + psi = util::na_value; message.clear(); } }; @@ -350,10 +389,17 @@ inline bool takes_shortcut(int par, const Settings& s) { // The implicit-function composite. Two perturbed evaluations per parameter, // neither of which re-solves the model: `dprofit` at the UNPERTURBED psi* gives // the mixed partial, and the outputs at that same psi* give the direct term. +// +// ONE composite for both paths. `dpsi_dtheta` null means "derive it by the +// implicit function theorem", which is right where the collar was solved for; +// non-null means the caller imposed the collar and knows how it moves. Keeping +// them in one function is what makes the two paths' agreement at psi* a real +// assertion rather than two implementations that happen to line up. inline void gradient_ift(Leaf& l, const double* theta, const Drivers& d, bool single, const int* pars, std::size_t npars, double psi_star, double H, const double* dY_dpsi, - const Settings& s, bool envelope, double* out) { + const Settings& s, const double* dpsi_dtheta, + bool envelope, double* M_out, double* out) { double th[n_pars]; double up[1 + n_outputs]; double dn[1 + n_outputs]; @@ -384,21 +430,24 @@ inline void gradient_ift(Leaf& l, const double* theta, const Drivers& d, dst[0] = l.dprofit_droot_collar_psi(psi_star); } // M = d2profit/dpsi dtheta, with psi held FIXED at psi*. - const double dpsi_dtheta = -((up[0] - dn[0]) / (2.0 * h)) / H; + const double M = (up[0] - dn[0]) / (2.0 * h); + M_out[k] = M; + const double d_psi = dpsi_dtheta == nullptr ? -M / H : dpsi_dtheta[k]; for (int j = 0; j < n_outputs; ++j) { const double direct = (up[1 + j] - dn[1 + j]) / (2.0 * h); - out[k * n_outputs + j] = direct + rounded(dY_dpsi[j] * dpsi_dtheta); + out[k * n_outputs + j] = direct + rounded(dY_dpsi[j] * d_psi); } // `collar` is not an output of the evaluation -- it IS psi*, held fixed, so // its direct term is zero by construction and the composite reduces to // dpsi*/dtheta. Set explicitly rather than left as the difference of two // identical numbers. - out[k * n_outputs + out_collar] = dpsi_dtheta; + out[k * n_outputs + out_collar] = d_psi; // The envelope theorem, ASSIGNED for the same reason `collar` is: profit's // indirect term is identically zero at a stationary point, so stating that - // beats multiplying a measured near-zero by dpsi*/dtheta. It is also immune - // to a non-finite dpsi*/dtheta, where `0 * x` would be NaN in this one column - // while the other four carried +-Inf. + // beats multiplying a measured near-zero by dpsi/dtheta. It is also immune + // to a non-finite `d_psi` -- which a CALLER supplies on the prescribed path + // -- where `0 * x` would be NaN in this one column while the other four + // carried +-Inf. if (envelope) { out[k * n_outputs + out_profit] = (up[1 + out_profit] - dn[1 + out_profit]) / (2.0 * h); @@ -445,13 +494,28 @@ inline void gradient_fd(Leaf& l, const double* theta, const Drivers& d, // when `batch` catches. See `batch` for why the batch does not propagate. inline void at(Leaf& l, const double* theta, const Drivers& d, bool single, const int* pars, std::size_t npars, const Settings& s, - Result& out) { + Result& out, const Prescribed* prescribed = nullptr) { out.reset(npars); apply(l, theta, d, single, -1, s.fast_stem_curve); - l.find_root_collar_psi(); + // Two ways in. The default SOLVES for the collar potential; `prescribed` + // IMPOSES one, which is what a caller tracking the optimum rather than + // finding it has (#88). + // + // ⚠️ `psi_star` keeps its name on both paths and is no longer always the + // argmax. It is "the collar the outputs were evaluated at", which is what + // every use of it below actually means. + bool clamped = false; + if (prescribed != nullptr) { + l.evaluate_root_collar_psi(prescribed->psi); + // Exact equality, for the reason `outputs_at` documents above. + clamped = !util::identical(l.opt_root_psi_, prescribed->psi); + } else { + l.find_root_collar_psi(); + } const double psi_star = l.opt_root_psi_; + out.psi = psi_star; outputs(l, out.value); // Is the composite's premise true HERE? Stationarity is what the whole @@ -466,7 +530,14 @@ inline void at(Leaf& l, const double* theta, const Drivers& d, bool single, // mildest pinned one 6.3e-06, so the 1e-08 default sits in an empty band four // orders wide on each side. const double h_psi = std::max(std::abs(psi_star), 1.0) * s.step; - const double resid = l.dprofit_droot_collar_psi(psi_star); + // ⚠️ WITH ITS FEASIBILITY, not bare. `dprofit_droot_collar_psi` returns a hard + // 0.0 SENTINEL on its shut-down and reversed-gradient exits, and a bare zero is + // indistinguishable from a stationary point. The solving path got away with the + // value alone because `H` collapses to zero too and `usable` catches the pair; + // the prescribed path does not divide by `H`, so it would adopt the sentinel as + // if it were dprofit/dpsi. + bool resid_feasible = false; + const double resid = l.dprofit_droot_collar_psi(psi_star, &resid_feasible); // Named halves: `f(a) - f(b)` has unspecified operand order in C++ and // left-to-right order in R, and `dprofit_droot_collar_psi` mutates the leaf. const double d_hi = l.dprofit_droot_collar_psi(psi_star + h_psi); @@ -480,15 +551,41 @@ inline void at(Leaf& l, const double* theta, const Drivers& d, bool single, out.H = H; out.stationarity = usable ? std::abs(resid / H) : std::numeric_limits::infinity(); - out.status = !usable ? Status::NoGradient - : (out.stationarity > s.stationarity_tol ? Status::Pinned - : Status::Interior); + if (prescribed != nullptr) { + // ⚠️ THE STATIONARITY TEST DOES NOT ROUTE HERE, AND IS STILL WORTH TAKING. + // What it decides on the solving path -- composite or fallback -- is + // meaningless at a collar the caller chose: there is no argmax to be pinned + // against, and differencing the solve would answer about the optimum + // instead. So `Method` is refused at the R boundary. + // + // The NUMBER keeps its meaning and gains a better one: it is how far the + // point handed over sits from the optimum, in MPa. Reported, and used for + // exactly one thing -- see the envelope below. + // + // ⚠️ NoGradient REACHES THIS PATH TOO, and it is not the solving path's + // condition. There `usable` also demands `H < 0` -- a MAXIMUM test, which a + // caller-chosen collar has no business satisfying: a prescribed psi away from + // the optimum may sit where profit is convex, and that is fine because + // nothing here divides by `H`. What disables the point is INFEASIBILITY. + // Almost every such point is already Clamped -- the shut-down state seats a + // collar of its own choosing -- but a caller can pass exactly that collar + // back, so "almost" is not a guarantee. + out.status = clamped ? Status::Clamped + : (resid_feasible ? Status::Prescribed : Status::NoGradient); + } else { + out.status = !usable ? Status::NoGradient + : (out.stationarity > s.stationarity_tol ? Status::Pinned + : Status::Interior); + } // `status` describes the POINT and is reported whichever route runs; // `use_ift` is the route. They differ only when the caller has forced one. - bool use_ift = s.method == Method::Auto ? out.status == Status::Interior - : s.method == Method::Ift; - if (use_ift && !usable) { + bool use_ift = prescribed != nullptr + ? (!clamped && resid_feasible) + : (s.method == Method::Auto + ? out.status == Status::Interior + : s.method == Method::Ift); + if (use_ift && prescribed == nullptr && !usable) { util::stop("leaf_gradient(): method = \"ift\" was asked for at a point with " "no usable curvature (H = " + util::to_string(H) + "), so -M/H " "has nothing to stand on. This is a shut-down or otherwise " @@ -514,14 +611,25 @@ inline void at(Leaf& l, const double* theta, const Drivers& d, bool single, const bool hi_ok = outputs_at(l, psi_star + h_psi, hi); const bool lo_ok = outputs_at(l, psi_star - h_psi, lo); if (!hi_ok || !lo_ok) { - if (s.method == Method::Ift) { + if (prescribed != nullptr) { + // The same reasoning as the clamp on psi itself, one step out: a psi + // inside the interval but within h_psi of an end cannot have dY/dpsi + // centred on it, and a one-sided difference over a shortened interval is + // exactly what this detector exists to refuse. No fallback to offer -- + // differencing the solve would answer about the optimum -- so the row is + // reported clamped. + clamped = true; + use_ift = false; + out.status = Status::Clamped; + } else if (s.method == Method::Ift) { util::stop("leaf_gradient(): method = \"ift\" was asked for at a point " "whose feasible collar interval is narrower than one step, " "so dY/dpsi cannot be centred on psi*. Use " "method = \"auto\"."); + } else { + use_ift = false; + out.status = Status::Pinned; } - use_ift = false; - out.status = Status::Pinned; } else { for (int j = 0; j < n_outputs; ++j) { dY_dpsi[j] = (hi[j] - lo[j]) / (2.0 * h_psi); @@ -529,14 +637,14 @@ inline void at(Leaf& l, const double* theta, const Drivers& d, bool single, // `dprofit_droot_collar_psi` is EXACT in psi -- forward AD plus the IFT at // the ci root-find -- so for profit alone there is something better than a // difference of the same quantity, and it is already computed. The other - // four have no such route and must be differenced. + // four have no such route and must be differenced. One rule, both paths, + // which is what keeps a prescribed psi* reproducing the solve exactly. // - // ⚠️ NOTHING CONSUMES THIS TODAY, and it is here rather than deleted for - // #88. The only reader is `gradient_ift` with `envelope` false, i.e. a - // FORCED Method::Ift at a pinned point -- and the 288-point grid test - // records that forcing it there throws ("narrower than one step") at all - // 42 pinned rows before this value is reached. So it is unexercised, not - // load-bearing: do not read a green suite as evidence about it. + // ⚠️ THIS IS THE READER #87 SAID DID NOT EXIST YET. Until the prescribed + // path landed, the only consumer was `gradient_ift` with `envelope` false -- + // a forced Method::Ift at a pinned point, which throws at all 42 pinned rows + // of the grid. A prescribed psi away from the optimum is not stationary, so + // it takes this branch for real, and the exactness now matters. dY_dpsi[out_profit] = resid; } } @@ -557,8 +665,44 @@ inline void at(Leaf& l, const double* theta, const Drivers& d, bool single, out.used_ift = use_ift; if (use_ift) { + for (int j = 0; j < n_outputs; ++j) { + out.dY_dpsi[j] = dY_dpsi[j]; + } + // ⚠️ NULL MEANS TWO DIFFERENT THINGS AND ONLY ONE OF THEM REACHES + // `gradient_ift`. There, null is "derive dpsi/dtheta by the implicit + // function theorem", which is right for a collar that was SOLVED for. A + // prescribed collar with no dpsi_dtheta means the opposite -- it does not + // move with theta, so the answer is the partial at fixed collar -- and + // passing the null straight through would silently give the caller the + // solving path's indirect term on top of it. Zeros are materialised here so + // that every C++ caller gets that, not just `batch`. + std::vector zeros; + const double* d_psi = nullptr; + if (prescribed != nullptr) { + d_psi = prescribed->dpsi_dtheta; + if (d_psi == nullptr) { + zeros.assign(npars, 0.0); + d_psi = zeros.data(); + } + } gradient_ift(l, theta, d, single, pars, npars, psi_star, H, dY_dpsi, s, - envelope, out.grad.data()); + d_psi, envelope, out.M.data(), out.grad.data()); + } else if (prescribed != nullptr) { + // ⚠️ A CLAMPED PRESCRIBED PSI GETS NO GRADIENT, RATHER THAN THE DIRECT TERM. + // Not a failure -- the outputs at the clamped collar are perfectly good, and + // TF24f relies on the clamp to pull an out-of-range tracked state back + // inside. It is that the derivative is not the one this can compute: the + // collar used is min(max(psi, a(theta)), b(theta)), so it moves with the + // BOUND, and dY/dtheta picks up the bound's derivative rather than the + // caller's dpsi_dtheta. That is the active-set problem arriving through the + // clamp instead of through the optimiser, and the direct term alone would be + // plausible and wrong in the documented way. + // + // `reset` already left grad, M and dY_dpsi as NA; value, H, stationarity and + // psi stay, because they describe the point rather than the derivative. + // + // The same holds for an INFEASIBLE one, where `dprofit` is a sentinel rather + // than a derivative -- see the status block above. } else { gradient_fd(l, theta, d, single, pars, npars, s, out.grad.data()); } @@ -584,21 +728,39 @@ inline void at(Leaf& l, const double* theta, const Drivers& d, bool single, // checks the status per row and not per cell would be reading a gradient with a // hole in it. Its `value`, `H` and `stationarity` are kept where they were // determined, because those describe the point rather than the derivative. +// +// `psi` is null to solve, or one collar potential per observation to impose +// (#88); `dpsi_dtheta` is then a COLUMN-MAJOR n x npars matrix, again as R hands +// one over, and null means zero throughout. inline std::vector batch(Leaf& l, const double* theta, std::size_t theta_nrow, const std::vector& drivers, bool single, const int* pars, - std::size_t npars, const Settings& s) { + std::size_t npars, const Settings& s, + const double* psi = nullptr, + const double* dpsi_dtheta = nullptr) { const std::size_t n = drivers.size(); std::vector out(n); double th[n_pars]; + std::vector dpsi(npars); for (std::size_t i = 0; i < n; ++i) { const std::size_t row = theta_nrow == 1 ? 0 : i; for (int j = 0; j < n_pars; ++j) { th[j] = theta[row + std::size_t(j) * theta_nrow]; } + Prescribed p; + if (psi != nullptr) { + p.psi = psi[i]; + if (dpsi_dtheta != nullptr) { + for (std::size_t k = 0; k < npars; ++k) { + dpsi[k] = dpsi_dtheta[i + k * n]; + } + p.dpsi_dtheta = dpsi.data(); + } + } try { - at(l, th, drivers[i], single, pars, npars, s, out[i]); + at(l, th, drivers[i], single, pars, npars, s, out[i], + psi == nullptr ? nullptr : &p); } catch (const std::exception& e) { out[i].status = Status::Error; out[i].used_ift = false; diff --git a/inst/include/phylloptim/leaf_model.hpp b/inst/include/phylloptim/leaf_model.hpp index 8c7eae0..c3e2e05 100644 --- a/inst/include/phylloptim/leaf_model.hpp +++ b/inst/include/phylloptim/leaf_model.hpp @@ -727,6 +727,12 @@ class Leaf { // indistinguishable; and defaulting to nullptr leaves every existing call site // and the generated R binding untouched. double dprofit_droot_collar_psi(double opt_root_psi, bool* feasible = nullptr); + // The same, with the flag as a return rather than an out-parameter, because + // R cannot reach a `bool*` -- RcppR6 has no form for one, so the generated + // binding drops it and every R-side composite silently gets the "a composite + // that ignores it inherits the bug" case this header warns about above. + // `{dprofit, feasible}`, in that order. + std::vector dprofit_droot_collar_psi_checked(double opt_root_psi); // Post-prepare body of dprofit_droot_collar_psi, with the same `feasible` // contract. Assumes the supply path's per-solve caches are already seated, so // the collar solve can share ONE supply_begin_solve across all ~10 of its @@ -2101,6 +2107,13 @@ inline double Leaf::dprofit_droot_collar_psi(double opt_root_psi, bool* feasible return dprofit_at_collar_psi(opt_root_psi, feasible); } +inline std::vector Leaf::dprofit_droot_collar_psi_checked( + double opt_root_psi) { + bool feasible = false; + const double d = dprofit_droot_collar_psi(opt_root_psi, &feasible); + return {d, feasible ? 1.0 : 0.0}; +} + inline double Leaf::dprofit_at_collar_psi(double opt_root_psi, bool* feasible) { using AD = xad::fwd::active_type; const double psi = opt_root_psi; diff --git a/man/leaf_gradient.Rd b/man/leaf_gradient.Rd index 65f8542..793d5f3 100644 --- a/man/leaf_gradient.Rd +++ b/man/leaf_gradient.Rd @@ -20,6 +20,8 @@ leaf_gradient( control = leaf_control(), supply = leaf_supply_multilayer(), pars = NULL, + psi = NULL, + dpsi_dtheta = NULL, step = 1e-06, stationarity_tol = 1e-08, method = c("auto", "ift", "fd"), @@ -86,6 +88,18 @@ differentiate at} scaling `beta_R_H` scales `r_R_H_min` by the same factor and needs no rebuild — but the two solves are still two solves.} +\item{psi}{a collar water potential to evaluate at, in MPa as a positive +magnitude, instead of solving for the profit-maximising one. `NULL` (the +default) solves. See the section above; `method` cannot be given with it.} + +\item{dpsi_dtheta}{how the prescribed `psi` itself responds to each +parameter, as one value per entry of `pars` (or one value recycled, or a +vector named by `pars`). Defaults to zero, which is the partial derivative +at fixed collar. Must be finite — an infinite value would return an +all-infinite gradient row with a finite, plausible-looking `profit` in it, +because the envelope theorem assigns that column rather than multiplying +through. Only meaningful with `psi`.} + \item{step}{relative step for the trait difference. The default `1e-06` is near the middle of the five decades over which the mixed partial was measured stable; it is also used, relative to the collar potential, for the @@ -120,11 +134,23 @@ A list with `profit` (also umol C m^-2 s^-1 per trait unit)} \item{`value`}{the solved outputs the gradient is taken at} \item{`method`}{`"ift"` if the implicit-function composite was used, - `"fd"` if the fallback was} - \item{`status`}{`"interior"`, `"pinned"` or `"no-gradient"`} - \item{`H`}{the curvature of profit in the collar potential at `psi*`} - \item{`stationarity`}{the implied Newton step, in MPa, that `method` was - decided on} + `"fd"` if the fallback was, `"prescribed"` if `psi` was given} + \item{`status`}{`"interior"`, `"pinned"` or `"no-gradient"` on the + solving path; `"prescribed"`, `"clamped"` or `"no-gradient"` with + `psi`} + \item{`H`}{the curvature of profit in the collar potential} + \item{`stationarity`}{the implied Newton step, in MPa: what `method` was + decided on, or with `psi` how far the collar you gave is from the + optimum} + \item{`M`}{the mixed partial `d2profit/dpsi dtheta`, one per entry of + `pars`, or `NA` where the fallback ran} + \item{`dY_dpsi`}{`dY/dpsi` at fixed traits, one per output. `profit`'s + entry is the analytic zero at a stationary point and the exact + `dprofit/dpsi` elsewhere -- see the envelope section} + \item{`psi`}{the collar the outputs were evaluated at: `psi*` on the + solving path, and on the prescribed path the value actually used, + which differs from the `psi` argument exactly when `status` is + `"clamped"`} } } \description{ @@ -195,6 +221,56 @@ indirect term is real. This applies the identity only where `status` is fallback as the other four. } +\section{A collar potential you supply, instead of the one this solves for}{ + +Pass `psi` and the outputs are evaluated at that collar rather than at the +argmax. This is for a caller whose model **tracks** the optimum instead of +finding it — plant's TF24f carries the collar potential as an ODE state, +`dpsi/dt = k * dprofit/dpsi`, so at finite gain its `dprofit/dpsi` is +deliberately non-zero and *is* the acclimation rate. + +There the derivation above simplifies rather than breaks. `psi` is exogenous, +so the indirect term is whatever the caller says it is: + +\deqn{dY/d\theta = \partial Y/\partial\theta|_\psi + + (\partial Y/\partial\psi)(d\psi/d\theta)} + +with `dpsi_dtheta` supplied — defaulting to zero, the partial at fixed +collar. Nothing is derived from `-M/H`, so nothing needs stationarity, and +`method` is refused because the two routes it chooses between are both about +a solved optimum. + +`M`, `H` and `dY_dpsi` come back in the result, because a caller with a +*dynamic* `psi` cannot supply `dpsi_dtheta` as a constant: for the +gradient-ascent law above it obeys `ds/dt = k(M + H s)`, and those are its +coefficients. Its fixed point is `-M/H`, which is what the solving path +returns — the two agree in the large-gain limit. + +⚠️ **`stationarity` is still reported, and is used for exactly one thing.** +It no longer routes anything; it now measures how far the collar you supplied +sits from the optimum. The one decision it makes is `profit`'s: at a +stationary point the envelope theorem applies and the analytic zero is used, +and away from one the exact `dprofit/dpsi` is used instead of differencing +it. So `psi = ` with `dpsi_dtheta = -M/H` reproduces the +solving path **bit-for-bit**, which the tests assert. + +⚠️ **An infeasible `psi` gets no gradient either**, reported as +`"no-gradient"`. `dprofit/dpsi` is a sentinel `0.0` on the shut-down and +reversed-gradient exits rather than a derivative, and this path never divides +by `H`, so it would otherwise adopt that zero as the real thing. + +⚠️ **A clamped `psi` gets no gradient, and that is the interesting case.** +The collar actually used is `psi` clamped into the feasible interval, so it +moves with the *bound* rather than with `dpsi_dtheta` — the active-set +problem, arriving through the clamp instead of through the optimiser. The +direct term alone would be plausible and wrong. `status` reports `"clamped"` +and the gradient is all `NA`; `psi` in the result is the collar that was +used. This is reported rather than thrown because a tracking model reaches +these points routinely — the clamp is how TF24f pulls an out-of-range state +back inside — and it fires for a `psi` within one step of an end too, since +`dY/dpsi` cannot be centred there. +} + \section{The active set, which is the reason this function is careful}{ Stationarity is the premise of the whole derivation, and it fails when the diff --git a/man/leaf_gradient_batch.Rd b/man/leaf_gradient_batch.Rd index 5b662a6..1d5f0d7 100644 --- a/man/leaf_gradient_batch.Rd +++ b/man/leaf_gradient_batch.Rd @@ -9,6 +9,8 @@ leaf_gradient_batch( traits = NULL, pars = NULL, theta = NULL, + psi = NULL, + dpsi_dtheta = NULL, step = 1e-06, stationarity_tol = 1e-08, method = c("auto", "ift", "fd"), @@ -29,6 +31,18 @@ the most expensive thing you can ask for.} \item{theta}{the model parameters as a matrix, for the case where they vary by observation. See the section above.} +\item{psi}{one collar water potential per observation to evaluate at +(positive magnitudes, MPa) instead of solving for the profit-maximising +one. See [leaf_gradient()]'s section on it. Recycling a single value across +the batch is deliberately not offered — rows sharing an imposed collar is a +coincidence rather than a design, and broadcasting one would make a length +mismatch look intentional.} + +\item{dpsi_dtheta}{how each prescribed `psi` responds to each parameter: +an `n` × `length(pars)` matrix, or one vector of `length(pars)` shared by +every observation. Defaults to zero. Must be finite, for the reason +[leaf_gradient()] records. Only meaningful with `psi`.} + \item{step}{relative step for the trait difference. The default `1e-06` is near the middle of the five decades over which the mixed partial was measured stable; it is also used, relative to the collar potential, for the @@ -61,10 +75,15 @@ A list with \item{`gradient`}{an `n` × `length(pars)` × 5 array, the last dimension being `A`, `gc`, `psi_stem`, `collar` and `profit`} \item{`value`}{an `n` × 5 matrix of the solved outputs} - \item{`method`}{`"ift"` or `"fd"` per observation} - \item{`status`}{`"interior"`, `"pinned"`, `"no-gradient"` or `"error"`} + \item{`method`}{`"ift"`, `"fd"` or `"prescribed"` per observation} + \item{`status`}{`"interior"`, `"pinned"`, `"no-gradient"`, `"error"`, + or with `psi` one of `"prescribed"`, `"clamped"`, `"no-gradient"`} \item{`H`, `stationarity`}{the curvature and the implied Newton step `method` was decided on, per observation} + \item{`M`}{an `n` × `length(pars)` matrix of mixed partials + `d2profit/dpsi dtheta`} + \item{`dY_dpsi`}{an `n` × 5 matrix of `dY/dpsi` at fixed traits} + \item{`psi`}{the collar each row was evaluated at} \item{`message`}{`""`, or why that row failed} } } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index de15051..07c576b 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -366,6 +366,18 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// Leaf__dprofit_droot_collar_psi_checked +std::vector Leaf__dprofit_droot_collar_psi_checked(phylloptim::RcppR6::RcppR6 obj_, double opt_root_psi); +RcppExport SEXP _phylloptim_Leaf__dprofit_droot_collar_psi_checked(SEXP obj_SEXP, SEXP opt_root_psiSEXP) { +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 opt_root_psi(opt_root_psiSEXP); + rcpp_result_gen = Rcpp::wrap(Leaf__dprofit_droot_collar_psi_checked(obj_, opt_root_psi)); + return rcpp_result_gen; +END_RCPP +} // Leaf__psi_stem_to_ci double Leaf__psi_stem_to_ci(phylloptim::RcppR6::RcppR6 obj_, double psi_stem, double psi_upstream); RcppExport SEXP _phylloptim_Leaf__psi_stem_to_ci(SEXP obj_SEXP, SEXP psi_stemSEXP, SEXP psi_upstreamSEXP) { @@ -2215,8 +2227,8 @@ BEGIN_RCPP END_RCPP } // gradient_batch_run -Rcpp::List gradient_batch_run(phylloptim::RcppR6::RcppR6 obj_, SEXP drivers, Rcpp::NumericMatrix theta, Rcpp::IntegerVector pars, double step, double stationarity_tol, std::string method, bool fast_stem_curve); -RcppExport SEXP _phylloptim_gradient_batch_run(SEXP obj_SEXP, SEXP driversSEXP, SEXP thetaSEXP, SEXP parsSEXP, SEXP stepSEXP, SEXP stationarity_tolSEXP, SEXP methodSEXP, SEXP fast_stem_curveSEXP) { +Rcpp::List gradient_batch_run(phylloptim::RcppR6::RcppR6 obj_, SEXP drivers, Rcpp::NumericMatrix theta, Rcpp::IntegerVector pars, double step, double stationarity_tol, std::string method, bool fast_stem_curve, SEXP psi, SEXP dpsi_dtheta); +RcppExport SEXP _phylloptim_gradient_batch_run(SEXP obj_SEXP, SEXP driversSEXP, SEXP thetaSEXP, SEXP parsSEXP, SEXP stepSEXP, SEXP stationarity_tolSEXP, SEXP methodSEXP, SEXP fast_stem_curveSEXP, SEXP psiSEXP, SEXP dpsi_dthetaSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -2228,7 +2240,9 @@ BEGIN_RCPP Rcpp::traits::input_parameter< double >::type stationarity_tol(stationarity_tolSEXP); Rcpp::traits::input_parameter< std::string >::type method(methodSEXP); Rcpp::traits::input_parameter< bool >::type fast_stem_curve(fast_stem_curveSEXP); - rcpp_result_gen = Rcpp::wrap(gradient_batch_run(obj_, drivers, theta, pars, step, stationarity_tol, method, fast_stem_curve)); + Rcpp::traits::input_parameter< SEXP >::type psi(psiSEXP); + Rcpp::traits::input_parameter< SEXP >::type dpsi_dtheta(dpsi_dthetaSEXP); + rcpp_result_gen = Rcpp::wrap(gradient_batch_run(obj_, drivers, theta, pars, step, stationarity_tol, method, fast_stem_curve, psi, dpsi_dtheta)); return rcpp_result_gen; END_RCPP } @@ -2274,6 +2288,7 @@ static const R_CallMethodDef CallEntries[] = { {"_phylloptim_Leaf__find_psi_stem_from_psi_root", (DL_FUNC) &_phylloptim_Leaf__find_psi_stem_from_psi_root, 3}, {"_phylloptim_Leaf__evaluate_root_collar_psi", (DL_FUNC) &_phylloptim_Leaf__evaluate_root_collar_psi, 2}, {"_phylloptim_Leaf__dprofit_droot_collar_psi", (DL_FUNC) &_phylloptim_Leaf__dprofit_droot_collar_psi, 2}, + {"_phylloptim_Leaf__dprofit_droot_collar_psi_checked", (DL_FUNC) &_phylloptim_Leaf__dprofit_droot_collar_psi_checked, 2}, {"_phylloptim_Leaf__psi_stem_to_ci", (DL_FUNC) &_phylloptim_Leaf__psi_stem_to_ci, 3}, {"_phylloptim_Leaf__hydraulic_cost_Sperry", (DL_FUNC) &_phylloptim_Leaf__hydraulic_cost_Sperry, 3}, {"_phylloptim_Leaf__hydraulic_cost_TF", (DL_FUNC) &_phylloptim_Leaf__hydraulic_cost_TF, 2}, @@ -2441,7 +2456,7 @@ static const R_CallMethodDef CallEntries[] = { {"_phylloptim_gradient_output_names", (DL_FUNC) &_phylloptim_gradient_output_names, 0}, {"_phylloptim_gradient_batch_prepare", (DL_FUNC) &_phylloptim_gradient_batch_prepare, 9}, {"_phylloptim_gradient_batch_check", (DL_FUNC) &_phylloptim_gradient_batch_check, 1}, - {"_phylloptim_gradient_batch_run", (DL_FUNC) &_phylloptim_gradient_batch_run, 8}, + {"_phylloptim_gradient_batch_run", (DL_FUNC) &_phylloptim_gradient_batch_run, 10}, {"_phylloptim_root_network_from_carbon", (DL_FUNC) &_phylloptim_root_network_from_carbon, 4}, {NULL, NULL, 0} }; diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp index f83f997..63140d4 100644 --- a/src/RcppR6.cpp +++ b/src/RcppR6.cpp @@ -108,6 +108,10 @@ double Leaf__dprofit_droot_collar_psi(phylloptim::RcppR6::RcppR6dprofit_droot_collar_psi(opt_root_psi); } // [[Rcpp::export]] +std::vector Leaf__dprofit_droot_collar_psi_checked(phylloptim::RcppR6::RcppR6 obj_, double opt_root_psi) { + return obj_->dprofit_droot_collar_psi_checked(opt_root_psi); +} +// [[Rcpp::export]] double Leaf__psi_stem_to_ci(phylloptim::RcppR6::RcppR6 obj_, double psi_stem, double psi_upstream) { return obj_->psi_stem_to_ci(psi_stem, psi_upstream); } diff --git a/src/gradient.cpp b/src/gradient.cpp index a52b926..0884dc8 100644 --- a/src/gradient.cpp +++ b/src/gradient.cpp @@ -146,7 +146,8 @@ Rcpp::List gradient_batch_run(phylloptim::RcppR6::RcppR6 obj_, SEXP drivers, Rcpp::NumericMatrix theta, Rcpp::IntegerVector pars, double step, double stationarity_tol, std::string method, - bool fast_stem_curve) { + bool fast_stem_curve, SEXP psi, + SEXP dpsi_dtheta) { DriverBatch* batch_drivers = checked(drivers); const std::size_t n = batch_drivers->size(); const std::size_t npars = static_cast(pars.size()); @@ -186,11 +187,36 @@ Rcpp::List gradient_batch_run(phylloptim::RcppR6::RcppR6 obj_, // one thing that has to stay simple here. const std::vector par_idx(pars.begin(), pars.end()); + // `psi` is NULL to solve, or one collar potential per observation to impose; + // `dpsi_dtheta` is then an n x npars matrix. R validates both -- shape, + // positivity and the `method` conflict -- so this only holds the shape it was + // promised, and says so rather than reading past the end if it was not. + Rcpp::NumericVector psi_v; + Rcpp::NumericVector dpsi_v; + const double* psi_p = nullptr; + const double* dpsi_p = nullptr; + if (psi != R_NilValue) { + psi_v = Rcpp::as(psi); + if (static_cast(psi_v.size()) != n) { + Rcpp::stop("gradient_batch_run(): `psi` must have one value per " + "observation"); + } + psi_p = psi_v.begin(); + if (dpsi_dtheta != R_NilValue) { + dpsi_v = Rcpp::as(dpsi_dtheta); + if (static_cast(dpsi_v.size()) != n * npars) { + Rcpp::stop("gradient_batch_run(): `dpsi_dtheta` must be n x npars"); + } + dpsi_p = dpsi_v.begin(); + } + } + const std::vector res = phylloptim::gradient::batch(*obj_, theta.begin(), theta_nrow, *batch_drivers, obj_->supply_kind_name() == "single", - par_idx.data(), npars, settings); + par_idx.data(), npars, settings, psi_p, + dpsi_p); const int n_out = phylloptim::gradient::n_outputs; Rcpp::NumericMatrix value(static_cast(n), n_out); @@ -201,12 +227,21 @@ Rcpp::List gradient_batch_run(phylloptim::RcppR6::RcppR6 obj_, Rcpp::CharacterVector message(static_cast(n)); Rcpp::NumericVector H(static_cast(n)); Rcpp::NumericVector stationarity(static_cast(n)); + Rcpp::NumericMatrix M(static_cast(n), static_cast(npars)); + Rcpp::NumericMatrix dY_dpsi(static_cast(n), n_out); + Rcpp::NumericVector psi_used(static_cast(n)); for (std::size_t i = 0; i < n; ++i) { const phylloptim::gradient::Result& r = res[i]; for (int j = 0; j < n_out; ++j) { value(static_cast(i), j) = r.value[j]; + dY_dpsi(static_cast(i), j) = r.dY_dpsi[j]; + } + for (std::size_t k = 0; k < npars; ++k) { + M(static_cast(i), static_cast(k)) = + k < r.M.size() ? r.M[k] : NA_REAL; } + psi_used[i] = r.psi; // R's array layout: [i, k, j] lives at i + n*k + n*npars*j. for (std::size_t k = 0; k < npars; ++k) { for (int j = 0; j < n_out; ++j) { @@ -216,7 +251,17 @@ Rcpp::List gradient_batch_run(phylloptim::RcppR6::RcppR6 obj_, } } status[i] = phylloptim::gradient::status_name(r.status); - method_used[i] = r.used_ift ? "ift" : "fd"; + // `method` names the route, and a prescribed collar is a third one: neither + // "ift" (nothing was derived from -M/H) nor "fd" (the solve was not + // differenced). + // + // ⚠️ `psi_p` is CALL-level, so this does not vary across `i` -- the whole + // batch is prescribed or none of it is. A row that threw therefore reports + // `status = "error"` with `method = "prescribed"`, which is right: `method` + // says which route was asked for, and `status` says what happened to the row. + method_used[i] = psi_p != nullptr + ? "prescribed" + : (r.used_ift ? "ift" : "fd"); message[i] = r.message; H[i] = r.H; stationarity[i] = r.stationarity; @@ -229,5 +274,6 @@ Rcpp::List gradient_batch_run(phylloptim::RcppR6::RcppR6 obj_, Rcpp::_["value"] = value, Rcpp::_["gradient"] = grad, Rcpp::_["status"] = status, Rcpp::_["method"] = method_used, Rcpp::_["H"] = H, Rcpp::_["stationarity"] = stationarity, - Rcpp::_["message"] = message); + Rcpp::_["M"] = M, Rcpp::_["dY_dpsi"] = dY_dpsi, + Rcpp::_["psi"] = psi_used, Rcpp::_["message"] = message); } diff --git a/tests/testthat/test-gradient-batch.R b/tests/testthat/test-gradient-batch.R index c0c77f2..52a4f74 100644 --- a/tests/testthat/test-gradient-batch.R +++ b/tests/testthat/test-gradient-batch.R @@ -548,3 +548,98 @@ test_that("a batch that has lost its C++ pointer says so rather than crashing", # dereferencing a double as a pointer is not an error the session survives. expect_error(phylloptim:::gradient_batch_check(1.0), "not a prepared batch") }) + +# --- a collar potential the caller supplies (#88) ----------------------------- + +test_that("the batch's prescribed psi agrees with leaf_gradient(), bit-for-bit", { + # The batch is a transcription and has to stay one, so the assertion is the + # same one the solving path gets: equality with the R reference, not agreement + # to a tolerance. Both routes of the new path are covered -- the collar it + # solved for handed back (which must reproduce the solve), and a collar + # deliberately off it. + pars <- c("vcmax_25", "stem_b", "cost_scale_TF24") + soils <- c(1.0, 2.0, 3.0) + b <- leaf_batch(psi_soil = soils, PPFD = 900) + a <- leaf_gradient_batch(b, pars = pars) + expect_true(all(a$status == "interior")) + + g <- leaf_gradient_batch(b, pars = pars, psi = a$psi, + dpsi_dtheta = -a$M / a$H) + expect_true(all(g$status == "prescribed")) + expect_true(all(g$method == "prescribed")) + expect_identical(g$gradient, a$gradient) + expect_identical(g$M, a$M) + expect_identical(g$dY_dpsi, a$dY_dpsi) + + off <- a$psi + 0.2 + h <- leaf_gradient_batch(b, pars = pars, psi = off) + for (i in seq_along(soils)) { + r <- leaf_gradient(psi_soil = soils[[i]], PPFD = 900, pars = pars, + psi = off[[i]]) + expect_identical(as.vector(h$gradient[i, , ]), as.vector(r$gradient), + label = paste("row", i)) + expect_identical(as.vector(h$M[i, ]), unname(r$M), label = paste("M", i)) + expect_identical(as.vector(h$dY_dpsi[i, ]), unname(r$dY_dpsi), + label = paste("dY_dpsi", i)) + expect_identical(h$psi[[i]], r$psi) + } +}) + +test_that("a clamped row costs that row and not the batch", { + # The same per-row discipline `status == "error"` has, for a case that is not + # an error: a tracking model will hand over collar potentials that have drifted + # outside the feasible interval, and taking out the whole likelihood evaluation + # for one of them would take out the draw. + pars <- c("vcmax_25", "stem_b") + soils <- c(1.0, 2.0, 3.0) + b <- leaf_batch(psi_soil = soils, PPFD = 900) + a <- leaf_gradient_batch(b, pars = pars) + + g <- leaf_gradient_batch(b, pars = pars, psi = c(99, a$psi[2:3])) + expect_identical(g$status, c("clamped", "prescribed", "prescribed")) + expect_true(all(is.na(g$gradient[1, , ]))) + expect_true(all(is.na(g$M[1, ]))) + # The clamped row still says where it was pulled to, and the rows either side + # are BIT-IDENTICAL to being asked alone -- the leaf must not carry the clamped + # row's state into its neighbours. + expect_lt(g$psi[[1]], 99) + expect_true(all(is.finite(g$value[1, ]))) + for (i in 2:3) { + solo <- leaf_gradient_batch(leaf_batch(psi_soil = soils[[i]], PPFD = 900), + pars = pars, psi = a$psi[[i]]) + expect_identical(as.vector(g$gradient[i, , ]), + as.vector(solo$gradient[1, , ]), + label = paste("row", i)) + } +}) + +test_that("the batch validates psi and dpsi_dtheta by shape", { + pars <- c("vcmax_25", "stem_b") + b <- leaf_batch(psi_soil = c(1.0, 2.0, 3.0), PPFD = 900) + expect_error(leaf_gradient_batch(b, pars = pars, psi = 2.0), + "one collar potential per observation") + expect_error(leaf_gradient_batch(b, pars = pars, psi = rep(-1, 3)), + "finite and positive") + expect_error(leaf_gradient_batch(b, pars = pars, psi = rep(2, 3), + method = "fd"), + "cannot be given with `psi`") + expect_error(leaf_gradient_batch(b, pars = pars, dpsi_dtheta = c(1, 2)), + "needs `psi`") + expect_error(leaf_gradient_batch(b, pars = pars, psi = rep(2, 3), + dpsi_dtheta = matrix(0, 2, 2)), + "3 x 2") + + # A vector is one value per parameter, shared by every observation -- and it + # must equal the matrix that spells that out, or the recycling is wrong. Each + # row is given its OWN psi* so none of them clamps; a shared psi across three + # soil potentials would, and then the columns compared here would be NA. + psi <- leaf_gradient_batch(b, pars = pars)$psi + shared <- leaf_gradient_batch(b, pars = pars, psi = psi, + dpsi_dtheta = c(0.1, 0.2)) + spelt <- leaf_gradient_batch(b, pars = pars, psi = psi, + dpsi_dtheta = matrix(c(0.1, 0.1, 0.1, + 0.2, 0.2, 0.2), 3, 2)) + expect_identical(shared$gradient, spelt$gradient) + expect_identical(unname(shared$gradient[, , "collar"]), + matrix(c(0.1, 0.1, 0.1, 0.2, 0.2, 0.2), 3, 2)) +}) diff --git a/tests/testthat/test-gradient.R b/tests/testthat/test-gradient.R index a3e982b..b531703 100644 --- a/tests/testthat/test-gradient.R +++ b/tests/testthat/test-gradient.R @@ -734,3 +734,227 @@ test_that("the setter's positional trait call cannot drift in arity", { expect_length(formals(leaf_model()$set_traits), 14L) expect_identical(names(leaf_traits()), names(formals(leaf_model()$set_traits))) }) + +# --- a collar potential the caller supplies (#88) ----------------------------- + +test_that("a prescribed psi at psi* reproduces the solving path bit-for-bit", { + # ⚠️ THE LOAD-BEARING TEST OF THE WHOLE FEATURE, and the reason `psi` reuses + # `.gradient_ift()` rather than getting a composite of its own. Give it back + # the collar it solved for, and the response of that collar that it derived, + # and the two paths are the same arithmetic -- so this is `identical()` and not + # a tolerance. Anything that made the prescribed path a second implementation + # would show up here as a last-bit difference rather than as a design note. + d <- grid_drivers(2.0) + pars <- c("vcmax_25", "stem_b", "cost_scale_TF24") + a <- do.call(leaf_gradient, c(d, list(pars = pars))) + expect_identical(a$status, "interior") + + b <- do.call(leaf_gradient, + c(d, list(pars = pars, psi = a$psi, dpsi_dtheta = -a$M / a$H))) + expect_identical(b$status, "prescribed") + expect_identical(b$method, "prescribed") + expect_identical(b$gradient, a$gradient) + expect_identical(b$M, a$M) + expect_identical(b$dY_dpsi, a$dY_dpsi) + expect_identical(b$psi, a$psi) + + # `collar` IS psi, so its column must be exactly what was supplied -- the one + # column whose value the caller controls outright. + expect_identical(unname(b$gradient[, "collar"]), unname(-a$M / a$H)) + + # And the equivalence is not an accident of the default: dpsi_dtheta = 0 gives + # the direct terms, which differ from the solving path by the indirect term. + z <- do.call(leaf_gradient, c(d, list(pars = pars, psi = a$psi))) + expect_identical(unname(z$gradient[, "collar"]), rep(0, length(pars))) + expect_true(all(abs(z$gradient[, "A"] - a$gradient[, "A"]) > 1e-12)) +}) + +test_that("a prescribed psi off the optimum matches a difference at fixed collar", { + # The claim the feature makes: at a collar the caller chose, dY/dtheta is the + # direct partial. Checked against a central difference computed the long way, + # at a psi deliberately away from psi* so the indirect term the solving path + # would add is large and its absence is visible. + d <- grid_drivers(2.0) + a <- do.call(leaf_gradient, c(d, list(pars = "vcmax_25"))) + psi_off <- a$psi + 0.2 + g <- do.call(leaf_gradient, c(d, list(pars = "vcmax_25", psi = psi_off))) + expect_identical(g$status, "prescribed") + # Far from stationary, which is the regime this exists for. + expect_gt(g$stationarity, 1e-3) + + l <- leaf_model(leaf_traits(), leaf_control(), leaf_supply_multilayer()) + v <- leaf_traits()$vcmax_25 + h <- max(abs(v), 1) * 1e-6 + at <- function(x, field) { + set_traits(l, leaf_traits(vcmax_25 = x)) + do.call(set_drivers, c(list(l), d)) + l$evaluate_root_collar_psi(psi_off) + l[[field]] + } + fields <- c(A = "assim_colimited_", profit = "profit_") + for (column in names(fields)) { + ref <- (at(v + h, fields[[column]]) - at(v - h, fields[[column]])) / (2 * h) + expect_equal(g$gradient["vcmax_25", column], ref, tolerance = 1e-7, + label = column) + } + + # profit's dY/dpsi is the EXACT dprofit/dpsi here, not a difference of it -- + # the same rule as the envelope, reaching the opposite answer because the point + # is not stationary. Compare against the primitive directly. + set_traits(l, leaf_traits()) + do.call(set_drivers, c(list(l), d)) + expect_identical(g$dY_dpsi[["profit"]], l$dprofit_droot_collar_psi(psi_off)) + expect_gt(abs(g$dY_dpsi[["profit"]]), 1) +}) + +test_that("a clamped psi reports itself and returns no gradient", { + # ⚠️ NOT AN ERROR, AND NOT THE DIRECT TERM EITHER. The collar used is `psi` + # clamped into the feasible interval, so it moves with the BOUND rather than + # with the caller's dpsi_dtheta -- the active-set problem arriving through the + # clamp. The direct term alone would be plausible and wrong, so it is withheld; + # a tracking model reaches these points routinely, so it is reported and not + # thrown. + d <- grid_drivers(2.0) + pars <- c("vcmax_25", "stem_b") + + g <- do.call(leaf_gradient, c(d, list(pars = pars, psi = 99))) + expect_identical(g$status, "clamped") + expect_identical(g$method, "prescribed") + expect_true(all(is.na(g$gradient))) + expect_true(all(is.na(g$M))) + expect_true(all(is.na(g$dY_dpsi))) + # What the point IS still comes back: `psi` is the collar actually used, and + # `value` describes it. That is what lets a caller see where it was pulled to. + expect_lt(g$psi, 99) + expect_true(all(is.finite(g$value))) + + # And it fires one step inside the end too, because dY/dpsi cannot be centred + # there -- a one-sided difference over a shortened interval is the same failure + # the solving path's clamp detector refuses. + a <- do.call(leaf_gradient, c(d, list(pars = pars))) + expect_identical( + do.call(leaf_gradient, c(d, list(pars = pars, psi = g$psi)))$status, + "clamped") + # ...while a psi comfortably inside it does not. + expect_identical( + do.call(leaf_gradient, c(d, list(pars = pars, psi = a$psi)))$status, + "prescribed") +}) + +test_that("an integer psi is not mistaken for a clamped one", { + # ⚠️ THE CLAMP TEST IS `identical()` AGAINST WHAT THE CALLER PASSED, and + # `identical(3, 3L)` is FALSE. So an integer `psi` -- `psi = 3L`, or the + # entirely ordinary `for (p in 2:5)`, since `2:5` yields integers -- was + # reported CLAMPED at a collar it had been given exactly, with an all-NA + # gradient and a `psi` in the result showing no movement, so there was nothing + # to notice. + # + # Neither sibling had it: `.gradient_check_psi_batch()` already coerced, and + # C++'s `util::identical` is `a == b`. R alone, and only because the checker's + # return was discarded. + d <- grid_drivers(2.0) + a <- do.call(leaf_gradient, c(d, list(pars = "vcmax_25"))) + whole <- 3 # inside the feasible interval at this point + dbl <- do.call(leaf_gradient, c(d, list(pars = "vcmax_25", psi = whole))) + int <- do.call(leaf_gradient, + c(d, list(pars = "vcmax_25", psi = as.integer(whole)))) + expect_identical(dbl$status, "prescribed") + expect_identical(int$status, "prescribed") + expect_identical(int$gradient, dbl$gradient) + expect_identical(int$psi, dbl$psi) + + # And the batch entry point agrees on the same input, which it did not before. + b <- leaf_batch(psi_soil = 2.0, PPFD = 900) + expect_identical( + leaf_gradient_batch(b, pars = "vcmax_25", + psi = as.integer(whole))$status, "prescribed") +}) + +test_that("an infeasible prescribed psi is no-gradient, not a sentinel zero", { + # ⚠️ `dprofit_droot_collar_psi` RETURNS A HARD 0.0 ON ITS SHUT-DOWN EXIT, and a + # bare zero is indistinguishable from a stationary point. The solving path got + # away with reading the value alone because `H` collapses to zero with it and + # `usable` catches the pair; the prescribed path never divides by `H`, so it + # would have adopted the sentinel AS dprofit/dpsi -- and with a non-zero + # `dpsi_dtheta`, which is the case this feature exists for, silently lost + # profit's indirect term at exactly the dry points a tracking model lives in. + # + # The fix reads the feasibility flag the header has always had and R could not + # reach. Most such points are caught as `clamped` first -- the shut-down state + # seats a collar of its own choosing -- so the test hands back exactly that + # collar, which is the case `clamped` cannot see. + d <- grid_drivers(6.0) # drier than psi_crit: shut down + solved <- do.call(leaf_gradient, c(d, list(pars = "vcmax_25"))) + expect_identical(solved$status, "no-gradient") + + g <- do.call(leaf_gradient, + c(d, list(pars = "vcmax_25", psi = solved$psi, + dpsi_dtheta = 1))) + # Unclamped -- it is the collar the shut-down state itself seated -- and still + # refused, which is the whole point of the test. + expect_identical(g$psi, solved$psi) + expect_identical(g$status, "no-gradient") + expect_true(all(is.na(g$gradient))) + expect_true(all(is.na(g$dY_dpsi))) + + b <- leaf_batch(psi_soil = 6.0, PPFD = 900) + expect_identical(leaf_gradient_batch(b, pars = "vcmax_25", + psi = solved$psi)$status, "no-gradient") +}) + +test_that("psi and dpsi_dtheta are validated together", { + d <- grid_drivers(2.0) + gr <- function(...) do.call(leaf_gradient, c(d, list(pars = "vcmax_25", ...))) + expect_error(gr(psi = -1), "positive") + expect_error(gr(psi = 0), "positive") + expect_error(gr(psi = c(1, 2)), "single") + expect_error(gr(psi = NA_real_), "finite") + # `method` is about a solved optimum, so it cannot come with a collar the + # caller chose -- "fd" especially, which would difference the solve and answer + # a question about the optimum instead. + expect_error(gr(psi = 2, method = "fd"), "cannot be given with `psi`") + expect_error(gr(psi = 2, method = "ift"), "cannot be given with `psi`") + # dpsi_dtheta without psi is meaningless: there the collar's response is + # derived, not supplied. + expect_error(gr(dpsi_dtheta = 1), "needs `psi`") + expect_error(gr(psi = 2, dpsi_dtheta = c(1, 2)), "length 1 or one value") + expect_error(gr(psi = 2, dpsi_dtheta = c(nope = 1)), "names must be exactly") + + # ⚠️ AN INFINITE `dpsi_dtheta` IS REFUSED, AND THE ENVELOPE IS WHY. The + # composite is `direct + dY_dpsi * dpsi_dtheta`, so four columns come back + # +-Inf -- but `profit` is ASSIGNED from the direct term at a stationary psi + # rather than multiplied through, so it came back finite and plausible beside + # them. Measured at `psi = psi*` with `dpsi_dtheta = Inf`: A, gc, psi_stem and + # collar all Inf, profit 0.0105. Reading `profit` alone is plant's own case + # (#87), so the one surviving column was the one most likely to be believed. + # + # `anyNA` already caught NA and NaN; `Inf` it did not, and `is.finite` is the + # single check that covers all three. + expect_error(gr(psi = 2, dpsi_dtheta = Inf), "finite") + expect_error(gr(psi = 2, dpsi_dtheta = -Inf), "finite") + expect_error(gr(psi = 2, dpsi_dtheta = NaN), "finite") + # At psi* specifically, which is where the envelope fires and where the + # asymmetry was, rather than only at an arbitrary prescribed collar. + solved <- do.call(leaf_gradient, c(d, list(pars = "vcmax_25"))) + expect_error(gr(psi = solved$psi, dpsi_dtheta = Inf), "finite") + + # The batch entry point validates the same way, on BOTH shapes it accepts -- + # the matrix form never reaches `.gradient_dpsi_dtheta()`, so one check there + # would have left it open. + b <- leaf_batch(psi_soil = c(2.0, 2.5), PPFD = 900) + expect_error(leaf_gradient_batch(b, pars = "vcmax_25", psi = c(3, 3), + dpsi_dtheta = Inf), "finite") + expect_error(leaf_gradient_batch(b, pars = "vcmax_25", psi = c(3, 3), + dpsi_dtheta = matrix(c(1, Inf), nrow = 2)), + "finite") + + # Named is matched by name, not by position -- the same discipline `pars` has, + # for the same reason. + pars <- c("vcmax_25", "stem_b") + by_name <- do.call(leaf_gradient, + c(d, list(pars = pars, psi = 3.0, + dpsi_dtheta = c(stem_b = 2, vcmax_25 = 1)))) + by_pos <- do.call(leaf_gradient, + c(d, list(pars = pars, psi = 3.0, dpsi_dtheta = c(1, 2)))) + expect_identical(by_name$gradient, by_pos$gradient) +}) diff --git a/vignettes/cpp-interface.Rmd b/vignettes/cpp-interface.Rmd index 9f95151..7ec5813 100644 --- a/vignettes/cpp-interface.Rmd +++ b/vignettes/cpp-interface.Rmd @@ -130,6 +130,26 @@ const std::vector out = // out[i].status is Interior / Pinned / NoGradient / Error ``` +To evaluate at a collar potential you are carrying rather than one to be solved +for — the case a model that *tracks* the optimum has — pass `psi`, one per +observation, and optionally an `n × npars` column-major `dpsi_dtheta`: + +```cpp +const double psi[2] = {2.7, 2.8}; // positive magnitudes, MPa +const std::vector out = + g::batch(l, theta, 1, {d, d}, /*single=*/false, pars, 1, g::Settings(), + psi, /*dpsi_dtheta=*/nullptr); // null here means ZERO... +// out[i].status is now Prescribed or Clamped +// out[i].M / out[i].H are the coefficients of ds/dt = k(M + H s) +``` + +⚠️ **A null `dpsi_dtheta` means zero here and "derive it" one layer down**, which +is the one place this API can bite: `gradient_ift`'s null is the implicit function +theorem, correct for a collar that was *solved* for, and passing it through would +hand a prescribed caller the solving path's indirect term on top of their partial. +`at()` materialises the zeros so no consumer has to know that. For one +observation, use `g::at()` with a `g::Prescribed`. + `theta` is `g::n_pars` = 16 doubles: the fourteen traits in `set_traits`' argument order, then `leaf_specific_conductance_max`, then the single-potential path's `resistance`. `g::par_names()` is that order and every position has a name @@ -162,6 +182,7 @@ The ingredients are still public, if you want to compose something else: | `set_traits(...)` | replace all fourteen traits, refreshing splines and the temperature cache | | `perturb_stem_b(b)` | move `stem_b` by rescaling the curve instead of rebuilding it | | `dprofit_droot_collar_psi(psi, &feasible)` | dprofit/dψ at a given collar potential | +| `dprofit_droot_collar_psi_checked(psi)` | the same as `{dprofit, feasible}`, for R — RcppR6 has no form for a `bool*`, so the generated binding drops the flag | | `evaluate_root_collar_psi(psi)` | solve the inner problem at a fixed collar | ⚠️ **Pass the `feasible` flag.** `dprofit_droot_collar_psi` returns `0.0` on its diff --git a/vignettes/phylloptim.Rmd b/vignettes/phylloptim.Rmd index b1a9946..0abd9da 100644 --- a/vignettes/phylloptim.Rmd +++ b/vignettes/phylloptim.Rmd @@ -386,6 +386,38 @@ theorem reaches — at an interior optimum its gradient is the direct partial alone — which `?leaf_gradient` explains along with the measurement that makes that 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, +pass it and the gradient is taken there: + +```{r} +g <- leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "vcmax_25", psi = 2.7) +c(status = g$status, method = g$method) +``` + +`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. `stationarity` +still comes back and now says how far the point you chose is from the optimum; +`M`, `H` and `dY_dpsi` come back 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 exactly: + +```{r} +a <- leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "vcmax_25") +b <- leaf_gradient(psi_soil = 2.0, PPFD = 900, pars = "vcmax_25", + psi = a$psi, dpsi_dtheta = -a$M / a$H) +identical(a$gradient, b$gradient) +``` + +⚠️ **A `psi` outside the feasible collar interval gets no gradient.** It is +clamped, so the collar actually used moves with the *bound* rather than with your +`dpsi_dtheta`, and the direct term alone would be plausible and wrong. `status` +says `"clamped"`, the gradient is `NA`, and `g$psi` tells you where it was pulled +to. + Naming something undifferentiable is an error rather than a silent omission, and it tells you what is available: