Prototype reverse-mode AD of SCM outputs

[dfalster]

Adds exact reverse-mode automatic differentiation of SCM outputs with respect
to traits and the birth-rate driver, for FF16, TF24 and TF24f. One backward
sweep returns the derivative of a scalar emergent output w.r.t. all inputs at
once — the shape of a calibration objective or a selection gradient — at
machine precision, with no per-trait re-run and no step-size lottery.

The C++ core gains a third template axis, the scalar type S, defaulting to
double. Every existing use and the whole R interface still instantiates
double and is bit-identical; the reference suite is untouched.

This is a prototype, not a finalised API: the access-point names and the engine
structure are expected to change. Every access point is either finite and
validated against a finite difference over the same reconstruction, or gated
with a clear error — never plausible-but-wrong. notes/ad-refactor-optimize-roadmap.md
has the follow-ups and the guide added here has the techniques and the limits.

Refs #472 (scope B), #537.

[dfalster]

@aornugent Could you provide a review on this before merging? We've demonstrated AD capability with one prototype implementation. We ended up implementing outside of the SCM object, but I'm not sure this is the right long-term solution. I'm thinking we land this PR and then re-implement, refine in a new round of work. I'm open to variants. I'm interested in your thoughts.

[aornugent]

This is huge!

The spike proves that this is possible — exact reverse-mode gradients through the SCM for the strategies I thought would be hardest, along with validation of that the results are accurate and flushing out the real difficulties on the way

I have opened #559 with a design roadmap for what comes next using this branch as the basis and keeping the tests you developed for regression.

Where #559 goes further is arguing that most of the surrounding machinery isn't needed. plant
already links odelia, which is where the ODE solver, the differentiable spline, and the XAD tape
live. Build on that and a lot falls away: the gradient becomes run_mutant/run differentiated
rather than three *_emergent.cpp engines; the R-side Rcpp::as<Environment> round-trip becomes
a native boundary; stand_stage_history collapses into one "record the adaptive nodes, replay
them fixed" step; <T,E,S> goes back to <T,E>; five hand-managed tapes become odelia's one.
The reasoning is all in the PR.

A useful side effect: because that machinery belongs in odelia rather than plant, it becomes
general AD infrastructure — reverse-mode gradients and Jacobians for any model built on odelia,
not just the SCM.

So the recommendation is to hold off merging the prototype and land the small odelia surface
first, using these Jacobians as the oracle, then port the plant path onto it. Merging as-is means
taking on a second tape lifecycle and scalar axis that the next step just unwinds. Would be good
to get your read on the direction before anyone starts building.

PS. I got here by working in a meta-repo with plant and odelia as submodules, which gave Claude access to both codebases and XAD at once (rather than working on plant with odelia installed).

[dfalster]

Thank you @aornugent ! I agree 100% with your proposed directions. I'll hold off on merging this as we develop a better long-term solution.

[dfalster]

Working detail, moved off the PR description so it does not land in git log. Under squash merge the description becomes the commit message verbatim; this comment is just as permanent and is linked from the squashed subject. See the family guideline.

Original description, unchanged:

---

Status: prototype

This is a prototype of automatic differentiation in plant — a working,
end-to-end-validated spike that establishes the capability and the techniques, not a
finalised API. It is deliberately additive and bit-identical to the reference model
(see below), so it is safe to land as a foundation, but the access-point names, the engine
structure, and several honest-scope gaps are expected to be refined later (see
"Deferred" and notes/ad-refactor-optimize-roadmap.md). Suggest we land in develop then refine in a new feature branch.

Summary

Adds exact reverse-mode automatic differentiation of plant's SCM outputs w.r.t.
traits and the birth-rate driver, for FF16, TF24 and TF24f (#472 scope B, #537). One
backward sweep returns the derivative of a scalar emergent output w.r.t. all inputs at
once — the "many inputs → one scalar" shape of a calibration objective or a selection
gradient — at machine precision, with no per-trait re-run and no step-size lottery.

The full narrative (how it was built, the techniques, the validation, the honest scope)
is in the guide added here: overstorey-staging/guides/autodiff-trait-gradients.qmd.
This description is the map; the guide is the territory.

The discipline: additive and bit-identical

The C++ core gains a third template axis — the scalar type S, defaulting to
double. Every differentiable piece is templated on S; AD only ever happens in C++ by
instantiating with an XAD active type (odelia's reverse-mode tape). Because S defaults
to double, every existing use and every line of the RcppR6/R interface still names the
double type and is bit-for-bit unchanged. The ordinary model, and the reference test
suite, are untouched — the AD is purely additive. The whole spike landed as a long series
of individually-verified, bit-identical steps rather than one rewrite.

Two recurring principles do all the work (both detailed in the guide):

1. Never differentiate a solver — differentiate its converged point. Root-finds via
   the implicit function theorem, optimisations via the envelope theorem, adaptive
   integrators via a frozen-schedule replay (a double discovery pass finds the
   schedule + resident light; a second pass replays it with the trait active).
2. The derivative of a built quantity must come from the same construction as the
   value — the same spline, the same harvested operating point.

Prototype access points

The prototype surface (names/signatures may change on refinement). All dispatch on
strategy; all are ordinary compiled package calls (no on-the-fly sourceCpp). Given a
run_scm(save_RK45_cache = TRUE) resident:

Call	Differentiates	Strategies
stand_gradient(scm, metrics, traits, species, feedback)	metrics × traits Jacobian of emergent stand outputs (canonical entry)	FF16, TF24f (all metrics); TF24 (offspring only)
offspring_production_gradient(scm, traits, species)	d(offspring_production)/d(trait) — the invasion / selection gradient	FF16, TF24, TF24f
stand_state_jacobian(scm, traits, species)	per-cohort state × trait Jacobian (escape hatch for any downstream metric)	FF16, TF24
grow_individual_to_size_gradient(individual, sizes, size_name, env, traits)	d(t*)/dθ and total d(state at t*)/dθ for one plant in a fixed env	FF16, TF24f
birth_rate_gradient(scm, metrics, species)	the birth-rate axis: resident d(census)/d(birth_rate) and d(net_reproduction_ratio)/d(birth_rate) = dR0/db (demographic equilibrium)	FF16

feedback = "frozen" (default) is the rare-mutant invasion gradient (canopy held fixed);
feedback = "resident" is the resident-total gradient via a coupled whole-stand replay
(every trait re-shades the canopy; on a multi-species stand, the cross-species total).

Validation

• Whole test suite green, unchanged reference behaviour: FAIL 0, PASS 2610 (make test / plain-R load_all); the AD-specific tests run against the installed DLL via
  make test-ad (FAIL 0, PASS 362).
• AD-vs-AD regression fixture (tests/testthat/fixtures/gradient-baseline.rds,
  scripts/gradient_fixture.R) pins every compiled gradient engine to its validated value
  — bit-identical for the relocation/frozen paths, noise-floor for the coupled paths.
• AD-vs-FD throughout: each engine is validated against a finite difference over the
  same reconstruction (the honest contract), plus the physics where available.
• FF16 emergent gradients are faithful to ~1e-13…1e-8; TF24/TF24f to the leaf-optimiser
  floor (~5e-7), as expected for a re-solved optimum.

Honest scope & known limits

Stated plainly in the guide's "Where it works, and where it breaks" section (swept
empirically across patch lifetime and species count):

• FF16 has no hard break in the tested envelope (patch lifetime → 200, → 5 species).
  Soft limits only: cost grows with cohort count (trait-count-free by design), and the
  frozen-grid reconstruction drifts gently (~1e-6 relative at long horizon).
• The coupled/resident gradients are the frozen-grid total — they freeze the ODE
  schedule + knot positions and re-evolve the stand. Against a full run_scm FD (which
  lets the adaptive grid re-respond) they differ ~16–30 % for geometry traits (~half is
  genuine grid response). This is the invasion/selection gradient's exact object; consumers
  wanting the fully-adaptive derivative should know the gap.
• TF24f resident census single-species gates at patch lifetime ≳ 5, and multi-species
  resident gates on a refined (clustered) schedule — both raise a clear error, never a
  NaN. Root cause is the frozen-step replay meeting a stiff log_density ↔ canopy
  feedback the live SCM only tames by adaptive stepping. Plausibly the same family as
  [TF24 hydraulics] SCM cohort-density blow-up under extreme seasonal drought #550 (the SCM size-density blow-up under drought) — same d(log_density)/dt = −dg/dh − mortality equation with a large dg/dh from the TF24 leaf — but a distinct
  manifestation: here run_scm stays healthy and the true gradient exists; only the replay
  diverges. [TF24 hydraulics] SCM cohort-density blow-up under extreme seasonal drought #550's root fixes (NSC [TF24 nsc] Enable resource storage, presumably NSCs, to buffer short-term variation in productivity #517 / hydraulic-optimum continuity) would help both.
• TF24 census/resident metrics are not available (only offspring_production); the
  census number density needs a leaf-optimiser cross-sensitivity the linearised harvest
  does not capture. TF24f (analytic tracked collar) gets census instead.

Every access point is finite-and-validated or gated-with-an-error — never
plausible-but-wrong, which is the property a downstream consumer needs.

Deferred (downstream / follow-ups)

Calibration, likelihoods, and adaptive-dynamics / equilibrium solvers live in downstream
packages (e.g. regnans) — plant's job is to expose the gradient cleanly. Follow-ups
noted in notes/ad-refactor-optimize-roadmap.md: TF24 census; the ResidentHarvest seam
for multi-patch/stochastic variants; standalone Leaf calibration; TF24f resident
birth-rate; cross-species net_reproduction_ratio; and (the one change that would retire
most gated cases) adaptive sub-stepping in the replay.

Notes for review

• Large but layered. ~17.5k net lines / 143 commits, but additive and bit-identical to
  the reference model at every step (the reference suite never changed). The C++ AD engines
  live in src/{ff16,tf24,tf24f}_emergent.cpp + inst/include/plant/gradient/ +
  inst/include/plant/models/*_production_kernel.h; the R API in R/emergent_gradient.R
  and R/tf24*_emergent_gradient.R.
• The foundation slice (<T,E,S> templating + exact AD growth-rate gradient) was PR [AutoDiff] Scalar-template the state/params + exact AD growth-rate gradient (#472 scope B, Milestone C foundations + #537 A1) #541,
  now closed and folded in here as branch ancestry.
• The guide is a staging copy destined for the Overstorey site; it builds no figures.

🤖 Generated with Claude Code