Take a gradient at a collar psi you supply

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 = <the solved 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

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

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)
       },

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