Differentiate profit, the output plant bills

.claude/CLAUDE.md

@@ -88,10 +88,12 @@ tests/testthat/gradient_golden.tsv
                                on that platform only, like tests/cpp/golden/ --
                                these are derivatives of argmax-evaluated outputs,
                                so they inherit its sqrt-amplified class and
-                               disagree cross-platform by up to 1.3e-3 -- ten
-                               times the solved outputs, because a finite
-                               difference divides the solver floor by the step
-                               and `R_d_25`'s step is the smallest here
+                               disagree cross-platform by up to 2.34e-3 (was
+                               1.3e-3 before #87 added `profit`) -- twenty times
+                               the solved outputs, because a finite difference
+                               divides the solver floor by the step, `R_d_25`'s
+                               step is the smallest here, and `profit` is the
+                               smallest magnitude carrying that same floor
 tests/validate/                R scripts comparing against plant (needs R)
 CMakeLists.txt                 the no-R build: C++ and Python consumers, and the
                                thing that makes "does not need R" runnable

NAMESPACE

@@ -3,6 +3,7 @@
 S3method(print,leaf_batch)
 export(Leaf)
 export(RootNetwork)
+export(gradient_output_names)
 export(gradient_par_names)
 export(leaf_batch)
 export(leaf_control)

NEWS.md

@@ -162,6 +162,66 @@ 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.
 
+## The gradient differentiates `profit`, which is what a demographic caller bills
+
+`leaf_gradient()` and `leaf_gradient_batch()` return a fifth column. The four
+that were there — `A`, `gc`, `psi_stem`, `collar` — are what a gas-exchange
+calibration observes, and they were chosen for the customer this feature was
+built for. They are **disjoint** from what `plant` reads off a solved leaf: its
+carbon budget is `leaf.profit_` (not `assim_colimited_`) and its water budget is
+`leaf.soil_consumption_`. So no trait gradient this package produced reached a
+demographic model at all, at the optimum or anywhere else
+([#87](https://github.com/traitecoevo/phylloptim/issues/87)).
+
+⚠️ **`profit` is the one output the envelope theorem reaches, and the only place
+this package uses it.** At an interior optimum `dprofit/dpsi = 0`, so the
+indirect term `(dprofit/dpsi)(dpsi*/dtheta)` vanishes identically and
+`dprofit/dtheta` is the direct partial at fixed ψ — no `dY/dpsi`, no `−M/H`. It
+is set from that term rather than computed through the composite, exactly as
+`collar` is set from `dpsi*/dtheta` rather than differenced.
+
+**That is a numerical decision, not a tidiness one, and the measurement is the
+reason.** The dropped term is *noise*, not an `h²` truncation: `profit` is the
+maximum, so it is flat, and a central difference of it divides the solve's ~1e-09
+floor by a ~1e-06 step. Over the golden grid's 136 interior rows —
+
+| | median | max |
+|---|---|---|
+| `\|dprofit/dpsi\|`, exact (forward AD) | 4.8e-15 | 5.4e-10 |
+| `\|dprofit/dpsi\|`, central difference | 7.8e-10 | 2.1e-04 |
+| relative move in `dprofit/dtheta` if kept | 2.7e-10 | 8.0e-05 |
+
+— eleven orders between the two instruments at the median, and the worst row sits
+in the band this repo calls a real difference rather than rounding. The identity
+is applied only where `status == "interior"`; at a pinned optimum `psi*` is a
+trait-dependent bound, `dprofit/dpsi` is not zero, and `profit` takes the same
+finite-difference fallback as the other four.
+
+The four existing columns are **bit-identical** — this is additive, like
+appending to `gradient_par_names()`. `tests/testthat/gradient_golden.tsv` gains a
+column and no existing cell moved, checked against master rather than against the
+branch point.
+
+⚠️ **The shut-down row's profit column is asserted against a closed form, not
+only recorded.** It is the one regime where `profit_` is written by a branch that
+leaves the other outputs alone (hazard 8), so a hex with nothing saying what it
+ought to be would pin a number rather than a fact. There `E = 0`, so `A = -R_d`
+exactly and the hydraulic cost does not depend on `R_d_25`: `dprofit/dR_d_25` is
+**−1**. And the shut-down collar is pinned at `psi_crit`, so
+`dcollar/dpsi_crit` is **1** — which is why `psi_crit` alone carries a non-zero
+profit gradient there.
+
+**`gradient_output_names()` is exported**, and R now *reads* the list rather than
+keeping a second copy. `gradient_par_names()` has to be duplicated-and-compared
+because R builds `theta` before any C++ call; the outputs have no such
+constraint, so adding one is a single edit.
+
+⚠️ **`uptake` was considered and is not here.** Every output must be a field R
+*copies* out of `operating_point_values()`; `uptake` is one R *computes*, by
+summing over the finite soil layers, so adding it means reproducing that
+summation and its order on the C++ side too. That is a separate decision from
+this one.
+
 ## An out-of-domain transport lookup says which spline, and which caller
 
 The stem curve is the only interpolator here built with extrapolation disabled, so

R/RcppExports.R

@@ -774,6 +774,26 @@ gradient_par_names <- function() {
     .Call('_phylloptim_gradient_par_names', PACKAGE = 'phylloptim')
 }
 
+#' The differentiated outputs, in the order C++ indexes them
+#'
+#' `A`, `gc`, `psi_stem`, `collar` and `profit`. Unlike [gradient_par_names()],
+#' which R holds its own copy of and a test compares, R **reads** this one — so
+#' the list exists once and adding an output is one edit rather than two that
+#' can disagree.
+#'
+#' The first four are what a gas-exchange calibration observes; `profit` is what
+#' a demographic consumer bills, and it is the one output the envelope theorem
+#' reaches. See [leaf_gradient()].
+#'
+#' @return A character vector of five names.
+#' @seealso [leaf_gradient()], [leaf_gradient_batch()]
+#' @examples
+#' gradient_output_names()
+#' @export
+gradient_output_names <- function() {
+    .Call('_phylloptim_gradient_output_names', PACKAGE = 'phylloptim')
+}
+
 gradient_batch_prepare <- function(root_network, psi_soil, soil_depth, PPFD, atm_vpd, ca, leaf_temp, atm_o2_kpa, atm_kpa) {
     .Call('_phylloptim_gradient_batch_prepare', PACKAGE = 'phylloptim', root_network, psi_soil, soil_depth, PPFD, atm_vpd, ca, leaf_temp, atm_o2_kpa, atm_kpa)
 }

R/gradient-batch.R

@@ -188,7 +188,7 @@ print.leaf_batch <- function(x, ...) {
 ##' observation.
 ##'
 ##' @section What it computes, and where the maths is written down:
-##' The same four derivatives at the same solved operating point, by the same two
+##' 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
 ##' [leaf_gradient()]'s documentation for all of it** — the implicit function
 ##' theorem, why the second term is not a correction, why the premise is tested
@@ -239,9 +239,9 @@ print.leaf_batch <- function(x, ...) {
 ##'
 ##' @return A list with
 ##'   \describe{
-##'     \item{`gradient`}{an `n` × `length(pars)` × 4 array, the last dimension
-##'       being `A`, `gc`, `psi_stem` and `collar`}
-##'     \item{`value`}{an `n` × 4 matrix of the solved outputs}
+##'     \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{`H`, `stationarity`}{the curvature and the implied Newton step
@@ -314,8 +314,8 @@ leaf_gradient_batch <- function(batch,
                             match(pars, par_names) - 1L, step,
                             stationarity_tol, method, fast_stem_curve)
 
-  dimnames(res$gradient) <- list(NULL, pars, .gradient_output_names)
-  dimnames(res$value) <- list(NULL, .gradient_output_names)
+  dimnames(res$gradient) <- list(NULL, pars, .gradient_output_names())
+  dimnames(res$value) <- list(NULL, .gradient_output_names())
   res
 }
 

R/gradient.R

@@ -96,8 +96,9 @@ set_traits <- function(x, traits) {
 
 ##' Trait gradients of a solved operating point
 ##'
-##' The derivatives of the solved outputs with respect to the traits: `dA/dtheta`,
-##' `dgc/dtheta`, `dpsi_stem/dtheta` and `dcollar/dtheta`, at one operating point.
+##' The derivatives of the solved outputs with respect to the traits:
+##' `dA/dtheta`, `dgc/dtheta`, `dpsi_stem/dtheta`, `dcollar/dtheta` and
+##' `dprofit/dtheta`, at one operating point.
 ##'
 ##' @section The maths, and why it is not just a finite difference:
 ##' The reported outputs are evaluated at the profit-maximising collar potential
@@ -120,6 +121,44 @@ set_traits <- function(x, traits) {
 ##' The second term is not a correction. For `cost_scale_TF24`, `beta2`, `stem_b`
 ##' and `stem_c` it is 100% of the answer, and for `vcmax_25` 52%.
 ##'
+##' @section profit, which is the one output the envelope theorem reaches:
+##' The first four outputs are what a gas-exchange calibration observes. `profit`
+##' is here because it is what a **demographic** consumer bills: plant's carbon
+##' budget reads `leaf.profit_`, not `assim_colimited_`, so before it was added
+##' the two sets were disjoint and no gradient from this package reached a
+##' demographic model at all.
+##'
+##' It is also the cheapest of the five, because it is the objective. At an
+##' interior optimum `dprofit/dpsi = 0`, so the indirect term above vanishes
+##' identically and
+##'
+##' \deqn{dprofit/d\theta = \partial profit/\partial\theta|_{\psi^*}}
+##'
+##' -- the direct partial alone, with no `dY/dpsi` and no `-M/H`. That is the
+##' envelope theorem, and it is the only place this package uses it.
+##'
+##' ⚠️ **The term this drops is noise, not an `h^2` truncation, and that is why
+##' dropping it matters.** `profit` is the maximum, so it is flat, and a central
+##' difference of it at `psi*` divides the solve's ~1e-09 floor by a ~1e-06 step.
+##' Over the golden grid's 136 interior rows:
+##'
+##' | | median | max |
+##' |---|---|---|
+##' | `\|dprofit/dpsi\|`, exact (forward AD) | 4.8e-15 | 5.4e-10 |
+##' | `\|dprofit/dpsi\|`, central difference | 7.8e-10 | 2.1e-04 |
+##' | relative move in `dprofit/dtheta` if kept | 2.7e-10 | 8.0e-05 |
+##'
+##' Eleven orders between the two instruments at the median, and the worst row
+##' sits in the band this package calls a real difference. So `dprofit/dtheta` is
+##' set from the direct term rather than computed through the composite, exactly
+##' as `collar` is set from `dpsi*/dtheta` rather than differenced.
+##'
+##' ⚠️ **The envelope does not survive the active set.** At a pinned optimum
+##' `psi*` is a trait-dependent bound and `dprofit/dpsi` is not zero there, so the
+##' indirect term is real. This applies the identity only where `status` is
+##' `"interior"`; everywhere else `profit` goes through the same finite-difference
+##' fallback as the other four.
+##'
 ##' @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
@@ -288,7 +327,8 @@ set_traits <- function(x, traits) {
 ##' @return A list with
 ##'   \describe{
 ##'     \item{`gradient`}{a matrix, one row per trait in `pars`, with columns
-##'       `A` (umol C m^-2 s^-1 per trait unit), `gc`, `psi_stem` and `collar`}
+##'       `A` (umol C m^-2 s^-1 per trait unit), `gc`, `psi_stem`, `collar` and
+##'       `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}
@@ -497,12 +537,40 @@ leaf_gradient <- function(psi_soil,
       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.
+      #
+      # ⚠️ 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.
+      dY_dpsi[["profit"]] <- resid
     }
   }
 
+  # ⚠️ THE ENVELOPE THEOREM, and the only place this package uses it. At a
+  # STATIONARY point dprofit/dpsi is analytically zero, so profit's indirect term
+  # vanishes identically and dprofit/dtheta is the direct partial alone. The
+  # composite is told to ASSIGN that column rather than reach it by multiplying a
+  # 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 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)
+                  fast_stem_curve, envelope = envelope)
   } else {
     .gradient_fd(l, reset, theta, pars, step, fast_stem_curve)
   }
@@ -517,11 +585,34 @@ leaf_gradient <- function(psi_soil,
 
 # --- internals ---------------------------------------------------------------
 
-# The four differentiated outputs, in one place so the composite and the fallback
-# cannot disagree about what they are. `collar` is psi* itself, which makes
-# dcollar/dtheta equal to dpsi*/dtheta -- so the two routes compute the same
-# quantity by different means, and a test can compare them.
-.gradient_output_names <- c("A", "gc", "psi_stem", "collar")
+# The differentiated outputs, READ from C++ rather than restated here.
+#
+# ⚠️ This is deliberately NOT the pattern `.gradient_par_names` uses. That one is a
+# second copy kept honest by a test comparing it with `gradient_par_names()`, and
+# it has to be, because R builds `theta` from `leaf_traits()` and so needs the
+# order before any C++ call. The OUTPUT list has no such constraint, so it is a
+# single definition and there is no hazard comment to write: adding an output is
+# one edit in `gradient.hpp`, not two that can disagree.
+#
+# `collar` is psi* itself, which makes dcollar/dtheta equal to dpsi*/dtheta -- so
+# the two routes compute the same quantity by different means and a test can
+# compare them. `profit` is what plant CONSUMES (`leaf.profit_`, not
+# `assim_colimited_`) and is the one output the envelope theorem reaches; see
+# ?leaf_gradient.
+#
+# ⚠️ Cached at FIRST CALL, not at build time, for the reason
+# `.gradient_outputs_idx()` below records -- and here there is a second reason:
+# calling into the shared library while the package is still being sourced is not
+# something to rely on.
+.gradient_output_names <- local({
+  nms <- NULL
+  function() {
+    if (is.null(nms)) {
+      nms <<- gradient_output_names()
+    }
+    nms
+  }
+})
 
 # ⚠️ ONE call, not four field reads. Every `l$field` is an R6 ACTIVE BINDING -- a
 # closure call wrapping a `.Call` -- and this function runs once per perturbation, so
@@ -541,15 +632,15 @@ leaf_gradient <- function(psi_soil,
   idx <- NULL
   function() {
     if (is.null(idx)) {
-      idx <<- match(.gradient_output_names, .operating_point_names)
+      idx <<- match(.gradient_output_names(), .operating_point_names)
     }
     idx
   }
 })
 
 .gradient_outputs <- function(l) {
   v <- l$operating_point_values()[.gradient_outputs_idx()]
-  names(v) <- .gradient_output_names
+  names(v) <- .gradient_output_names()
   v
 }
 
@@ -815,7 +906,7 @@ 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) {
+                          fast_stem_curve = TRUE, envelope = FALSE) {
   seat <- .gradient_reseat_base(reset, theta, pars, fast_stem_curve)
   out <- t(vapply(seq_along(pars), function(k) {
     p <- pars[[k]]
@@ -846,8 +937,16 @@ leaf_gradient <- function(psi_soil,
     # difference of two identical numbers.
     g <- direct + dY_dpsi * dpsi_dtheta
     g[["collar"]] <- dpsi_dtheta
-    g[.gradient_output_names]
-  }, numeric(length(.gradient_output_names))))
+    # 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.
+    if (envelope) {
+      g[["profit"]] <- direct[["profit"]]
+    }
+    g[.gradient_output_names()]
+  }, numeric(length(.gradient_output_names()))))
   reset(theta)
   rownames(out) <- pars
   out
@@ -869,8 +968,8 @@ leaf_gradient <- function(psi_soil,
       l$find_root_collar_psi()
       .gradient_outputs(l)
     }
-    ((side(1) - side(-1)) / (2 * h))[.gradient_output_names]
-  }, numeric(length(.gradient_output_names))))
+    ((side(1) - side(-1)) / (2 * h))[.gradient_output_names()]
+  }, numeric(length(.gradient_output_names()))))
   rownames(out) <- pars
   reset(theta)
   out

README.md

@@ -324,10 +324,15 @@ traits, which is what a gradient-based optimiser or a Hamiltonian sampler wants:
 ```r
 g <- leaf_gradient(psi_soil = 2.0, PPFD = 900,
                    pars = c("vcmax_25", "stem_b", "cost_scale_TF24"))
-g$gradient   # rows: parameters.  columns: A, gc, psi_stem, collar
+g$gradient   # rows: parameters.  columns: A, gc, psi_stem, collar, profit
 g$method     # "ift" or "fd" -- see below
 ```
 
+The first four columns are what a gas-exchange calibration observes. `profit` is
+the objective, and it is there for a demographic consumer: `plant` bills carbon
+from the leaf's profit rather than its assimilation, so without that column no
+gradient from this package reached a demographic model at all.
+
 `pars` is not restricted to traits: `leaf_specific_conductance_max` and, on the
 single-potential path, `resistance` are differentiable too, because a calibration
 fits them and nothing in the derivation cares whether a parameter is a trait.
@@ -345,6 +350,11 @@ moving that optimum — and for `cost_scale_TF24`, `beta2`, `stem_b` and `stem_c
 the second route is **100%** of the answer. Differentiating the optimality
 condition rather than the solved output gets both terms exactly.
 
+`profit` is the exception, and it is the cheapest column for the reason it is the
+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.
+
 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