Add a general discrete-event queue

[dfalster]

Adds a general (time, action) event queue. Node introductions move onto it,
and rainfall pulses, thinning and heat damage join them as event types.

An event carries when it happens, its type, its target (environment, patch or
one species) and its values. What each one actually did -- as against what was
asked -- is readable as scm$event_log.

Runs supplying no events are unchanged: FF16, K93 and TF24 are identical() on
ODE step times, fitness and final state.

Closes #628. Unblocks #627 (thinning) and #601 (shared queue). Design notes in
notes/plan-events.md.

[dfalster]

What this builds, and what it deliberately does not

The mechanism already existed; only the payload was missing

SCM::run_next_impl has always stopped the integrator, mutated the patch, re-read the ODE state and resumed — that is exactly the pattern a rainfall pulse needs. What it lacked was a typed payload: every event was a node introduction with a bare species index. This PR gives the queue a type tag, a target and a parameter vector, and routes non-introduction events through Patch::apply_event().

Verification

Bit-identity for event-free runs, captured on develop before the first code commit and compared with identical() on raw arrays (not by eye):

	ODE step times	fitness	final state
FF16 (207 steps)	identical	identical	identical
K93 (240 steps)	identical	identical	identical
TF24 (1055 steps)	identical	identical	first 793 identical, one new slot holding exactly 0

Expressing the same introductions through the new events path is also identical() to the default path, on both a one- and a two-species run — the latter to exercise tied introduction times. Full suite: FAIL 0 | WARN 0 | SKIP 4 | PASS 3023.

The fifth accumulator is pulse-only on purpose. odelia's controller takes its error norm over every state component, so an accumulator with a non-zero rate would join the step-size decision and could move TF24. Held at rate zero it contributes exactly zero. That is why no scientific_version bump and no scenario re-bless are needed — and it was measured, not assumed.

Three things found along the way

• r_set_max_time() read events.back() on an empty list. Undefined behaviour, on the normal path — both make_node_schedule() and node_schedule_default() set max_time before adding any times. It had always read harmless garbage; changing the event's member layout turned it into a segfault. Fixed, with a test.
• A pulse suppresses the continuous infiltration that follows it. Wetting layer 0 raises the saturation-excess term, so a pulsed run gains strictly less infiltration than the pulse delivered (0.0585 of a 0.06 m pulse, in the test). Real behaviour, worth knowing before reading a water budget.
• Pulsed columns end drier than unpulsed ones. K(θ) rises as θ^16.14, so the extra water drains fast and the wetter interval costs more drainage than it stores. Final storage is the wrong thing to test a pulse with; the cumulative fluxes are the right thing.

Two things a reviewer should not read as bugs

• An event is a stop time. Adding one changes the adaptive step sequence, so a run with events differs from one without at solver tolerance, even away from the events. Bit-identity is claimable only when the set of stop times and actions is unchanged.
• Harvest and partial disturbance are not separate types. They are thinning with different selectivity, and harvest() / partial_disturbance() are R names that read better at their own call sites. One action, three ways of asking. [thinning] Implement thinning as an event #627's size-class thinning is the same action again.

Where this departs from #628

• No per-cohort target. A cohort has no stable address across a run: nodes are appended and never removed, and refine_schedule() changes how many exist, so "cohort 7" in a schedule written up front is not well defined. Selecting cohorts is a height band in the action's parameters instead — well defined, and what [thinning] Implement thinning as an event #627 actually asks for.
• heat_damage is the weakest of the actions, deliberately. This model runs its TF24 leaf at a constant 25 °C, so there is no physiological response for a temperature to drive. It is a mortality that rises with heat exposure, computed by sub-integrating over the event's nominal duration at half-hourly steps with demography frozen. The hook is the part worth keeping: the real response belongs in [TF24 hydraulics] Leaf thermal damage, repair & acclimation (ATLS) #566's damage state.

Not in scope

Runtime insertion (designed for, unexposed — #601 needs it for stochastic deaths); state-triggered events, which need rootfinding odelia does not have; the stochastic tower, which still refuses any event but an introduction; and disturbance as a patch-clearing event — Disturbance_Regime is a pure survival weighting that never touches state, genuinely a different mechanism from a partial-disturbance event.

Two hazards to weigh before the next step

1. Step-size inheritance. step_size_last survives set_state_from_system(), so a leg starting after a discrete change begins at whatever step the previous quiet leg grew to. Design pass: the soil cascade needs no change #608 measured this as the trigger for TF24 blow-ups, and pulses create far more leg boundaries than introductions do. Nothing here addresses it.
2. odelia's domain hooks are not adopted. 0.3.1 has ode_state_valid() and util::DomainError; plant pins 0.2.1. A capped pulse still parks layer 0 at exactly θ_sat, where K'(θ) is ≈1.5e4 yr⁻¹.

Both are items from #608, which is still open and unmerged — its design doc lives only on notes/soil-redistribution-design.

Build note

plant develop does not compile against the phylloptim currently installed here (0.6.0 dropped hydraulic_cost_Sperry and changed the Leaf constructor). Everything above was built and measured against the versions plant's own Remotes pins — phylloptim@037673b7 (0.2.0) and odelia@v0.2.1.

[dfalster]

Follow-up: the generic layer now reads generically

@dfalster pointed out a layering error in the first pass: the queue, Patch and the Environment base are shared by every strategy and environment — including, in principle, size-structured animals — but their vocabulary assumed plants, water and heat. Only TF24_* should know about rain.

was	now
RainfallPulse	ResourcePulse
HeatDamage	ClimateExtreme
Thinning	Harvest
height_min / height_max	size_min / size_max

The pulse now rides an abstraction plant already had. Environments declare n_resources(); TF24's are its soil layers. So a pulse names its pool through target_index — exactly as a species-targeted event names a species — and Environment gains one virtual, add_resource_pulse(i, amount). FF16 and K93 declare no resources, so a pulse aimed at them is refused with a message that says how many they have. That is a smaller interface than the water-specific hook it replaces, and it made the "no per-cohort target" argument fall out symmetrically: target_index means which one, within the target, whatever the target is.

Model-specific names stay where they are accurate. rainfall_pulse(time, depth) and TF24_Environment$add_water_pulse(depth) are the same action under the name that reads correctly for TF24; the R wrapper lives in its own file, R/tf24_events.R, so the layering is visible in the file structure and not only in the prose. It builds its rows directly rather than delegating, so a length mismatch is reported against depth — the argument the caller typed — rather than the generic amount.

One consolidation. harvest now covers what was three names: leave the size band at its defaults and it is an across-the-board knock-down, set size_min and it takes everything above a size, set both and it thins one size class (#627). One implementation, and thinning() / partial_disturbance() are gone rather than kept as aliases — three names for one action was clutter.

Verification unchanged

Re-ran the whole ladder after the rename: FF16, K93 and TF24 still identical() on ODE step times, fitness and state (TF24's one extra slot still exactly 0); both roxygen examples run; full suite FAIL 0 | WARN 0 | SKIP 4 | PASS 3023.

[dfalster]

Rebased onto #633, and the evidence re-earned rather than re-asserted

#633 moved TF24 by +0.36%, so every number in this PR's earlier comments was measured against dependencies that no longer exist. Re-ran the whole ladder against the new develop (phylloptim 0.6.0, odelia 0.3.1) rather than carrying the old figures forward.

The claim is unchanged, because it was never a claim about absolute values — it is that this change moves nothing:

	ODE step times	fitness	offspring	final state
FF16 (207 steps)	identical	identical	identical	identical
K93 (240 steps)	identical	identical	identical	identical
TF24 (1065 steps)	identical	identical	identical	first 793 identical, new slot exactly 0

TF24 is 1065 steps now, not the 1055 quoted earlier — that difference is #633's, not this PR's.

Expressing the same node introductions through the events path is still identical() to the default path on both one- and two-species runs, and the schedule still round-trips through Events (141 in, 141 out). Full suite: FAIL 0 | WARN 0 | SKIP 4 | PASS 3015.

Rebase mechanics. Only one file conflicted, R/RcppR6.R, which is generated. Rather than hand-merge it I took one side to get the rebase through and then regenerated with make RcppR6 — which did change the file, confirming the side I took was wrong and that hand-merging generated output would have been a mistake. Both NEWS entries survived intact. develop is now an ancestor of this branch, so it fast-forwards.

[dfalster]

Adopted odelia's domain checks — and a finding about #599

What changed

Two conditions that used to kill a run outright are now rejected steps:

1. A non-finite environment state. Design pass: the soil cascade needs no change #608 measured every observed TF24 soil excursion to be explicit-integrator overshoot — the step a leg inherits across a discrete change, not a defect in the water balance — and the same case integrates cleanly from a smaller step. So it is handed to the stepper to shrink and retry.
2. An infeasible leaf probe out of phylloptim's collar root-find.

⚠️ The second needed a bridge. phylloptim's infeasible_error and odelia's DomainError are both std::runtime_error siblings, not related by inheritance, so odelia's catch (const util::DomainError&) cannot see phylloptim's throw — it escapes and kills the solve having taken zero steps, which is exactly what #608 measured. solve_leaf() now translates it. Deliberately only that one type: a util::stop() from phylloptim still propagates, so a bug stays a bug instead of becoming step-shrinking until "Cannot achieve the desired accuracy".

Two things deliberately left alone:

• A runaway cohort density stays fatal. That divergence is in the equations rather than the stepper, so shrinking cannot recover it and trying would only burn steps before failing anyway with a less specific message.
• ode_state_valid() checks only the environment block. A node's log_density is legitimately -Inf (a cohort that never established), so a blanket finiteness test over the whole ODE vector would reject valid states and stall the solver at its minimum step.

#599 no longer reproduces — so this is insurance, not a demonstrated fix

I ran #599's own recipe (run_stochastic_collect(), TF24, patch_area = 1, default parameters) and measured a control as well, because the develop #599 was written against no longer exists — #633 has landed since.

build	seeds 1–10 that throw
#599's documented develop	6 (seeds 1, 4, 5, 7, 9, 10 of the 17/40)
current develop + this PR, without the domain hooks	0
current develop + this PR, with them	0

So the hooks did not fix #599 — something already had, before them. I am not claiming otherwise, and #599 looks closeable; worth someone confirming across the full 40 seeds and on Windows, since the issue notes it was a platform lottery near the boundary.

That leaves these hooks as insurance against a hazard that is real but not currently biting: #608's step-size inheritance is untouched, and every event is a leg boundary, so pulses create far more opportunities for it than introductions ever did.

Note the scope limit: ode_state_valid() is on Patch. The stochastic tower uses StochasticPatch, a separate class, so only the solve_leaf() bridge (shared through the Strategy) reaches stochastic runs. Extending it is part of #601's dedup.

Verification

Inert where nothing is wrong — FF16, K93 and TF24 still identical() on step times, fitness and final state against a develop baseline captured under the current dependencies. Full suite FAIL 0 | WARN 0 | SKIP 4 | PASS 3038.

The new test pins the contract directly rather than relying on a failure to provoke it: a sane state is accepted, NaN/Inf in any environment slot is refused, and a -Inf log-density in the species block is not.

[dfalster]

Hi @yangsophieee @itowers1 @elijahmagistrado

Here's the event structure we discussed. Feedback, welcome. Also, for general guidance on GitHub workflows, feature branches and reviewing PRs, see new page on overstorey https://traitecoevo.github.io/overstorey/contributing/how-we-work.html

[dfalster]

Rebased onto develop again (#607, #608, #619, #605, #635)

Re-measured rather than re-asserted, since TF24 moved substantially on develop in the meantime: the reference run is now 7350 ODE steps, up from 1065, with TF24@v9. That is #619 and #635, not this PR.

	ODE step times	fitness	offspring	final state
FF16 (207 steps)	identical	identical	identical	identical
K93 (240 steps)	identical	identical	identical	identical
TF24 (7350 steps)	identical	identical	identical	first 793 identical, new slot exactly 0

Events path still identical() to the default path on all three, schedule still round-trips 141 in / 141 out. Full suite FAIL 0 | WARN 0 | SKIP 4 | PASS 3100.

One merge that needed care, and would have compiled wrong

#619 added Patch::ode_state_valid() too — checking something different: that a strategy's declared non_negative_states() (TF24's storage) stay ≥ 0 across the node block. Mine checks the environment block for finiteness. Complementary, but the same function.

git auto-merged patch.h without reporting a conflict, and the result was wrong twice over: my inline definition survived alongside #619's declaration and out-of-line definition — a redefinition that would not have compiled — and my environment check had been silently dropped from the merged body. Resolved by hand: one definition, both checks, environment finiteness first so it is not short-circuited by #619's bounded.empty() early return.

Two more resolved by hand rather than by side-picking:

• Seat TF24 on the TF24_floor cost curve #635 changed TF24_Strategy::solve_leaf() to call leaf.optimise(), not find_root_collar_psi(), and its comment explains why: the latter is hard-wired to CostCurve::TF24 and would silently ignore a seated TF24_floor. Kept Seat TF24 on the TF24_floor cost curve #635's body inside my infeasible_error translation.
• [TF24 hydraulics] Move the TF24 strategies onto the TF24_floor cost curve, with the water price at zero by default #634/Seat TF24 on the TF24_floor cost curve #635 already split TF24f_Strategy::solve_leaf() into a curve dispatcher over solve_leaf_for<K>(), which makes the solve_leaf_impl() split I had added redundant. Dropped mine and wrapped their dispatcher instead.

R/RcppR6.R conflicted again and was resolved the same way as last time — take a side to get through, then make RcppR6, which changed it. Generated files get regenerated, never merged.

Also refreshed

notes/plan-events.md claimed #608's design doc lived only on an unmerged branch. #608 has landed, so it now links to notes/plan-tf24-soil-redistribution.md in the repo, and records that hazard 2 (out-of-domain probes) is now addressed while hazard 1 (inherited step size) is not.

[dfalster]

CI fix: #637 caught a real bug in one of my tests

Error ('test-events.R:129:5'): the events path reproduces the default path exactly
Error: Unknown trait name for K93_Strategy: lma (did you mean eta?)

The events-path guardrail loops over FF16 and K93 but asked for lma in both. K93 has no such parameter — the trait was silently ignored, so the K93 half of that test had been running two identical default strategies rather than the two distinct ones it claimed to compare. #637 started refusing an unknown trait name, and this is exactly the class of thing it was written to catch.

Fixed by perturbing each strategy's own default (lma for FF16, eta for K93) rather than naming a trait literally, so it also survives the defaults moving.

It passed locally only because my build predated #637.

Also merged develop in, rather than rebasing again

@itowers1 had pulled before my last force-push and pushed a merge resolving the collision. That merge's tree is byte-identical to my rebased tip — it added no content, only history. Rather than force-push over it a second time I merged develop in and pushed a fast-forward, so nobody else's checkout gets rewritten again.

This branch now carries #637, #640, #635, #605, #619, #608 and #607.

Evidence re-earned against this develop

	ODE step times	fitness	offspring	final state
FF16 (216 steps)	identical	identical	identical	identical
K93 (240 steps)	identical	identical	identical	identical
TF24 (1219 steps)	identical	identical	identical	first 793 identical, new slot exactly 0

Events path still identical() to the default path on all three; schedule round-trips 141 in / 141 out. Full suite FAIL 0 | WARN 0 | SKIP 4 | PASS 3131.

Step counts differ from my previous comment for two unrelated reasons: TF24 moved again on develop, and the FF16 figure changed because the harness now uses the strategy's default lma instead of the meaningless lma = 1 it had been passing.

One thing worth knowing about the local pin

develop now pins phylloptim (== 0.8.0) — an exact version, not a floor. That is a good tightening given how much drift this branch has seen, and it caught a stale local install immediately (package 'phylloptim' 0.9.0 was found, but == 0.8.0 is required). Anyone building this branch needs phylloptim at exactly 845390c.

[elijahmagistrado]

• I am happy with the general structure of the event system
• the default settings - if rainfall events are enabled, then ensure the user gets a warning that the continuous rainfall is not automatically disabled. or make it so that if rainfall events are enabled, continuous rainfall is disabled automatically
• i think the first implementation of the rainfall pulse should be my proposed system first, then the richards equation comparison can be done later. moisture cascade should be the default system.
• need to keep the actual implementation of how soil water works as a to-do in issue [events] RE rainfall event design and pilot #629, current default of targeting the first layer is a placeholder and likely to chang
• current tests don't include a test of the simultaneous event ordering but it's currently works as expected

There were some issues in the implementation found by GPT 5.6 Sol which should be reviewed by Claude:

---

Event information is lost when collect = TRUE

The event log itself works when run_scm() returns an SCM object:

scm <- run_scm(p, events = ev)
scm$events
scm$event_log

However, run_scm(..., collect = TRUE) returns a different, tidied results object rather than the SCM. That returned object currently contains neither the event schedule nor the event log.

I reproduced this on the current PR with:

p <- scm_base_parameters("FF16")
p$max_patch_lifetime <- 2
p$node_schedule_times <- list()

p <- add_strategies(
  p,
  trait_matrix(1, "lma")
)

ev <- events(
  events_default(p),
  harvest(time = 1, fraction = 0.2)
)

out <- run_scm(
  p,
  events = ev,
  collect = TRUE
)

names(out)

The result was:

steps
n_spp
species
env
offspring_production
net_reproduction_ratios
p

There is no events or event_log element:

out$events
# NULL

out$event_log
# NULL

This matters because events are deliberately supplied separately from p. Therefore, the returned p does not contain the rainfall, harvest or climate schedule. If only the collected result is saved, it no longer records either:

• what events were requested; or
• what those events actually applied.

For rainfall, this includes scientifically important information such as how much water was accepted into the soil and how much was rejected as runoff.

I suggest returning both objects with collected results:

results[["events"]] <- scm$events
results[["event_log"]] <- scm$event_log

Validation should confirm that:

out <- run_scm(p, events = ev, collect = TRUE)

expect_identical(out$events, scm$events)
expect_identical(out$event_log, scm$event_log)

At minimum, the test should verify that:

• the complete requested schedule is returned;
• the log contains every applied non-introduction event;
• event order is preserved;
• requested and applied values are retained;
• reset or repeated runs do not leak records from a previous run;
• an event-free run returns a valid empty log rather than omitting the field.

The current event-log structure can remain unchanged for this PR. The immediate issue is ensuring that the log is not lost from the main collected-output workflow.

Events after max_patch_lifetime fail late inside the solver

The event validator currently checks that event times are finite and non-negative, but it does not check them against the simulation horizon:

p$max_patch_lifetime

I reproduced the problem with a two-year simulation and a harvest scheduled at year 3:

p <- scm_base_parameters("FF16")
p$max_patch_lifetime <- 2
p$node_schedule_times <- list()

p <- add_strategies(
  p,
  trait_matrix(1, "lma")
)

ev <- events(
  events_default(p),
  harvest(time = 3, fraction = 0.5)
)

run_scm(p, events = ev)

The event object is accepted, but the run later fails with:

Error: time_max must be greater than (or equal to) current time

The underlying schedule is inconsistent:

simulation horizon = 2 years
final event time   = 3 years

The event queue advances toward the year-3 event even though the run should finish at year 2. Once the event at year 3 becomes the current event, the queue still uses year 2 as the final endpoint. This effectively leaves the solver with an invalid interval running from year 3 back to year 2, producing the low-level time-ordering error.

This should instead be rejected when the SCM schedule is constructed, before any integration occurs. The required condition is:

[
0 \leq t_{\mathrm{event}} \leq \texttt{max_patch_lifetime}.
]

The error should identify the offending event and the permitted horizon, for example:

Event 142 (harvest) occurs at time 3, after
max_patch_lifetime = 2

I suggest tests covering all three boundaries:

# Before the horizon: valid
harvest(time = 1.9, fraction = 0.5)

# Exactly at the horizon: valid and applied
harvest(time = 2, fraction = 0.5)

# After the horizon: rejected before integration
harvest(time = 2.1, fraction = 0.5)

The after-horizon test should verify that:

• construction fails before the solver runs;
• the error names the event type and time;
• the error reports max_patch_lifetime;
• the same validation applies to every event type;
• vectorised constructors identify the specific offending row;
• simultaneous events exactly at the horizon remain valid.

I would not silently discard events outside the horizon. An out-of-range event may indicate a time-unit error, a truncated simulation or an incorrectly reused forcing record. Ignoring it could produce a plausible but scientifically incomplete result.

More minor issues

• Event-target validation is incomplete. Invalid combinations such as an environment-targeted harvest or a patch-targeted resource pulse are accepted and then silently interpreted as valid actions. Each event type should enforce its permitted targets.

• Simultaneous ordering only partly works as intended. Different event types are correctly ordered—resource pulse, climate extreme, harvest, then node introduction—but events of the same type and time are applied in reverse input order. For rainfall pulses, total water is unchanged, but accepted and rejected water can be attributed to the wrong event record. Input order should be preserved or simultaneous pulses should be combined.

• Event-specific parameters need earlier validation. For example, intensity = NaN becomes a zero-damage climate event, negative sensitivity can produce negative applied mortality, and size_min > size_max silently produces an empty harvest. These should be rejected before integration.

• Adding events in the middle of the run_scm() arguments breaks existing positional calls. Moving it to the end of the function signature would preserve backward compatibility.

• The event-log documentation overstates what harvest records. It reports the requested fraction and number of numerical cohort nodes affected, not the actual density, biomass or number of individuals removed. Either record the removed quantity or narrow the wording.

---

[dfalster]

Thanks @elijahmagistrado — these were all real, and two of them were bugs I would not have found without a user of the interface. Addressed in a6a70f0.

The two significant ones

Event record lost when collect = TRUE. Fixed: results$events and results$event_log are returned. You are right that this is the workflow that matters — events are deliberately supplied separately from p, so a saved collected result recorded neither the schedule nor what it applied, including how much of each pulse the soil took and how much it shed. Tested for the full schedule, the applied log, order, requested-vs-applied, an empty-but-present log on an event-free run, and no leakage across re-runs.

Events past max_patch_lifetime. Now refused when the schedule is built, before any integration, naming the row, its type, its time and the horizon. Tested at all three boundaries — 1.9 valid, 2 valid and applied, 2.1 rejected — for every type, and with the offending row identified among many. I agree with your reasoning for rejecting rather than dropping: an out-of-range event usually means a units slip or a reused forcing record, and silently discarding it yields a plausible answer missing the intervention.

The minor ones

Target validation. Each type now declares every target it accepts, so an environment-aimed harvest or a patch-aimed pulse is refused, with the message naming what the type will take.

Simultaneous ordering. Real bug, exactly as you describe: same-type same-time events were applied in reverse input order, so for two pulses at one instant the accepted and shed water was attributed to the wrong records. The queue's tie-break is stable now. This changed the tie order for same-time introductions too, which two test-node-schedule.R assertions pinned — updated, and the SCM identity checks confirm nothing else moved.

While fixing it I found a related trap: events() sorted by time only, so the R object could list a different order from the one that runs. It now sorts by time and type, and there is a test asserting the R object matches scm$events — a guard against R's idea of the order drifting from the C++ enum.

Parameter validation. intensity = NaN, negative sensitivity or duration, size_min > size_max, negative pulse depth and fraction outside [0, 1) are all refused at construction now.

Argument position. events moved to the end of run_scm()'s signature.

Log wording. Took the "record the removed quantity" option rather than narrowing the wording: harvest and climate extremes now report {fraction_applied, nodes_affected, density_removed}. nodes_affected is documented as bookkeeping about the discretisation; density_removed is what the intervention actually took out.

On the science points

Rainfall pulses alongside continuous rainfall. I have gone with the warning rather than auto-disabling, because silently changing a driver the user set seemed worse than telling them: rainfall_pulse(..., env = env) warns if that environment still carries a non-zero rainfall driver, and says how to zero it. Happy to switch to auto-disable if you would rather — it is a one-line change, and you have the better claim on which is less surprising.

Moisture cascade first, Richards later; first-layer targeting is a placeholder. Agreed, and not attempted here. rainfall_pulse()'s layer argument now carries an explicit ⚠️ pointing at #629 and saying the answer there is expected to change it.

Simultaneous ordering untested. Now tested, both across types and within one type.

Full suite FAIL 0 | WARN 0 | SKIP 4 | PASS 3166; FF16, K93 and TF24 still identical() to develop on an event-free run.

[elijahmagistrado]

I am happy with the most recent changes.