--- title: "Differentiating the operating point" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Differentiating the operating point} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") ``` A model that contains this leaf and is itself differentiated, such as plant's stand gradient or a C++ fit of the model to observations, needs the derivative of the leaf's operating point with respect to its traits and its soil. This vignette is how the leaf supplies it: the three calls a consumer makes, what each operating point does, why the derivative comes from the implicit function theorem rather than from recording the solve, and what a consumer must not get wrong. It contains no runnable R; the C++ below is compiled and run by `tests/cpp/test_leaf.cpp`. ## One surface, two consumers The package differentiates the operating point one way, and two kinds of caller read it. | | `leaf_gradient()` and `gradient::batch` | a C++ model reading the surface | |---|---|---| | caller | R, or C++ calibrating a leaf against data | a C++ model that differentiates a run containing the leaf, such as plant's stand gradient | | scalar | a forward scalar with one lane per parameter and one for the decision variable; derivatives are returned as numbers | the caller's: `double`, a tangent, or an adjoint on the caller's tape | | method | implicit function theorem, with every row from automatic differentiation of closed forms | the same | | curves and routes | all nine, on both routes and both supplies | the same | | outputs | A, gc, psi_leaf, collar, profit | profit, carbon, each layer's water uptake and the residual (cuticular) flux; with `physical_outputs_at`, also A, E, psi_leaf, ci and gc | | classification | a status, derived from the kind | a kind (`interior`, `boundary-soil`, ...) written by the branch the solve took | `gradient::at` is a thin reader of the surface: it seats the parameters, solves, makes one pass through the three calls below, and reads the requested lanes. Its status is the kind, coarsened: | kind (C++) | status (R) | |---|---| | `interior` | `interior`, or `no-gradient` where the curvature is not below `-interior_curvature_floor` | | `boundary-soil`, `boundary-dry`, `hydraulic-shutdown` | `pinned` | | `solver-refused`, `non-finite-gradient`, `unsolved` | `no-gradient` | | `prescribed` | not reached on the solving path | ProfitMax with no objective (no light, or no stream to normalise over) is `no-gradient` whatever its kind. `prescribed`, `refused` and `no-gradient` are the statuses of the path that takes `psi`: `refused` means the collar fell outside the feasible interval, and `no-gradient` that the leaf potential it implies lies in the stem curve's tail, where every row would be rounding. The rest of this vignette is about the surface itself. ## Terms - **row**: one output's derivatives against a list of inputs. - **supplied row**: a row computed here and handed to the caller as data, with nothing recorded behind it. The caller's tape gets one statement per output, however long the leaf's arithmetic is. - **kind**: `Leaf::OperatingPointKind`, which condition placed the operating point. There are eight. - **draw**: `SupplyDraw`, the soil's water draw at one collar potential: the total uptake `E_up` and its slope in the collar, and each layer's uptake. It records the collar it was taken at. - **the three scalars**: the leaf's derivative surface is a set of member templates on a scalar `S`. At `double` they are an ordinary evaluation; at `odelia::ode::tangent_scalar` they carry a directional derivative; at `odelia::ode::adjoint_scalar` they record onto the caller's reverse tape. One source, not three copies. - **the pack**: `leaf_pars`, the differentiable parameters as one array of `n_leaf_pars` indexed by `par_*`. All but the last are what a calibration fits (`gradient::n_theta`); the last, `par_PPFD`, is the light the leaf stands in, seated per observation. ## The three calls After the ordinary solve at `double`, a consumer at scalar `S` makes three calls. `K` is the cost curve the solve used. ```cpp l.find_root_collar_psi(); // the solve, at double phylloptim::leaf_pars pars = /* the seated traits, at S, with directions set as needed */; phylloptim::SupplyInputs supply{psi_soil, r_R_H_min, r_R_V_sum, pars[phylloptim::par_root_P50], pars[phylloptim::par_root_c]}; auto draw = l.supply_draw_at(S(l.opt_root_psi_), supply); phylloptim::Leaf::HeldCoords held; S psi = l.operating_collar(draw, pars, curvature, &held); auto out = l.billed_outputs_at(psi, draw, pars, &held); // out.profit, out.carbon, out.uptake[j] and out.residual_flux carry their // derivatives in pars and supply; physical_outputs_at returns the same and // out.A, out.E, out.psi_leaf, out.ci and out.gc ``` `curvature` is `l.profit_curvature_at_collar()` at an interior point and is not read at any other kind. `profit_row_at_current_point` in `tests/cpp/test_leaf.cpp` is this sequence in full, compiled. `profit` is the objective that placed the point. `carbon` is what a plant grows on: `A` less the part of the curve's cost that is carbon spent, which each curve's `Curve` states (decision D6). Only TF24's hydraulic cost is carbon, so `carbon` is `profit` on TF24 and `A` on every curve whose cost is a price or a risk term. `physical_outputs_at` adds the point's own outputs at the live collar, for a model fitted to what a leaf is observed to do; it costs more, because on a curve whose price reads only the flux (TF24, TF24_floor, CF77) `billed_outputs_at` bills `carbon` from `profit` and the price at the soil's flux, and never forms the collar's channel at an interior point. `held` is optional. At an interior point, `operating_collar` solves the leaf potential and the intercellular CO2 at the collar on its way to the marginal, and `held` carries them to `billed_outputs_at`, which otherwise solves them again. Passing it moves no number: the forward tangent is bit-identical and the reverse rows agree to rounding. `held` lives on the caller's tape, like the draw, so it must not outlive the recording it was made on. At double and at every other kind `operating_collar` leaves it unplaced and `billed_outputs_at` solves its own. The calls assume, and mostly do not check: - the draw is taken at the solved collar (`operating_collar` refuses another); - `held`, if passed, comes from `operating_collar` at this point (`billed_outputs_at` refuses one placed at another collar or kind); - nothing moves the leaf between the solve and `billed_outputs_at`; - `pars` and `supply` hold, in value, the leaf's seated traits (`passive_pars()`) and its seated soil state and root network; - `supply` has the shape `Leaf::held_supply()` gives on the leaf's supply path: each layer's potential and resistances on the multi-layer path, or the one soil potential and the series resistance (`r_R_V_sum[0]`, the pack's `par_resistance`) on the single-potential path; - on ProfitMax, whose normalisers draw the supply again at its reference collar, the draw keeps its own copy of `supply` for that; no other curve's does; - with the thermal-damage layer on, its inputs are the pack's slots `par_A` to `par_c_repair_flux`: the capacities on the pack apply its shifted optimum, its functional fraction and its maintenance costs, as the solve's members do. With the energy balance on, the surface follows it on every curve: leaf temperature moves with the stem's flux and with the light that heats the leaf, and the capacities, the deficit, the cuticle and the first-order condition's temperature term carry its rows. ProfitMax's |A|max, which can then peak inside the stream rather than at `psi_crit`, takes its rows at that peak. ## The solve path, and the eight kinds The solve runs at `double` and is never differentiated. It brackets the collar between the wet bound, where uptake from the soil equals the cuticle's residual flux (zero at `g_min = 0`), and the dry bound, where the stem reaches its critical potential or the root reaches its own; then it root-finds the first-order condition `dprofit/dcollar = 0` inside the bracket, or reports that a bound binds. ``` prepare_collar_solve() kind := unsolved | +-- no interval to optimise over --------> hydraulic-shutdown +-- interval closed or inverted ---------> hydraulic-shutdown v root-find on dprofit/dcollar = 0 plain double, off any tape | +--> interior profit stationary in the collar +--> boundary-soil the wet bound binds +--> boundary-dry the stem's dry end binds +--> solver-refused the stationary point is a minimum +--> non-finite-gradient no usable gradient at an endpoint v supply_draw_at -> operating_collar -> billed_outputs_at ``` What `operating_collar` does at each kind: | kind | the collar is placed by | rows | |---|---|---| | `interior` | stationarity of profit in the collar | yes | | `boundary-soil` | the wet bound: uptake equals the residual flux | yes | | `boundary-dry` | the stem reaching its dry end: `psi_crit` for the curves bounded there (CF77, SOX, JW26, LeastCost), the end of the stem curve for the rest, where its potential is held (#183) | yes | | `hydraulic-shutdown` | the same wet bound, with the stem at the collar | yes, as `boundary-soil` | | `unsolved`, `prescribed`, `solver-refused`, `non-finite-gradient` | refused | none; `operating_collar` throws | A dry pin is decided at the end itself. The solve probes the marginal a step inside the dry end and, where it rises there, solves the end to collar precision and reads the marginal there: falling, the maximum is interior however close to the end it sits, and its curvature is then taken from a difference in the logarithm of its distance to the end, where the marginal is close to linear; rising, the leaf is pinned. The two answers are one continuous function of the drivers, which a consumer differencing a whole run depends on: a pin decided a step inside the end sat past its maximum, and a run crossing that decision jumped there. A narrow interval is polished, and a closed one is shut. When the feasible interval is narrower than `bracket_tol_at(bound_b)`, both of its ends are solved to collar precision before anything is decided on it; one no wider than `collar_root_tol` is `hydraulic-shutdown`, as is a soil from which the stem can carry nothing (past the end of its curve). Every open leaf therefore takes its stem from continuity. The collar route loses precision as the root resistance goes to zero (#176). **The leaf route** is the same surface with a different decision variable. There the stem's upstream end is held at the soil, so the leaf potential is what the solve places, and `decision_psi()` returns it; `supply_draw_at` carries the soil potential as an input of its own and walks no supply. The solve records the same kinds, read off which candidate the one solver returned: an optimum at the wet end, where the stem carries the residual flux from the soil, is `boundary-soil`; one at the stem's dry end is `boundary-dry`; a root of `dprofit/dpsi_leaf = 0` is `interior`; the one feasible potential, and ProfitMax with no objective, are `hydraulic-shutdown`; an argmax resolved only to a width is `solver-refused`. `operating_collar` places each by its condition in the leaf potential, and the coordinates follow the leaf potential one for one. On the stem curve's far tail the leaf route's profit is flat to ~1e-17 per MPa, so its argmax is undetermined there; `operating_collar` refuses such a point as flatter than `interior_curvature_floor`. A shut leaf is the ordinary evaluator at its wet bound (#182): `E` is the cuticle's residual flux, the conductance to CO2 is the cuticle's `rho g_min`, and `ci` solves the same supply-demand balance as an open leaf, above `ca` where the leaf is a net source. Nothing is written by hand, so the rows differentiate the value the solve reports and the shut leaf is continuous with the open one next to it. At `g_min = 0` there is no conductance and the shut leaf takes the `g_min -> 0` limit: the root of `A(ci) = 0`, or `A` as `ci -> infinity` where there is none. Light and heat shut no leaf: a leaf in the dark or too hot to fix carbon is an ordinary optimum, and the kind `shade-death` is gone. The shut stem sits at the collar rather than being held at a critical potential, so it keeps following a drying soil past `psi_crit` with no kink in the stem or root traits. ## Why the implicit function theorem Putting a trait on the tape and running the root-find live gives a derivative back through every iteration: the sensitivity of where that particular sequence of iterations stopped. It depends on the starting guess, the tolerance and the iteration count, converges to the right answer only as the tolerance tightens, and its error belongs to the solver's settings rather than to the model. The implicit function theorem gives the model's derivative directly. With the point `p*` defined by a residual `F(p*, theta) = 0`, ``` dp*/dtheta = -(dF/dp)^-1 dF/dtheta at p* ``` Both partials are of the residual, which is a closed form, so neither needs anything the solver did. The kind decides which residual: stationarity of profit at an interior point, zero flux at the wet bound, stem continuity at the stem's critical potential. One statement covers all of them. For profit at the optimum, ``` dPi*/du = dPi/du |_p - sum_j mu_j dc_j/du |_p ``` At an interior optimum every multiplier `mu_j` is zero and this is the envelope theorem: differentiate the objective holding the optimum fixed. At a bound the multiplier term is exactly the correction the envelope theorem leaves out, and `billed_outputs_at` supplies it by letting profit see the collar move through the bound's own implicit function. Only the objective is stationary in the collar, so only `profit` is taken with the collar held: `carbon` and the physical outputs carry the collar's movement at every kind. ## Why the kind is recorded, not inferred `dprofit_at_collar_psi` returns a literal `0.0` on its shutdown and reversed-gradient exits. A test of the obvious form, marginal profit near zero means an interior optimum, reads a shut leaf as stationary, and a curvature formed by differencing the same function agrees with it, so the two checks confirm each other and are both wrong. The kind is therefore written by the branch the solve exited and never inferred afterwards. `prepare_collar_solve` resets it to `unsolved` at the start of every solve, so a branch that forgets reports `unsolved` rather than the previous solve's kind. A consumer branches on `operating_point_kind()` and never on a returned value. A consumer that records a solve and replays it later puts the point back with `replay_operating_point(collar, kind)`, which hands over exactly the rows the solve did: a pinned solve places its collar at a polished bound and corrects its rows to that bound's root, and the replay re-takes both. ## Inside the surface **The draw is one walk.** The soil layers are in parallel, coupled only through the collar potential, so the draw's slope in the collar is a sum of per-layer terms and falls out of the same loop that accumulates the draw. The collar is stored passive in the draw, because the draw is the supply at a point; where the point moves is the supplied row's business, and taking the collar live would record the search that placed it. **Second derivatives without a nested tape.** `kernel_slope_at` returns a kernel's slope as its value and the kernel's mixed second partials as its rows. It evaluates the kernel once at a tangent-over-tangent scalar, both levels at `double`, so nothing reaches a tape; `supply_derivatives` then puts the result on the consumer's tape as one statement. An adjoint nested inside a tangent would cost many times the kernel's own statements, and odelia refuses that nesting at compile time. **Value from the table, rows from the curve.** Both vulnerability curves read their value from a pre-integrated table, because the solve ran on the table, and their trait rows from the closed form, because a table carries none. `closed_form_rows.hpp` states which query slope each uses. The two differ by about 1e-4 relative at the default 400 knots, which bounds how far a row can sit from a finite difference of the table. ## Extending it **`for_each_active` is a silent-zero hazard.** odelia walks an aggregate's active members through its `for_each_active`, and does nothing for an aggregate that has none. `PhotoCapacity`, `SupplyDraw`, `HeldCoords` and `Normalisers` declare it. A member added to either without being visited contributes no rows: its columns read exactly zero, every number stays finite, and nothing is raised. Visit every scalar member. **Another curve** is its `Curve`'s scalar forms (the `SurfaceForms` concept: its factor, its cost slope, the carbon part of its cost and, where they apply, its upstream slope and its price), beside the double forms the solve calls and held to them bit for bit by a test. Nothing in the surface is written per curve; a curve without its forms is refused at compile time. A curve that seeds normalisers before its search, as ProfitMax does, also needs them on the pack (`normalisers_at`). **A new fitted parameter** goes before `par_PPFD`, because `gradient::n_theta` is `par_PPFD`; a slot appended last falls outside every caller's `theta`. ## What is tested - `test_supplied_rows`: each layer's supplied collar slope against a central difference of the uptake. - `test_transpose`: a forward tangent against the reverse sweep through the same rows, on every curve at every kind it reaches, for every output of both calls. It cannot see a row that is wrong both ways. - `test_leaf`: each curve's scalar forms against its double forms, bit for bit; finite rows at every answered kind on a grid; a replay hands over its solve's rows bit for bit; TF24_floor at a zero price is TF24 on the surface, bit for bit, in every value and row; every output of both calls equal to the solve's own numbers at `double`; profit continuous, and every open leaf transpiring what its roots supply, across hydraulic shutdown and the light-compensation point; the collar, profit and output rows of every curve with the energy balance on, against a difference of the solve, at every kind each reaches; a second-difference sweep over soil potential, light, deficit and air temperature, with the energy balance off and on, one and three layers and the thermal-damage layer, finding no slope jump in profit, `A`, `E`, `psi_leaf`, `lambda`, or (with the energy balance off) two derivative rows at `g_min > 0`; the water budget `E = E_s + r` and `E_up = E`; the residual conductance's rows against a difference of the solve; the collar, profit, carbon, `A`, `E`, `psi_leaf`, `ci`, `gc`, uptake and residual-flux rows of every curve at every kind it reaches against a difference of the solve, and a shut leaf whose stem is seated at its collar billing its root traits nothing; the root curve's rows far along it against the closed form. Not yet tested: `boundary-dry` driven through `operating_collar` on a tape beyond `test_transpose`'s one fixture. ## Deliberately not built A fully analytic block, supplying every row with nothing recorded anywhere, is possible in principle. It needs eight closed forms that do not exist: four mixed second partials of the assimilation kernel and four of the cost. The supply side stays recorded either way, so the statements walked per leaf would fall by about half, not to the handful the kernels alone cost, and the transpose test that referees the current rows could not referee those, because a wrong closed form makes both routes wrong in the same way.