Flexible allometry: leaf and sapwood as ODE states

NEWS.md

@@ -7,6 +7,45 @@ entry gives the `old -> new` migration; the `plant-update-interface` skill
 (`.claude/skills/plant-update-interface/`) reads this section to migrate
 products using plant.
 
+* **The density coordinate is now chosen per model, and TF24 defaults to birth
+  date (#516).** `Control`'s boolean flag becomes a three-valued string, so that
+  a model can carry a default a passed `Control` does not silently override.
+  Migration:
+
+  * `control(node_density_in_birth_date = TRUE)`  -> `control(node_density_coordinate = "birth_date")`
+  * `control(node_density_in_birth_date = FALSE)` -> `control(node_density_coordinate = "height")`
+  * `ctrl$node_density_in_birth_date`             -> `ctrl$node_density_coordinate` (`"auto"`/`"birth_date"`/`"height"`)
+  * *(resolved value, per species/node)*          -> `species$density_in_birth_date` / `node$density_in_birth_date`
+
+  ⚠️ **This changes TF24 and TF24f results for identical inputs**, by several
+  times on offspring production, and bumps `TF24@v9 -> v10` (`TF24f@v9.1 ->
+  v10.1`). `"auto"` is the new default and resolves to `"birth_date"` for
+  TF24/TF24f and `"height"` for FF16/K93, so **FF16 and K93 are bit-identical**
+  — `auto` resolves to the same `false` the old default carried, and their
+  reference baselines are untouched.
+
+  Why it is not merely a numerical preference: the height coordinate's density
+  rate carries a compression term that equals `d(growth)/d(height)` only when
+  growth is a function of size alone. TF24's reserve pool feeds back into
+  growth, so the finite-difference probe moves height at fixed *absolute*
+  carbon and shifts the reserve fraction, whereas a cohort actually grows at
+  roughly constant reserve fraction — a different derivative, not a coarser one
+  (#590). Refinement confirms it: the birth-date answers are converged at the
+  default schedule while the height ones keep climbing and the exclusion ratio
+  widens instead of closing.
+
+  An unrecognised value is refused rather than read as one of the two
+  coordinates, since a typo would otherwise return a plausible number. The
+  empty string is accepted as a synonym for `"auto"`, matching how
+  `shading_model` already spells "let the model decide".
+
+  ⚠️ **A `Species` built outside a `Patch` cannot use the birth-date coordinate.**
+  The no-argument `introduce_new_node()` leaves every node carrying the same
+  birth date, which spans zero width, so the competition integral collapses to
+  exactly `0` and `log_densities` is `NA`. A real run cannot reach this — `Patch`
+  checks the birth dates are distinct and errors — but bare-`Species` diagnostic
+  code should set `node_density_coordinate = "height"` explicitly.
+
 * **An unknown trait name is now an error (#636).** `generate_strategy()` (and
   `add_strategies()` / `add_mutant()`, which route through it) refuse a trait
   whose name is not a parameter of the model, naming the offenders and their
@@ -539,6 +578,15 @@ were not previously recorded here:
   the two. The depth profile itself is not exposed — a fixed-width aux slot
   cannot carry a per-node vector — so the second half of #625 (reporting each
   leaf output against crown depth) remains open.
+* **TF24 sapwood can track the Huber value that maximises growth (`a_sw`, #516).** `d(log A_s departure)/dt` gains `a_sw * R_s`, an integral controller on the marginal *growth* return of sapwood area. It converges onto `R_s = 0` rather than tracking with an offset, and because integrating is averaging, a small `a_sw` makes the stem follow the long-run mean of a signal that swings with the weather — no extra tracked state.
+
+  The objective is growth, not net production, and the distinction is the whole design: against production extra sapwood is nearly always worth building, because it raises `k_max` and nothing charges the plant for the leaf area it did not build instead. `R_s` subtracts that price, which moves the optimum from far out (production's) in to where height growth actually peaks. That optimum is **not a single number**: it rises with height (0.87x the pipe-model ratio at 5 m, 1.05x at 10 m, 1.21x at 16 m, 1.30x at 20 m, moist soil) and rises again as soil dries, so no single fixed `theta` reproduces it. Measured after #617; the levels are the hydraulics' to set.
+
+  ⚠️ `d(profit)/d(kmax)` costs a second full leaf solve, so `a_sw > 0` roughly doubles the leaf work per rate evaluation. The envelope-theorem shortcut — re-evaluate at the collar potential already found and the indirect term vanishes — **does not apply**: `opt_root_psi_` is a root-find (supply matching demand), not a maximisation, so the derivative of profit with respect to it is not zero. Holding it fixed under-reports the sensitivity by 15–19%. The default `a_sw = 0` skips the whole block.
+
+  **`a_sw = 0` is the default and the model is then bit-exact** — the controller block is skipped entirely, and the one expression that changed unconditionally (`theta` becomes `theta * exp(psi)` in the conductance) is `theta * 1.0` at rest.
+
+  ⚠️ **Refused with the deep-crown shading model.** The sensitivity is measured at a single radiation and deep-crown integrates profit over crown positions, so the difference is not its derivative. `prepare_strategy()` throws rather than skipping it, because skipping leaves the controller with only its negative cost term and the stem shrinks without bound while the run looks plausible the whole way down.
 
 * **The NSC storage pool is bounded by the shape of its own flow (`TF24@v9`,
   `TF24f@v9.1`).** `dS/dt` was `net_flux > 0 ? net_flux : floor_gate * net_flux`

R/RcppExports.R

@@ -1521,6 +1521,10 @@ Node___FF16__FF16_Env__growth_rate_at_birth__get <- function(obj_) {
     .Call('_plant_Node___FF16__FF16_Env__growth_rate_at_birth__get', PACKAGE = 'plant', obj_)
 }
 
+Node___FF16__FF16_Env__density_in_birth_date__get <- function(obj_) {
+    .Call('_plant_Node___FF16__FF16_Env__density_in_birth_date__get', PACKAGE = 'plant', obj_)
+}
+
 Node___FF16__FF16_Env__fecundity__get <- function(obj_) {
     .Call('_plant_Node___FF16__FF16_Env__fecundity__get', PACKAGE = 'plant', obj_)
 }
@@ -1585,6 +1589,10 @@ Node___TF24__TF24_Env__growth_rate_at_birth__get <- function(obj_) {
     .Call('_plant_Node___TF24__TF24_Env__growth_rate_at_birth__get', PACKAGE = 'plant', obj_)
 }
 
+Node___TF24__TF24_Env__density_in_birth_date__get <- function(obj_) {
+    .Call('_plant_Node___TF24__TF24_Env__density_in_birth_date__get', PACKAGE = 'plant', obj_)
+}
+
 Node___TF24__TF24_Env__fecundity__get <- function(obj_) {
     .Call('_plant_Node___TF24__TF24_Env__fecundity__get', PACKAGE = 'plant', obj_)
 }
@@ -1649,6 +1657,10 @@ Node___TF24f__TF24_Env__growth_rate_at_birth__get <- function(obj_) {
     .Call('_plant_Node___TF24f__TF24_Env__growth_rate_at_birth__get', PACKAGE = 'plant', obj_)
 }
 
+Node___TF24f__TF24_Env__density_in_birth_date__get <- function(obj_) {
+    .Call('_plant_Node___TF24f__TF24_Env__density_in_birth_date__get', PACKAGE = 'plant', obj_)
+}
+
 Node___TF24f__TF24_Env__fecundity__get <- function(obj_) {
     .Call('_plant_Node___TF24f__TF24_Env__fecundity__get', PACKAGE = 'plant', obj_)
 }
@@ -1713,6 +1725,10 @@ Node___K93__K93_Env__growth_rate_at_birth__get <- function(obj_) {
     .Call('_plant_Node___K93__K93_Env__growth_rate_at_birth__get', PACKAGE = 'plant', obj_)
 }
 
+Node___K93__K93_Env__density_in_birth_date__get <- function(obj_) {
+    .Call('_plant_Node___K93__K93_Env__density_in_birth_date__get', PACKAGE = 'plant', obj_)
+}
+
 Node___K93__K93_Env__fecundity__get <- function(obj_) {
     .Call('_plant_Node___K93__K93_Env__fecundity__get', PACKAGE = 'plant', obj_)
 }

R/RcppR6.R

@@ -1,6 +1,6 @@
 ## Generated by RcppR6: do not edit by hand
 ## Version: 0.2.4
-## Hash: 2d29268c82e8129b38fbce4d1c832ee0
+## Hash: ece3258b408ecfd871ba1d78ada61c97
 
 ##' @importFrom Rcpp evalCpp
 ##' @importFrom R6 R6Class
@@ -1974,6 +1974,13 @@ Node <- function(T, E) {
           stop("Node<FF16,FF16_Env>$growth_rate_at_birth is read-only")
         }
       },
+      density_in_birth_date = function(value) {
+        if (missing(value)) {
+          Node___FF16__FF16_Env__density_in_birth_date__get(self)
+        } else {
+          stop("Node<FF16,FF16_Env>$density_in_birth_date is read-only")
+        }
+      },
       fecundity = function(value) {
         if (missing(value)) {
           Node___FF16__FF16_Env__fecundity__get(self)
@@ -2072,6 +2079,13 @@ Node <- function(T, E) {
           stop("Node<TF24,TF24_Env>$growth_rate_at_birth is read-only")
         }
       },
+      density_in_birth_date = function(value) {
+        if (missing(value)) {
+          Node___TF24__TF24_Env__density_in_birth_date__get(self)
+        } else {
+          stop("Node<TF24,TF24_Env>$density_in_birth_date is read-only")
+        }
+      },
       fecundity = function(value) {
         if (missing(value)) {
           Node___TF24__TF24_Env__fecundity__get(self)
@@ -2170,6 +2184,13 @@ Node <- function(T, E) {
           stop("Node<TF24f,TF24_Env>$growth_rate_at_birth is read-only")
         }
       },
+      density_in_birth_date = function(value) {
+        if (missing(value)) {
+          Node___TF24f__TF24_Env__density_in_birth_date__get(self)
+        } else {
+          stop("Node<TF24f,TF24_Env>$density_in_birth_date is read-only")
+        }
+      },
       fecundity = function(value) {
         if (missing(value)) {
           Node___TF24f__TF24_Env__fecundity__get(self)
@@ -2268,6 +2289,13 @@ Node <- function(T, E) {
           stop("Node<K93,K93_Env>$growth_rate_at_birth is read-only")
         }
       },
+      density_in_birth_date = function(value) {
+        if (missing(value)) {
+          Node___K93__K93_Env__density_in_birth_date__get(self)
+        } else {
+          stop("Node<K93,K93_Env>$density_in_birth_date is read-only")
+        }
+      },
       fecundity = function(value) {
         if (missing(value)) {
           Node___K93__K93_Env__fecundity__get(self)

R/scenario_eval.R

@@ -109,7 +109,7 @@ scenario_to_config <- function(row, mapping) {
 }
 
 ##' @return \code{scenario_control} returns the \code{Control} the gateway runs
-##'   under: \code{control()} with \code{node_density_in_birth_date = TRUE}.
+##'   under: \code{control()} with \code{node_density_coordinate = "birth_date"}.
 ##'
 ##'   Every scenario here is TF24, and TF24 is the model the two density
 ##'   coordinates genuinely disagree on (#590). The transport equation's
@@ -129,7 +129,10 @@ scenario_to_config <- function(row, mapping) {
 ##' @rdname scenario_eval
 ##' @export
 scenario_control <- function() {
-  control(node_density_in_birth_date = TRUE)
+  ## Set explicitly rather than left on "auto" (which TF24 now resolves to the
+  ## same thing), so the coordinate recorded in a stored scorecard's metadata is
+  ## the resolved one and not a value a reader has to resolve themselves.
+  control(node_density_coordinate = "birth_date")
 }
 
 ##' @param config A config list from \code{scenario_to_config}.
@@ -411,7 +414,10 @@ run_scenarios <- function(scenarios = read_scenario_table(),
   ## verdict -- so a stored scorecard that does not say which produced it cannot
   ## honestly be compared against another, and the blessed baseline is exactly
   ## such a stored scorecard.
-  meta$node_density_in_birth_date <- ctrl$node_density_in_birth_date
+  ## The requested coordinate. scenario_control() sets it explicitly rather than
+  ## leaving it on "auto", so what is stored here is the coordinate the run
+  ## actually used and not one the reader has to resolve against the model.
+  meta$node_density_coordinate <- ctrl$node_density_coordinate
   meta$max_patch_lifetime <- max_patch_lifetime
   attr(scorecard, "metadata") <- meta
   attr(scorecard, "max_patch_lifetime") <- max_patch_lifetime

inst/RcppR6_classes.yml

@@ -375,7 +375,7 @@ Control:
     - node_gradient_direction: int
     - node_gradient_richardson: bool
     - node_gradient_richardson_depth: size_t
-    - node_density_in_birth_date: bool
+    - node_density_coordinate: std::string
     - ode_step_size_initial: double
     - ode_step_size_min: double
     - ode_step_size_max: double
@@ -628,6 +628,11 @@ Node:
     ## Exact and current for the boundary node; frozen at its own birth for an
     ## introduced one.
     growth_rate_at_birth: {type: double, access: member, readonly: true}
+    ## Which coordinate this node's density is carried in, with Control's
+    ## "auto" already resolved against the model's own requirement. Read this
+    ## rather than restating which models prefer which coordinate: the density
+    ## rate and the boundary condition both differ between them.
+    density_in_birth_date: {type: bool, access: member, readonly: true}
     fecundity: {type: double, access: member}
     ode_size: {type: size_t, access: member}
     ode_state: {type: "std::vector<double>", access: function, name_cpp: "odelia::ode::r_ode_state", name_cpp_set: "odelia::ode::r_set_ode_state"}
@@ -672,8 +677,10 @@ Species:
     ## |dh/dtau| per node, boundary node last; the conversion between the two.
     height_jacobian: {type: "std::vector<double>", access: member, readonly: true}
     node_times: {type: "std::vector<double>", access: member}
-    ## Which coordinate this species' size distribution is carried in. Set from
-    ## the patch's control, not the strategy's (Patch::add_strategies).
+    ## Which coordinate this species' size distribution is carried in: the
+    ## patch's control (not the strategy's -- Patch::add_strategies), with
+    ## "auto" resolved against the model's own requirement. This is the resolved
+    ## value, so it is what a run should be reported against.
     density_in_birth_date: {type: bool, access: member, readonly: true}
     patch_densities: {type: "std::vector<double>", access: member, name_cpp: r_patch_densities}
     pr_patch_survival_at_birth: {type: "std::vector<double>", access: member, name_cpp: r_pr_patch_survival_at_birth}
@@ -1343,6 +1350,12 @@ TF24_Pars:
     - a_st1: double
     - a_st2: double
     - a_st3: double
+    - a_pl0: double
+    - a_pl1: double
+    - a_pl2: double
+    - a_pl3: double
+    - a_pl4: double
+    - a_sw: double
     - k_I: double
     # ⚠️ The leaf traits carry phylloptim's names (#634). One name per quantity,
     # spelled the same on both sides of the boundary; see TF24_Pars in

inst/include/plant/RcppR6_post.hpp

@@ -430,7 +430,7 @@ template <> inline SEXP wrap(const plant::Control& x) {
   ret["node_gradient_direction"] = Rcpp::wrap(x.node_gradient_direction);
   ret["node_gradient_richardson"] = Rcpp::wrap(x.node_gradient_richardson);
   ret["node_gradient_richardson_depth"] = Rcpp::wrap(x.node_gradient_richardson_depth);
-  ret["node_density_in_birth_date"] = Rcpp::wrap(x.node_density_in_birth_date);
+  ret["node_density_coordinate"] = Rcpp::wrap(x.node_density_coordinate);
   ret["ode_step_size_initial"] = Rcpp::wrap(x.ode_step_size_initial);
   ret["ode_step_size_min"] = Rcpp::wrap(x.ode_step_size_min);
   ret["ode_step_size_max"] = Rcpp::wrap(x.ode_step_size_max);
@@ -479,8 +479,8 @@ template <> inline plant::Control as(SEXP x) {
   ret.node_gradient_richardson = Rcpp::as<bool >(xl["node_gradient_richardson"]);
   // ret.node_gradient_richardson_depth = Rcpp::as<decltype(retnode_gradient_richardson_depth) >(xl["node_gradient_richardson_depth"]);
   ret.node_gradient_richardson_depth = Rcpp::as<size_t >(xl["node_gradient_richardson_depth"]);
-  // ret.node_density_in_birth_date = Rcpp::as<decltype(retnode_density_in_birth_date) >(xl["node_density_in_birth_date"]);
-  ret.node_density_in_birth_date = Rcpp::as<bool >(xl["node_density_in_birth_date"]);
+  // ret.node_density_coordinate = Rcpp::as<decltype(retnode_density_coordinate) >(xl["node_density_coordinate"]);
+  ret.node_density_coordinate = Rcpp::as<std::string >(xl["node_density_coordinate"]);
   // ret.ode_step_size_initial = Rcpp::as<decltype(retode_step_size_initial) >(xl["ode_step_size_initial"]);
   ret.ode_step_size_initial = Rcpp::as<double >(xl["ode_step_size_initial"]);
   // ret.ode_step_size_min = Rcpp::as<decltype(retode_step_size_min) >(xl["ode_step_size_min"]);
@@ -1373,6 +1373,12 @@ template <> inline SEXP wrap(const plant::TF24_Pars& x) {
   ret["a_st1"] = Rcpp::wrap(x.a_st1);
   ret["a_st2"] = Rcpp::wrap(x.a_st2);
   ret["a_st3"] = Rcpp::wrap(x.a_st3);
+  ret["a_pl0"] = Rcpp::wrap(x.a_pl0);
+  ret["a_pl1"] = Rcpp::wrap(x.a_pl1);
+  ret["a_pl2"] = Rcpp::wrap(x.a_pl2);
+  ret["a_pl3"] = Rcpp::wrap(x.a_pl3);
+  ret["a_pl4"] = Rcpp::wrap(x.a_pl4);
+  ret["a_sw"] = Rcpp::wrap(x.a_sw);
   ret["k_I"] = Rcpp::wrap(x.k_I);
   ret["vcmax_25"] = Rcpp::wrap(x.vcmax_25);
   ret["stem_P50"] = Rcpp::wrap(x.stem_P50);
@@ -1483,6 +1489,18 @@ template <> inline plant::TF24_Pars as(SEXP x) {
   ret.a_st2 = Rcpp::as<double >(xl["a_st2"]);
   // ret.a_st3 = Rcpp::as<decltype(reta_st3) >(xl["a_st3"]);
   ret.a_st3 = Rcpp::as<double >(xl["a_st3"]);
+  // ret.a_pl0 = Rcpp::as<decltype(reta_pl0) >(xl["a_pl0"]);
+  ret.a_pl0 = Rcpp::as<double >(xl["a_pl0"]);
+  // ret.a_pl1 = Rcpp::as<decltype(reta_pl1) >(xl["a_pl1"]);
+  ret.a_pl1 = Rcpp::as<double >(xl["a_pl1"]);
+  // ret.a_pl2 = Rcpp::as<decltype(reta_pl2) >(xl["a_pl2"]);
+  ret.a_pl2 = Rcpp::as<double >(xl["a_pl2"]);
+  // ret.a_pl3 = Rcpp::as<decltype(reta_pl3) >(xl["a_pl3"]);
+  ret.a_pl3 = Rcpp::as<double >(xl["a_pl3"]);
+  // ret.a_pl4 = Rcpp::as<decltype(reta_pl4) >(xl["a_pl4"]);
+  ret.a_pl4 = Rcpp::as<double >(xl["a_pl4"]);
+  // ret.a_sw = Rcpp::as<decltype(reta_sw) >(xl["a_sw"]);
+  ret.a_sw = Rcpp::as<double >(xl["a_sw"]);
   // ret.k_I = Rcpp::as<decltype(retk_I) >(xl["k_I"]);
   ret.k_I = Rcpp::as<double >(xl["k_I"]);
   // ret.vcmax_25 = Rcpp::as<decltype(retvcmax_25) >(xl["vcmax_25"]);

inst/include/plant/control.h

@@ -60,10 +60,23 @@ struct Control {
   bool   node_gradient_richardson;
   size_t node_gradient_richardson_depth;
 
-  // Carry the size distribution as a density in birth date rather than in
-  // height. The two resource integrals then run over introduction times, and
-  // the density rate is mortality alone.
-  bool node_density_in_birth_date;
+  // Which coordinate the size distribution's density is carried in:
+  //
+  //   "auto"       let the model decide (the default)
+  //   "birth_date" density in birth date; the two resource integrals run over
+  //                introduction times and the density rate is mortality alone
+  //   "height"     density in height; the density rate additionally carries a
+  //                compression term
+  //
+  // The height coordinate's compression term is the total derivative of the
+  // growth rate along a cohort's own trajectory, which equals the partial
+  // d(growth)/d(height) only when growth is a function of size alone. A model
+  // carrying state that feeds back into growth breaks that, so the coordinate
+  // is not a free numerical choice for every model -- hence "auto", resolved
+  // per model against Strategy::density_in_birth_date_default. An explicit
+  // value is an override and always wins, so both coordinates stay available
+  // for the comparison tests. See notes/516-flexible-allometry.md.
+  std::string node_density_coordinate;
 
   double ode_step_size_initial;
   double ode_step_size_min;
@@ -95,6 +108,28 @@ struct Control {
   double ci_niter;
 };
 
+// Resolve Control::node_density_coordinate for a model whose own preference is
+// `model_prefers_birth_date` (Strategy::density_in_birth_date_default). Only
+// "auto" consults the model; an explicit setting wins. An unrecognised value is
+// refused rather than silently treated as one of the two coordinates, because
+// the coordinates do not agree for a model with carried state -- a typo would
+// otherwise read as a real result.
+inline bool density_in_birth_date(const Control& control,
+                                  bool model_prefers_birth_date) {
+  const std::string& coord = control.node_density_coordinate;
+  if (coord == "auto" || coord.empty()) {
+    return model_prefers_birth_date;
+  }
+  if (coord == "birth_date") {
+    return true;
+  }
+  if (coord == "height") {
+    return false;
+  }
+  util::stop("Unknown node_density_coordinate \"" + coord +
+             "\": expected \"auto\", \"birth_date\" or \"height\"");
+}
+
 inline odelia::ode::OdeControl make_ode_control(const Control& control) {
   return odelia::ode::OdeControl(control.ode_tol_abs,
                          control.ode_tol_rel,