Forward-model correctness fixes in TF24's leaf, soil and demographic paths

[aornugent]

Ten fixes to TF24's forward model, plus one that reaches FF16 and K93. Six are bit-identical corrections of latent defects; four move simulated numbers and need baselines re-blessed.

The common thread is state leaking between plants that should be independent. Every cohort of a species shares one Strategy, and so one Leaf. Several of these are the same defect at different depths: a buffer sized but not cleared, a member refreshed by one caller and read by another, a rate computed once and read after the state it described changed. Individually small; together they were enough that re-running one plant's rates from its own inputs did not reproduce what the solver had just computed. Solving this now helps with reverse mode where we want to compute from rates backwards.

The changes

1. Uptake charged to the wrong plant. Leaf::set_physiology sized its per-layer uptake buffer with resize, whose fill reaches only new elements, and the solve writes only as deep as the plant has roots — so layers below rooting depth kept the previously solved plant's values, billed to the water balance. (assign, matching the sibling buffers two lines above.) 33.8% of plant-time records affected, 7.3% holding a stale non-zero value; the gap is the early run, separable only by solving the same states in different orders. The shut-down exits leave the same fields stale (zero incidence at the default driver, margin 27x tolerance); and the profit derivative reads a soil vector only the solve refreshes, differing by 100% when called stale. Addresses [TF24] soil_consumption_ carries over between cohorts #577 and [TF24] Shut-down exits leave the leaf's uptake untouched #578.
2. Uptake integrated from the wrong lower endpoint. Water started at the smallest existing plant and returned zero for a lone plant, so a transpiring recruit drew nothing — at 0.70% of output times, the establishment window. It now starts at the newborn, as the light field already does, and the one-plant case cannot arise. Reaches the recruitment channel, since establishment reads the soil this depletes.
3. Rates not recomputed at introduction. Introduction widens the state and rebuilds the field but left the rates describing the pre-introduction state, which the integrator then used as the derivative of the post-introduction state. Wrong by more than a plant's own magnitude at 51 of 141 introductions (the transport term amplifies). Family-wide, and the only change moving FF16 and K93.
4. Three canopy profiles collapsed to one. TF24 wrote its own cumulative form, density, inverse and a second inlined copy of the cumulative form; it now uses the shared CanopyShape FF16 and K93 already do, and the crown-centre coordinate (three copies) becomes one static function. Compared before switching over 15 087 real pairs: every difference bounded by 2 units in the last place, so the forward move is the step controller, not the profile. develop's exponent-specialised multiplication chains are untouched — an earlier draft rewrote them for readability and moved two models for nothing; they supply values, carry no exponent term, and were never the route to a derivative. One pow stays because the root-mass distribution calls the cumulative form at a different exponent.
5. Leaf-area density finite at the crown base. The density divided by height, giving 0/0 at the base for every height — the field's lowest sample point. It now branches to the limit there: zero above exponent 1, 2/height at 1. The density and cumulative form are an exact derivative pair, verified to hold to 4e-12–6e-11 wherever a difference can check it.
6. Soil moisture inverse converts its units. Missing the MPa→Pa factor, so it returned values scaled by 8.19, above saturation. R-facing utility only, now with a round-trip test.
7. Resource vector sized by resource count. Sized by ODE width (nine) and filled with NA, so four NAs per plant per stage were integrated to NaN and discarded — 1 979 of 2 000 calls, none after. Now an explicit n_resources() on the environment, one source of truth.
8. Newborn leaf solved once per stage, not twice. Setup computed the rates, then establishment recomputed the same quantity at identical arguments (a newborn is at birth size already). The stored value is passed in instead — bit-identical, one fewer solve per stage per species, the total difference equal to the stage count. FF16 fixed the same way; K93 reads no carbon budget.

Measurements, for the whole change

One species, five soil layers, default driver, full patch lifetime, refine_schedule = FALSE, built -O2 -DNDEBUG -g0. Both numbers from one session on one machine — a value gate that does not name its flags measures the compiler, since this tree at -O0 differs by 0.145% in offspring from the step controller alone.

	offspring production	accepted steps
develop	42.140173575095666	5 055
this branch	42.176246845059751	5 105

+0.086% in offspring, +50 accepted steps. The four value-moving fixes measured individually shift offspring by +0.14% to +0.79%; composed, +0.086%. They do not add — most of each individual figure is the adaptive controller landing on a different accepted sequence, not a systematic change. Treat the composite as the number, and distrust any per-fix figure below about 0.15%.

FF16 and K93 move only through change 3: FF16 +0.0034% (two more steps), K93 +0.000038% (same step count), both at full lifetime.

Runtime unchanged. Change 3 adds 141 rate evaluations against ~30 000; change 8 removes one solve per stage per species.

Tests: 1 007 assertions pass across 15 files. Three fail, and they are the re-blessing this asks for — all in FF16's "offspring arrival", on offspring production, an ODE time and a step count, against testthat's default ~1.5e-08 tolerance given a 3.4e-05 movement. Nothing under tests/testthat/FF16_reference/ was regenerated, and FF16's reference-comparison file still passes because its own tolerance absorbs the shift — so the hard-coded assertions are the tripwire, not the reference files. Two pre-existing failures (mutant-fixture tests) are unrelated and identical on develop.

[aornugent]

Flagging an alternative: 49b03bb adds a compute_rates() call to introduce_new_nodes, which is correct and unblocked. But it increases the number of state-mutating paths in Patch that must remember to recompute — set_ode_state (both overloads) and now introduce_new_nodes — because ode_rates() is const and only copies out whatever the last set_ode_state computed. Any future path that reaches a state by another route reintroduces the same bug.

The structural fix is to move compute_rates() out of both set_ode_state overloads and into ode_rates(), which is then no longer const. ode_rates today returns whatever the last set_ode_state computed, so any caller that arrives at a state by another route reads rates belonging to a different one. Patch::introduce_new_nodes widens the state and rebuilds the light field without recomputing rates, and set_state_from_system — the second blocked signature below — then seeds dydt_in from those stale rates and marks them clean. Under first-same-as-last that vector becomes k1. With ode_rates computing, that cannot arise — the rates are always those of the state just set.

No arithmetic moves and no reordering: the forward path calls set_ode_state then ode_rates, so the same work happens in the same sequence. The invariant would hold by construction — rates are always those of the state just set — and the compute_rates() added here becomes unnecessary rather than correct.

It also separates loading a state from evaluating at it, which is what a reverse pass wants: re-establish a stage's state and environment without paying for a rate evaluation the recordings will redo. See aornugent#65 as an example.

It is blocked in odelia, which is why it isn't proposed as part of this PR.
Two signatures take the System by const reference and would need mutable ones:

site	what it does
odelia/ode_interface.hpp — r_ode_rates(const T& obj)	the R accessor
odelia/ode_solver_internal.hpp — set_state_from_system(const System& system)	seeds dydt_in from ode_rates and marks it clean

Happy to raise a paired PR if this is preferred.

---

Generated by Claude Code

[aornugent]

One more in 903b8292.

The mechanism. Patch::reset() returns the patch to t = 0 and calls environment.clear(). For TF24 that reaches clear_environment(), which clears the light profile and never touches vars.states — the five soil moisture states and the four cumulative flux accumulators. So a reset patch is at time zero with the previous run's soil.

It is reachable, and by the one caller that matters. run() begins with reset(), and refine_schedule() calls run() in a loop. So every refinement iteration after the first computes its error signal against a run that started from the previous iteration's depleted soil, and the schedule is then chosen on that signal.

Measured, on the branch before the fix, one species, five layers, default driver, full lifetime:

	offspring production	accepted steps
first run()	42.176246845059751	5 105
second run() on the same SCM	42.296696440471194	5 062

and the state after run() then reset() begins 0.284128746037276 … against 0.214 for a freshly built environment, with the flux accumulators at 81.2 rather than zero. After the fix the second run reproduces the first exactly and the post-reset state is bitwise equal to a fresh one.

The forward run is unmoved — 42.176246845059751 at 5 105 accepted steps, the same figure this branch already reports. That is the safety property rather than a null result: the fix restores an initial condition, it does not change one. Nothing to re-bless.

The starting state is snapshotted, not recomputed, and that is the one design choice here. The constructor's last line is the only initialiser, so recomputing it would reproduce the default — but set_soil_water_state is part of the R interface and is exercised (a test hands run_scm a one-layer soil at 0.1; scenario_eval sets it too). Recomputing the default at reset() would overwrite the caller's choice and move the first run's initial condition, which is exactly what this must not do. So the state is recorded where it is set and restored where it is cleared, and it cannot drift from either.

Generated by Claude Code

[dfalster]

Excellent work catching these @aornugent . Have flagged on for @itowers1 to comment on

[dfalster]

@itowers1 Can you check this?

[itowers1]

Yep, technically true if psi_soil was passed in as MPa, which it emerges from the soil moisture model as. The soil_moist_from_psi functions are not actually in use at the moment from what I can see. Might be worth raising a seperate issue that we need to swap all potential values and parameters to MPa for consistency.

[dfalster]

Thanks @itowers1 — and confirmed on "not in use": the only references anywhere are the RcppR6 binding and the new round-trip test, with no C++ caller, so the fix corrects an exported utility without moving any simulated number. The forward direction on the line above already divides by 1e6 with an explicit Pa -> MPa comment, so this is just its exact inverse.

Agreed the real fix is the consistency pass — that's #587.

---

Drafted with Claude Code

[dfalster]

Approving on the substance, but three things before this merges — the first is blocking.

CI is red and needs the re-blessing committed. The three FF16 "offspring arrival" assertions fail identically on macOS and Windows. I reproduced them locally and bisected: deleting the compute_rates() added to introduce_new_nodes restores develop's numbers to the last digit, so the attribution in the description holds and everything else in the branch really is bit-identical for FF16. I also traced why that fix is needed rather than merely tidy — set_state_from_system copies ode_rates() into dydt_in and sets dydt_in_is_clean = true, so setup_dydt_in never re-evaluates and that vector becomes k1 under first-same-as-last. So: happy for the numbers to move, but please commit them (4.215899, c(11.995204, 16.474988), 293) rather than merging red. The assertion at line 215 still passes on its 1e-4 tolerance while having drifted to 16.8895016065 — worth recording at the same time.

The FF16 figure in the description has the wrong sign. It says +0.0034%. Offspring actually moves down: 0.00074% and 0.00122% in the two-species run, and up 0.00024% in the one-species run. The 3.4e-5 looks like the ode_times[100] shift (3.7e-5), not offspring. Worth fixing so nobody chases the sign later.

Species::consumption_rate assumes a monotone grid. It now integrates over r_heights_rev() with new_node prepended, which is the same assumption #574 had to abandon for the competition profile, where inverted heights are handled by sorting in compute_competition_unordered. The uptake integral has no equivalent fallback, so an inverted grid gives a partially-cancelling trapezium. Not a regression — develop integrated over the same unsorted grid — and CI won't catch it while the dry-start TF24 test is skipped. But this PR extends that function, so it's the natural moment to make the two paths consistent. Happy for it to be a follow-up issue if you'd rather keep this PR closed.

---

Drafted with Claude Code

[dfalster]

Agreed on the direction — ode_rates() computing is the right shape, and yes please to the paired odelia PR. Three things before you cut it.

It won't hold by construction as written. set_ode_state(it, time) builds the field with compute_environment(true); introduce_new_nodes uses (false). Moving compute_rates() into ode_rates() guarantees the rates match the state just set, but not the environment they were computed against — so a future path that reaches a state without refreshing the field reintroduces the same bug one level up. If you want the strong invariant, the environment refresh has to move with it: one sync point on Patch doing both, called by every mutating path.

Please demonstrate the bit-identity rather than assert it. Under the refactor introduce_new_nodes drops its compute_rates() and set_state_from_system → ode_rates() computes at the same point against the same rescale = false field, so the result should be bitwise equal to this branch's head. That's the difference between the FF16 re-blessing surviving the refactor and needing a second one.

There is a cheaper option in odelia alone, and the reason to reject it is the argument you should be leading with. Have set_state_from_system stop copying rates and leave dydt_in_is_clean = false; setup_dydt_in then calls derivs, which re-sets the state and recomputes from it, so consistency holds however the caller arrived. Same class of bug closed, no plant change at all — and no signature change either, because set_state_from_system would no longer call ode_rates and could stay const&. Half the blocked-in-odelia table goes away. But it routes through set_ode_state(it, time) and so rebuilds the field with rescale = true, which moves the numbers again, whereas yours keeps the field introduce_new_nodes deliberately built. That's a stronger argument for your version than the reverse-mode one you're currently leading with.

Cost isn't an objection, for the record: there are five ode_rates call sites in odelia, three of them inside derivs where set_ode_state immediately precedes, so nothing in the RK stages recomputes twice. The extra work is one evaluation per set_state_from_system — the same order as the 141 introductions this branch already pays for.

Two small things while you're in there. Patch::r_compute_rates() already exists as an explicit compute-from-R entry point, so after this there'd be two ways to do it — collapse them. And since we're already pinned to odelia master for the mutant cache hooks (#564), fold the signature change into that pin rather than making it a separate bump.

---

Drafted with Claude Code

[dfalster]

Follow-up opened: #589, the structural version of change 3, plus traitecoevo/odelia#46 for the two signatures it needs. Both drafted against the questions raised above.

The bit-identity is demonstrated, not asserted. #589 is bitwise equal to this branch on FF16 (one and two species), K93, TF24 at full patch lifetime and FF16 schedule refinement — same offspring production to all 17 digits, same accepted step counts. So the re-blessing you land here stays valid through the refactor, which was the thing worth checking. Runtime is unchanged too: 44.26 s vs 44.36 s on TF24 full lifetime over six alternating reps, against a 1.5 s spread within each build.

One thing worth knowing before this merges. Moving the compute into ode_rates() made patch$ode_rates segfault on a fresh R-constructed patch, and the cause is already on develop: Patch holds environment_type* environment_ptr at its own environment member with default copy semantics, so every copy's pointer refers to the source's members. scm.h does patch = sys and RcppR6 hands copies to R. It survives today only because the sole dereference is compute_rates(), which every live path reaches on an object whose source is still alive — the first new caller to dereference it crashes. #589's first commit replaces it with an index. Nothing for this PR to do, but it is a live trap for anything else that touches that pointer.

Also filed #588: Node::compute_initial_conditions computes rates and then seeds mortality and log-density, so the rates it stores describe the state before the seeding — the same shape one level down. Harmless through the solver once #589 lands (the patch recomputes first); it is what makes test-patch.R:86 merely equal rather than identical, at ~1e-22.

---

Drafted with Claude Code