Changelog • odelia

odelia 0.11.1

A step, a replay and a sweep hand a System the same stage times. Each explicit stepper computed a stage’s time as time + c * h in every compiled copy of a step: the run’s, a replay’s and the sweep’s. A compiler free to fuse a multiply-add can fuse one copy and not another. g++ 15 does this on Ubuntu 26.04’s x86-64-v3 default, so a sweep handed a System a time an ulp from any the run handed it, and a System reading a record of the run by time found nothing there (traitecoevo/plant#767). Every explicit stepper’s stage time, and a substep’s time and fraction within a step, now comes from util::time_in_step. It rounds the product on its own through a volatile, so no compiler can fuse it anywhere. A new standalone check, built at -O2 under the compiler’s own contraction default, runs a spring for thousands of steps under both explicit steppers. It requires the run, a replay and the sweep to hand the System every stage time rounded that way. It fails without the fix both on ubuntu-latest and under Apple clang on arm64.

A patch bump. On an x86-64 baseline build, which cannot fuse, no number moves. Where a compiler fused some stage times, as Apple clang does on arm64, those times move by an ulp and an adaptive run’s later steps with them: the standalone spring’s state moves by up to about 4e-12 relative. plant’s runs and gradients on macOS are bit-identical to 0.11.0’s.

odelia 0.11.0

A sweep holds one active System across every width it passes through. Solver::sweep_adjoint rebound the System afresh for each segment and again for each insertion, so a System that widens at most steps paid its own construction twice a step. It now rebinds once and re-seats the same active System at each width (adjoint_system::reseat), where the System lets it. A System that declares assign_from(passive) (AssignsFrom) is written in place; one that does not gets a fresh active System at each width, as before. assign_from must leave the System equal to a rebind in every member a recording reads. A new standalone test sweeps a widening System across two insertions with and without the hook and finds every adjoint bit-identical.

An output adjoint that is a multiple of an earlier one costs no sweep. A sweep is linear in the adjoints it starts from, so sweep_adjoint now carries such a row down as the earlier one and rebuilds it at the end (adjoint_multiples, at 1e-14 of the row’s largest entry). Every other row is what it would be swept alone, bit for bit; a rebuilt row agrees with its own sweep to rounding. All-zero rows are swept as given.

A spline’s span coefficients are double whatever the knot data’s scalar is. A read is linear in the four knot values and slopes around it, so at an active scalar it is now taken in double and the knot data’s partials attached in one expression: one statement a read, inside the knots or past either end, and none a build. The value is the passive spline’s bit for bit. set_data takes its vectors by value, so a caller done with them moves them in rather than recording a copy of each. for_each_active walks only the knot data. A quintic spline at an active scalar is refused at compile time, since its curvatures are not kept.

The tape stores an operation in 12 bytes. XAD_REDUCED_MEMORY is set in the vendored Config.hpp, and OperationsContainer writes a statement’s usual one to four operations without a memmove call (tools/patches/0003-xad-reduced-memory-copy-into.patch). The setting fixes the tape’s layout, so a consumer must be built against these headers and linked against this release’s library.

Together, with plant’s own changes to use them, these take a gradient of plant’s FF16 stand over every trait column (lifetime 60, birth-date coordinate) from 0.211 to 0.142 s for one metric and from 0.374 to 0.233 s for five, and of its TF24 stand at lifetime 6 from 1.65 to 1.54 s and 2.31 to 2.01 s. The runs are bit-identical, and the gradients agree with 0.10.5’s to 1e-13 relative (traitecoevo/plant#677).

New tests: the two standalone sweeps above, and in test-spline-ad.R a build that records nothing, a slope read differentiated against central differences, and read_linear’s weights equal to an active read’s partials exactly.

A minor bump: consumers re-pin, and gradients taken through a spline at an active scalar, or as a multiple of another output, move at rounding. No value moves.

odelia 0.10.5

The reverse-mode scalar has a header of its own, and implicit_node.hpp no longer includes the System contract. adjoint_scalar<T>, adjoint_tape<T>, visit_active and count_active_slots move from ode_interface.hpp to adjoint_scalar.hpp, beside tangent.hpp. implicit_node.hpp includes that header in place of ode_interface.hpp, so a consumer that records partials on a tape and takes no solver, as phylloptim does, now includes nothing an edit to the solver’s interface touches. ode_interface.hpp includes adjoint_scalar.hpp, so every name is still reached where it was. A file that reached ode_interface.hpp’s other names through implicit_node.hpp alone must now include it itself.

A patch bump: no name changes, and no number moves.

odelia 0.10.4

implicit_value holds and restores its inputs’ adjoints without a tape call per slot. The form that sweeps its residual locally reads each input’s adjoint, zeroes it, and puts it back afterwards. It did that through Tape::derivative(slot), an out-of-line call with a range check, four times per slot, and kept the slots and held values in two vectors allocated on every call. It now fetches the tape’s adjoint array once on each side of the sweep, keeps up to 64 slots on the stack, and pushes its partials as one run. A residual that calls implicit_value itself gets its own buffers. On plant’s TF24 drought stand, 40 years and 108 introductions at ode_tol 1e-3, the reverse sweep’s CPU fell from 6.25 s to 6.03 s, median of six interleaved runs, with every gradient entry bit-identical (traitecoevo/odelia#81).

The adjoint array comes from Tape::getDerivativesData(), a local patch to the vendored XAD (tools/patches/0002-xad-derivatives-data.patch), compiled into odelia’s library like the rest of the tape runtime. A consumer built against these headers needs odelia 0.10.4’s library at run time. The three-argument implicit_value also no longer allocates a vector for its one correction. New standalone tests hold any count of inputs, either side of the stack buffer’s size, and a residual nesting a second implicit_value, to exact equality.

A patch bump: the tape receives the same operations in the same order, and no number moves.

odelia 0.10.3

A solver’s parameters can be read and set on the solver. A Solver steps its own copy of the System it was built from, so System$set_params() after building a solver never reached that solver, and nothing said so. The Lorenz and leaf-thermal solvers gain $pars() and $set_params(params), through Solver_params_impl and Solver_set_params_impl in solver_interface.hpp, which read and write the copy’s parameters in the order ad_parameters() lists them. set_params re-seeds the solver at its current state and time, as set_state does, so the rates it carries are the new parameters’ and the step size starts afresh; call $reset() as well to start again from the initial state. The copy is now stated at the Solver class, at Solver_new, and in the solver constructors’ R documentation; the fitting vignette integrates at the fitted values on the solver that fitted them instead of building a second one. A test pins both halves: parameters set on the System after the solver is built do not move its trajectory, and parameters set on the solver give the trajectory of a solver built at those values.

A patch bump: additions only, and no number moves.

odelia 0.10.2

A sum of many spline reads can be recorded as one tape statement. A model that integrates a function of a field over quadrature points reads the field at every point, and each read at an active scalar records statements that every sweep walks again. Three additions let such a model compute the sum in double and supply its partials once:

hermite_spline<S>::read_linear(u) returns a cubic read as the knot data it is linear in: the read’s value and slope, the four knot values and slopes it depends on (input, pointers into the spline), and the weight on each. The value and slope equal eval(u) and slope(u) bit for bit, beyond either end too, in a build without floating-point contraction; with it, the end line can differ from an active spline’s eval by an ulp.
value_and_partials(f, x...) evaluates f once at a forward-mode scalar carrying one direction per argument, and returns the value and every partial, so a small integrand written once over its scalar needs no hand-written derivative.
supply_derivatives also takes its partials as two parallel arrays, input pointers and derivatives. A caller that supplies a hundred partials a call can fill storage it reuses, which a span of input_and_derivative, holding references, cannot be.

supply_derivatives now pushes partials to the tape from a stack buffer in batches of 64, rather than one call per partial. The tape receives the same operations in the same order; on plant’s FF16 gradient, one call per partial cost 4% more instructions.

plant uses these to record each crown integral as one statement (traitecoevo/plant#677, traitecoevo/plant#692). New tests in test-spline-ad.R hold read_linear’s value and slope to eval’s bit for bit inside the knots, on a knot and past each end, built with -ffp-contract=off, and its weights to the tape’s gradient to 1e-14; check value_and_partials against a function’s exact partials; and, in test-implicit-value.R, supply 150 partials through the parallel form, with zeros and a passive input among them, and refuse them all when one past the first batch is not finite.

A patch bump: additions only, and no number moves.

odelia 0.10.1

A System can scale each component’s tolerance. A System that declares tolerance_scales(time, scales) (ScalesTolerances) fills scales with one positive, finite multiplier per component for an adaptive step starting at time, and the step multiplies each component’s error level by it, so a scale above 1 loosens that component. The scales are read once at the step’s start and held across its retries. OdeControl::error_ratio() and adjust_step_size() take them as an optional span, empty scaling nothing, so they apply under both step-size rules and in the test of whether the stepper’s error on a subsystem set the step: a loosened subsystem is not substepped for an error the step no longer reads. The scales decide which steps are taken and nothing else: a pinned run, a replay and the sweep take the steps as recorded. A scale that is not positive and finite, or a vector that is not one per component, is refused.

A subsystem’s substeps take their tolerance from the control. OdeControl::substep_tol (default 1e-9, the fixed value it replaces) is the relative tolerance that chooses a subsystem’s substeps; the absolute is 1e-4 of it. Loosening it saves little, because the substeps evaluate only the subsystem’s rates: on TF24’s soil, 1e-7 halved the substeps and left run time unchanged.

A patch bump: a System that declares no scales steps as before, bit for bit, and the default substep tolerance is the value it replaces.

odelia 0.10.0

The reverse-mode names say what they name (#66). A System on another scalar is rebind<U>(), not rebind_from<U>(). The adjoint scalar and the System lifted onto it are adjoint_scalar and adjoint_system, beside tangent_scalar. The sweep is Solver::sweep_adjoint(), not solve_adjoint(), and the rows it sweeps at one state width are a segment: state_at_segment(), and AdjointOverflowError (util::stop_adjoint_overflow()) for an adjoint that left the range of a double. A row of a schedule is a schedule_row, not an instruction, and schedule_from() builds one from a recording, replacing program_from(). Solver::schedule() is steps_taken(), because it leaves out the insertion rows. The two advance_recorded() overloads are advance_schedule(schedule), which steps a schedule and solves its stages afresh, and replay(recording), which loads what the recording’s stages solved for. record_with_derivatives() is supply_derivatives(), returning a supply_report, and value_with_slope<T> is with_slope<T>, in with_slope.hpp. In R, $fit() is $loss_and_gradient() on both example solvers, its ic argument is initial_state, and $set_target() takes schedule rather than times, since it is the step times replayed and not the times observed. Nothing else changes, and no number moves.

Adjoints, directions and partials have one word each (#66). adjoint_rows is adjoint_batch: one adjoint per output, size() adjoints of width() entries, with adjoint_batch::one(), adjoint_batch::identity() and to_vectors() for one_row(), all_rows() and to_rows(). A “row” is now only a row of a schedule or a recording. seed_direction() is set_direction(), beside derivative_along(), and “seed” is gone from the library: what a sweep starts from is each output’s adjoint, and what a tangent carries is a direction. What supply_derivatives() attaches to a value is one partial per input. “Replay” now means only replay(recording), which loads what a recording’s stages solved for; a run on a schedule’s steps that solves its stages afresh is a pinned run, and the messages that told a caller to “replay the run’s schedule” now say to step it. An evaluation’s place in its step’s row is its stage: end_slot is end_stage, and the slot argument of store_solved, load_solved and derivs is the stage. Cash–Karp’s stage at t + h, whose inputs a substepped step reads, is the last stage: last_stage, inputs_at_last_stage.

A minor bump: every consumer re-pins once for the names. No number moves.

odelia 0.9.1

A Dormand-Prince step substeps a System’s subsystem as a Cash-Karp step does. The substepped step is SubsystemSubsteps (ode_substeps.hpp), which both explicit steppers derive from on their own tableau: method = "dopri" takes, replays and sweeps substeps as "rkck" does, and is past the stability limit on the subsystem where it times the step exceeds 3.306, the end of Dormand-Prince’s real stability interval. Dense output inside a substepped step is the stages’ for the rest of the state and no better than the predictor’s for the subsystem, whose substeps are not kept. A System with a subsystem run under "dopri" or "rodas" read the Cash-Karp stepper’s empty error vectors after every step; the solver now asks the stepper the method names, and "rodas" takes the subsystem with the rest.

Whether the stepper’s error on the subsystem set the step is asked in the control’s norm. The subsystem’s part of the error, the rest’s and substepping’s coupling error compare as the rule reads them (OdeControl::error_ratio: the largest ratio under gsl, the RMS under hairer), and the subsystem set the step where its part exceeds both. Under gsl that is the component-wise test as before, bit for bit; under hairer the sums of squares compare, so a subsystem of a few components among many sets the step only where its error is a larger part of the whole than the rest’s.

A patch bump: error_ratio and ode_substeps.hpp are additions, and a run under gsl, or of a System without a subsystem, steps as before.

odelia 0.9.0

Each rate evaluation is handed its slot of the step’s row. store_solved and load_solved take the slot as a second argument: 0 to 4 for the stages after the first, and end_slot (5) for the evaluation at the state the step ends at. A System that reads another run’s recording finds the same evaluation there by its slot and its time. The time alone does not identify it, because Cash–Karp’s stage at t + h and the step’s end share a time but not a state. derivs takes the row and the slot rather than one slot’s values, and carried_rates holds the row.

A minor bump: every System that declares SolvesForValues must take the slot. No number moves.

odelia 0.8.1

A Cash-Karp step can integrate a System’s subsystem in substeps of its own. A System that names a subsystem (HasSubsystem: a range of its state whose rates read the rest only through one input per component) has it substepped on any step that spans one of its forcing’s knots, that is past Cash-Karp’s stability limit on it, or that follows a step whose size Cash-Karp’s error on the subsystem set and substepping’s coupling error would not have. The stiffness is read off the step before, as the subsystem’s rates’ change over its state’s change between the stage at t + h and the end (Step::subsystem_stiffness), and the step is past the limit where it times the step exceeds 3.735, the end of Cash-Karp’s real stability interval; a step whose two states coincide keeps the last measure. The coupling error is the one a substepped step measures and, on a step that did not substep, half the step times the inputs’ miss at t + h of the line along the input slope. OdeControl::substep_subsystem (default true) turns off all but the knots’ substepping. The System no longer reports whether its subsystem is weakly coupled: subsystem_weakly_coupled() is gone from HasSubsystem. The rest takes Cash-Karp’s step; the subsystem is integrated over substeps chosen at 1e-9 under its inputs extrapolated along their slope since the last step (the predictor, which the stages read), then again under the line from the start’s inputs to those at the stage at t + h (the corrector, which is the step’s end). The larger of the passes’ gap and the stage inputs’ defect against the line is the subsystem’s error. The substeps land on every stage fraction and every knot, so a step need not stop at a knot to read the forcing between its stages.

A step that substepped records its input slope and substep ends on its row (instruction::subsystem), and every other pass takes them as recorded: a program’s replay, a recording’s replay, a replay at another scalar and the sweep. So a replay at other parameters takes the run’s substeps, and the sweep differentiates the step the run took. A grid of times alone (step_to) is refused across a step that would substep, and a recorded row without substeps across a knot of the forcing, since neither holds the substeps and the subsystem taken with the rest can run past its stability limit, finitely. Only method = "rkck" substeps; "rodas" and "dopri" refuse a row that substepped. drivers::Function::knots(from, to) and Drivers::knots give where a driver changes from one cubic to the next.

A patch bump: instruction gains a defaulted field, step_by takes the row it steps by, and a System without a subsystem steps as before, bit for bit.

odelia 0.8.0

Hairer’s step-size rule is the default, stabilised as dopri5.f has it (#64). OdeControl now defaults to Controller::hairer ("hairer" in ode_solve() and OdeControl$set_controller()): the error in the RMS norm over components, accept at or below 1, and rescale after every accepted step, damped by the previous step’s error with dopri5.f’s beta = 0.04 (set_stabilisation(), R: OdeControl$set_stabilisation()). The gsl rule’s dead band repeated a step that landed near the edge until it was rejected, and never grew a step whose error sat inside the band. On Lorenz to t = 100 at 1e-6 under Dormand–Prince, the gsl rule rejects 799 attempts for 32000 evaluations; the default rejects 210 for 27600, against deSolve ode45’s 88 and 28200.

The stabilisation matters at an explicit stepper’s stability limit. Without it (set_stabilisation(0), the classical rule and deSolve’s rk), the step grows past the limit, is rejected and grows again: 310 rejections where the stabilised rule makes 1 on a stiff scalar problem under Cash–Karp at 1e-3. It can also grow until the stages cancel and a step returns to where it began, with an error estimate small enough to accept, so the solution stalls short of the equilibrium it should reach. The classical rule takes fewer steps on smooth problems whose step is set by accuracy rather than stability, and plant uses it for that reason.

The rule’s memory of the last step is the solution’s alone. reset(), and so set_state(), clears it, so a solver reused for a second run steps as a new one does. A step the solver refuses after the rule accepted it, or one cut short to land on a requested time, leaves it as it was. A state with no components reads as zero error, so a run that starts empty and widens at an insertion steps on from it. An accepted step, which the stabilisation can shrink, stops at step_size_min. A rejected step shrinks by 0.9 r^(-1/ord), slightly harder than dopri5.f’s stabilised exponent; the two agree at set_stabilisation(0).

A minor bump, because every adaptive step sequence changes. set_controller(Controller::gsl) restores the old sequences exactly. Because the RMS norm averages the error over components, a System with a large state whose error sits in a few components gets a looser bound on each component than the same tolerance gave under gsl, and should tighten its tolerance to keep the same accuracy.

odelia 0.7.5

A cubic read at an active position is one tape statement. hermite_spline<S>::eval at an active query took the value and the slope through the span’s Horner loops and then added the query’s derivative, and every named intermediate on that path was a statement: eight per read. A model that integrates over a field reads it at every quadrature point, so in plant’s FF16 those reads were a large share of a sweep’s tape (plant#677). Inside the knots a cubic’s read is now one expression (Span::value_at_query), with the same operations in the same order, so the value and its derivative with respect to the position do not move; the quintic and the end extensions are unchanged. A knot’s partial is now one product taken from the root of that expression rather than a chain through the Horner steps, so an adjoint reaching the knot values agrees with before to rounding and not bit for bit (up to 1e-11 relative through natural_slopes); a unit seed hides this, since then the two products are equal. A test in test-spline-ad.R reads at an active position inside a span, on an interior knot, on each end knot and past each end, under a seed other than one. It holds the value and the position derivative to the composed read bit for bit and the knot adjoints to it within rounding, checks the knot adjoints against central differences, and counts one statement for the interior read.

A patch bump: no signature changes, and no value moves; a gradient through the knot values moves at rounding.

odelia 0.7.4

A Dormand-Prince run can now be recorded, replayed and swept (#65). Only method = "rkck" recorded a run, and "dopri" refused a sweep and a replay, so a run that wanted dense output and a gradient could not take both from one stepper. Dormand-Prince is explicit like Cash-Karp, and everything the sweep does transfers; what it lacked was its own row and transpose. DopriStep::step now takes the row a step solves into, or loads from where it is handed a const one, and DopriStep::step_adjoint records one step and sweeps it once per seed on the Dormand-Prince tableau. The forward step and that recording build each stage through one set of kernels, so the sweep transposes the arithmetic the step performed, and a dopri run’s forward numbers are bit-identical to before. The two agree to the last place where the compiler does not contract a * b + c into a fused multiply-add; under clang’s default it does so for doubles and not through the adjoint scalar’s operators, and a stage state can then differ by an ulp. Cash-Karp’s kernels have always shared this, and the header comments now say so.

The row stays six entries wide, not seven. Dormand-Prince makes the same six new evaluations a step as Cash-Karp: five stages, and the one at the state the step ends at, which first-same-as-last hands the next step as its first stage. That first stage is carried or derived afresh and never has an entry of its own under either stepper, so step_record keeps its shape and solved[5] still means the evaluation at the step’s end. The seventh evaluation is not recorded by a sweep, because the state the step ends at does not depend on it; only the error estimate and the dense output do.

Every standalone sweep and replay test now runs under both steppers: a sweep against a central difference of the same steps, across an insertion, from an empty start and over empty steps; a subdivided row refused; a replay that keeps the rows it loaded; a carried first stage loaded at its own instant; and a re-seed that derives it afresh. A further standalone test replays a recorded adaptive Lorenz run bit for bit, sweeps it against a central difference in a parameter and an initial condition, counts six recorded rates a step, and checks that keeping the record leaves dense output unchanged bit for bit. The R tests of one step’s transpose, of the stage states the step and its recording build, and of whole-run sweeps loop over both steppers.

A recording replays only under the stepper that made it. The two rows have one shape but not one meaning: entry i is each tableau’s own stage i + 2, and Dormand-Prince puts stages 4 to 6 at other times and states than Cash-Karp. A Cash-Karp recording replayed under Dormand-Prince, or the reverse, loaded half of every row at the wrong stage with every number finite. Each step_record now names its Method, and advance_recorded refuses a row the solver’s own stepper did not take. A program (schedule()) carries no rows and replays under either. Solver::method() reports the stepper a solver was built with.

$fit() from R now replays under the solver’s own method, so a Dormand-Prince solver’s loss and gradient are Dormand-Prince’s. A "rodas" solver’s schedule is still replayed under "rkck", since the Rosenbrock stepper cannot be swept. The R test that holds each stage state the recording builds to the one the step built is compiled with -ffp-contract=off and holds both steppers bit for bit: under contraction the Cash-Karp half passed only because its one input happened to show none.

A patch bump. solved_row and carried_rates move from Step to namespace scope beside step_record, and Method moves from ode_solver_internal.hpp to ode_interface.hpp, which every consumer reaches through the solver; step_record gains a defaulted field, and no Solver signature changes. "rodas" still refuses a sweep and a replay.

odelia 0.7.3

A sweep crosses steps that have no state. A run that starts empty and steps before its first insertion – plant’s invader, whose first cohort comes after the replay starts – recorded steps of width zero below the insertion, and the sweep transposed them, refusing the empty seed (“a seed must have at least one entry”). A step with no state has no output, so it carries no adjoint and adds nothing to any parameter’s; a range at width zero is now passed over, as a range with no step in it already was, and so is an insertion that leaves the state empty (a replay records the resident’s introductions for an invader with none of its own there). A standalone test steps twice with no state around an empty insertion, inserts a cell, and checks the sweep against a central difference. A patch bump: no sweep that ran before moves.

odelia 0.7.2

A replay that keeps its states now keeps the rows it loaded. advance_recorded(rec) hands each stage the row an earlier run solved for, and a replay asked to keep its states recorded its own states beside an empty row, so a sweep of the replay’s recording re-derived at every stage what the replay had been told. The replay’s recording now carries the rows it loaded, and a sweep of it loads them the same way: the loaded values are exogenous to that sweep and carry no derivative. This is what a frozen-environment derivative is – plant’s invader against a recorded resident. The copy is the price of a recording that stands alone, and only a replay asked to keep its states pays it. A standalone test replays a run through a different parameter and checks that the replay’s recording holds the original’s rows and that its sweep matches a central difference of the replay to 1e-8.

And a sweep now takes each step’s first rates the way the run took them. First-same-as-last hands a step the sixth evaluation of the step before, and the sweep re-derived it with no row, loading rows for the other five stages only; a System that places an inner solve from its row, or stands in a loaded field, re-solved or rebuilt at one stage in six, and the replay test above came back 0.9% low. Each step row now records where its first rates came from (step_record::k1_time): the time of the evaluation they were carried from, or NaN where the run derived them afresh – at its start, after an insertion, and after any set_state_from_system() between steps, which keeps no row. Where they were carried, the sweep loads the step below’s sixth row at that time, which is fl(t + h) and can be an ulp from the time the row below records (a replay, or a final step set onto its interval’s end); a System that checks the instant of what it loads, as plant’s Patch does, would otherwise refuse it. Where they were derived afresh the sweep derives them afresh, with no row, so a System that stands in a loaded value across an insertion supplies it there itself. Two more standalone tests check both: a carried row loaded at its own instant on a step that lands an ulp off, and a re-seed between steps that loads no row.

A pinned interval whose sum rounds short of its end is one step again (#72). advance_fixed asked whether a sub-step was its interval’s last by adding it to the time, and 0.2 + (0.9 - 0.2) is 0.8999999999999999: the interval was crossed as a step short of 0.9 and an ulp-wide one onto it, recorded subdivided, and refused by any sweep or replay although nothing had been refused. A sub-step as wide as what remains of its interval is now final whatever its sum rounds to. An interval this affected now takes one step rather than two, so its values move in the last digits; a run whose sub-steps were refused is subdivided as before. A standalone test pins {0, 0.2, 0.9, 1} and checks that no row is subdivided, the sweep matches a central difference, and a replay lands where the run did.

A patch bump. step_record gains a defaulted field, Step::step_adjoint takes the carried rates in place of a bare row (the one-step form without them is unchanged), and a replay that keeps no states is untouched. Every sweep changes where it takes a carried k1: in the row, for a System that solves for values, and at fl(t + h) rather than t for any System, which moves a time-dependent System’s gradient in the last digits.

odelia 0.7.1

A run that starts with no state at all can now be recorded and swept. The solver fills the run’s first recorded row from the next seeding while that row’s state is empty, and a System whose state is empty for real – plant’s FF16 patch before its first cohort, with no environment state of its own – was re-filled from the seeding its first insertion makes. The run’s start was then recorded one cohort wide, the insertion row duplicated it, and solve_adjoint could not reconcile a one-cohort state with the zero cohorts the schedule implies at the start. Whether the first row has been filled is now a flag rather than a test of its length. state_and_parameter_adjoints then accepts an empty state, since the transpose of that first insertion is a map from the parameters alone; it refuses only where the state is empty and the System lists no parameters either. A standalone test sweeps a run from an empty start against a pinned central difference.

A patch bump: no signature changes, and a run whose first state is non-empty records and sweeps exactly as before.

odelia 0.7.0

A recorded run can now be differentiated in reverse mode, so one solve yields the derivative with respect to every parameter at once. Forward mode answers for one parameter per solve; reverse mode answers for all of them, but it must first record every operation the solve performed, and a run of several thousand adaptive steps performs far more of them than memory will hold. The solver stores its state at each accepted step and replays one step’s arithmetic at a time while the derivative is taken, so what is held is bounded by a single step rather than by the run. Memory then grows with the number of steps, which is cheap, rather than with the arithmetic, which is not.

A value a submodel solved for rather than computed is lifted onto the record carrying a derivative obtained by other means — implicit_value and record_with_derivatives — which keeps an iterative solver’s own iterations out of the answer. What selects rather than moves (a step size, a knot position, an arm) stays a plain double and is replayed, because differentiating a selector manufactures a discontinuity the model does not have.

A forward pass can now REPLAY a recording instead of only taking one. The store/load channel already picked its direction by the constness of what a walk handed over, and step_adjoint already handed a const row — but every forward entry point funnelled through one place that zeroed its scratch and passed it mutable, so nothing on the way forward could load. Step::step’s row parameter is templated, so a const row loads and a mutable one stores with no change to the body; advance_recorded() gains an overload taking a recording rather than a program, pairing each step with its own row so the two cannot be crossed. No new concept and no new name. method='rodas' refuses rather than silently re-deriving: the Rosenbrock stepper keeps no per-stage row.

What this is for: a pass that must take the run’s answers rather than its own. The first consumer is plant’s invasion run, where an invader stands in a resident’s recorded field — exogenous to it, so its derivative is zero rather than severed.

⚠️ The recorded row is six long, not five, and that changes step_record. Five of the six are a step’s stages; the sixth is the evaluation at the state the step ends at, which first-same-as-last hands the next step as its own k1. A sweep re-derives that one at the state it was handed and still reads only the first five. A forward replay cannot — re-deriving is the thing it replays to avoid — and a step whose k1 was re-derived is wrong at first order in h.

The forward-mode scalar lives in tangent.hpp, apart from the reverse-mode machinery in adjoint.hpp, so a consumer wanting a directional derivative and no record is never handed vocabulary for one. tangent.hpp rejects a forward scalar nested above a reverse one at compile time: three kernels costing 31 statements flat cost 566 nested.

ode_fit.hpp is removed; compute_gradient has no drop-in, and vector_jacobian_product plus sweep.hpp replace it in C++. The R fitting interface is kept on the new machinery: $set_target(times, target, obs_indices) and $fit(ic, params) return the same least-squares loss and gradient as before, now by one reverse sweep over a replay of times, for both the Lorenz and the leaf thermal examples. What goes is the active argument on the Solver_* bindings and the separate AD solver it made, since every solver can now fit. A System needs set_recorded_state() and for_each_active() to be swept. The System contract becomes C++20 concepts (HasOdeTime, SolvesForValues, ChecksState, Rebindable), and rebind() becomes rebind_from().

A pinned step the domain refused is crossed in several sub-steps (0.4.0) and recorded as one row of the interval, holding the last sub-step’s solved values; one Runge–Kutta step of that interval is not the step the run took, and a sweep of it gave a derivative forty orders of magnitude off with every number finite. The row now carries subdivided, and solve_adjoint and a replay refuse it by name. solve_adjoint also refuses one batch handed as both seed and accumulator, which the width check let through whenever the state and parameter counts matched. The variadic implicit_value returned rows that included an input’s pre-existing adjoint; they are now the theorem’s alone.

Only method = "rkck" records a run: "dopri" and "rodas" refuse a sweep or a replay rather than re-deriving one.

A consumer that LinkingTo odelia now compiles as C++20 (CXX_STD = CXX20) and must set the same two XAD defines as src/Makevars, -DXAD_NO_THREADLOCAL -DXAD_USE_STRONG_INLINE: the storage class of the active tape does not change the symbol’s name, so a translation unit built without them reaches the same tape through the other storage class.

A minor bump: rebind() becomes rebind_from(), the trait probes become concepts, ode_fit.hpp and the active argument on the Solver_* bindings go, and step_record changes shape, so a System or a consumer that touched any of those recompiles differently. 0.6.0, 0.6.1 and 0.6.2 shipped between this change’s first draft and its merge; the reverse-mode machinery sits on top of all three.

odelia 0.6.2

Resizing a solver no longer costs the square of the state length when RODAS is not the stepper (plant#656). SolverInternal::resize() resizes every stepper the solver holds, and RodasStep::resize() zero-filled its two dense n×n buffers (J, W) each time, whatever the method. A consumer that changes its state length often paid that on every change: plant grows its state at each cohort introduction, so a run over n unknowns cost O(n²) per introduction for a stepper it never calls. The buffers are now sized on RODAS’s first step at a new length, and W without a fill, since every step overwrites it. Results are bit-identical; only allocation moved. In plant, run_mutant() over 51 FF16 mutants against one resident fell from 45 s to 0.72 s (63×), now cheaper per mutant than running them one at a time, and a three-species FF16 resident run from 0.76 s to 0.37 s.

A patch bump: no signature changes, and the same results for every method.

odelia 0.6.1

A fixed point, its stability, and how it moves with every parameter, without integrating to it (#39). SteadyState<System> in ode_steady_state.hpp solves f(y*) = 0 by damped Newton with an Armijo line search, reusing the implicit stepper’s Jacobian and dense LU; sensitivity() is the implicit-function-theorem dy*/dθ = −(df/dy)⁻¹ (df/dθ) from the same factorisation, one linear solve per parameter and no pass through the transient; eigenvalues() and spectral_abscissa() read attraction off df/dy (Hessenberg reduction, then Francis double-shift QR, in ode_linalg.hpp). df/dθ is Jacobian::compute_params(), the forward sweep that already gives df/dy with its unit seed on a parameter instead of a state component; the System names those parameters through ad_parameters(), a vector of pointers in a fixed order, which is the hook the reverse-mode work (#59) uses for the same purpose. Newton finds any root, so solve_with_warmup() is the call that seeks an attractor: it integrates the transient first, RODAS for a stiff system, when Newton from the given guess does not converge or converges to a root that is not attracting, the trivial equilibrium of a demographic model being the usual case. A Newton step that is not finite or does not reduce the residual at any fraction down to min_lambda stops the solve unconverged rather than being taken; sensitivity() and the stability queries refuse until a solve has converged. Fixed-point attractors of autonomous systems only; a system that declares ode_autonomous() is not asked for df/dt, and time_dependence() reports it for one that does not.

The same from R: ode_steady_state(). A right-hand side in the shape ode_solve() takes, func(t, y, parms), is solved to its fixed point by the same Newton on the callback System of 0.6.0, with the Jacobian from jacfunc or forward differences; the result carries the Jacobian, its eigenvalues and the stability verdict, and warmup = <span> integrates the transient past a repelling root. There is no AD through R code, so the sensitivity dy*/dp is a central difference of func in each numeric element of parms, two evaluations per parameter, put through the implicit function theorem with the returned Jacobian in R. A function written for deSolve or rootSolve::steady() runs unchanged.

The ad_parameters() contract is that the rates read each parameter live. A quantity the constructor derived from a parameter carries no tangent, so that parameter’s sensitivity column is zero with no error. check_parameters() compares the AD columns with a central difference taken through rebind(), which rebuilds the System and so does see the derived quantity; a breach reads as a discrepancy of order one, an honest System as finite-difference noise.

Asserted on a two-state system with a closed-form equilibrium (dy*/dθ within 1e-8 of the closed form and 1e-6 of central differences of the re-solved equilibrium, eigenvalues within 1e-9 of their analytic values), on a logistic system whose trivial root Newton reaches from a small guess and whose interior attractor the warm start then finds, on a cached-parameter System that check_parameters() flags, and on the QR routine itself across complex pairs, a companion matrix, a Jordan block and 4×4 and 6×6 cyclic permutations, all to 1e-14.

A patch bump: a new header, additive, no solver signature changed, so nothing that LinkingTo odelia recompiles differently.

odelia 0.6.0

A right-hand side written in R can be solved by the same steppers (#62). ode_solve() is shaped like deSolve::ode(), so a function written for it runs unchanged, with any stepper, and returns a matrix with a time column; OdeSolver underneath is driven a step at a time and lets a caller change the state between steps, at a different length if need be, which is what a consumer with events needs. domain_error() called anywhere below a callback leaves it and has the step rejected and retried smaller, the #55 bargain reached from R; any other error reaches R as itself. ode_counts() reports what a solve cost in right-hand-side evaluations and Jacobian formations, which for a right-hand side that is itself a large computation is the whole cost of the integration.

The thing that steps a run-time callable is CallbackSystem in the R-free core (ode_callback_system.hpp), a System over std::functions, exercised as plain C++ in tests/standalone/; the R adapter in src/ is the few lines that wrap Rcpp::Functions, and a Python binding would be the same lines over a py::function. Evaluation is lazy, so the three places the solver sets a state without wanting its rates cost nothing.

A third stepper, Dormand–Prince 5(4), with dense output of its own order (#24). method = "dopri" is the pair behind deSolve’s ode45, and what it has over Cash–Karp is a free fourth-order continuous extension from its own stages: Solver::advance_collect(times, dense = TRUE) (R: OdeSolver$advance_collect(), and the default in ode_solve() under "dopri") lets the controller choose its steps and reads the requested times off the step that spans each, so ten thousand output rows cost the same integration as two, at the step’s own accuracy — asserted on a quartic, which the interpolant reproduces to rounding, and on Lorenz, where dense rows sit within 5e-7 of landed rows at 1e-8. Under the other two steppers the interpolant is cubic Hermite on the step’s endpoints, one order short (40–100× the integration error on Lorenz, measured), so with those ode_solve() defaults to dense = FALSE, which lands a step on every requested time as ode45 does. The solver also exposes interpolate(), get_step_size() / set_step_size(), a bounded single step(time_max), a rejection counter and a mid-step flag.

Three additions to the core open the implicit stepper to a callback system and cost every RODAS user less. A System may supply its own Jacobian through an optional ode_jacobian(y, t, dydt, J) hook, taken in preference to the AD route when both exist (fd_jacobian() is the one-line finite-difference body for a system with no rebind()); a System may declare ode_autonomous() and is then not asked for a df/dt term, which was one extra evaluation per step; and the Jacobian and df/dt formed at a step’s start are kept across a retry of the same step rather than evaluated again for the same numbers. RODAS also now declares first-same-as-last, since the derivative it computes at the new point is exactly the next step’s start — it had been recomputed on every step. A singular W is a step rejection rather than a fatal. All of this is additive and bit-identical for a system that declares nothing new, asserted by test-rodas.R unchanged; per accepted RODAS step the budget is six evaluations, plus one for df/dt unless autonomous, plus the Jacobian once.

Measured against deSolve with the same R right-hand side, installed build, median of seven. Per call the two are the same: on Lorenz 0.95 µs here, 1.03 µs there, the R function itself 0.6 µs of each — what remains is R’s. What differs is how many calls a fine output grid costs. Lorenz to t = 100 at 1e-6 with 10001 rows: ode_solve() 27 ms for 32011 evaluations, deSolve::ode(method = "ode45") 57 ms for 60002, lsoda 31 ms for 22729, so 2.1× ode45 and level with lsoda; against a 1e-12 reference over the first two time units the three err by 5.8e-5, 5.6e-6 (its landed steps over-resolve) and 1.3e-4. With 1001 rows all three take 25–27 ms: the integration itself, not the output, is then the cost. Compiled Lorenz through Lorenz_Solver takes 1 ms, a hundred times faster than any of them. Stiff Van der Pol (ε = 1e-4, t = 2): RODAS through R 823 steps and 7449 evaluations, the explicit steppers above 10000, deSolve radau 364 and 3279; odelia and radau agree to 5e-6.

The step-size rule is a switch, with today’s rule the default (#64). OdeControl gains set_controller(): "gsl" is the rule odelia has always had (GSL’s standard control: accept below 1.1, grow only below 0.5, otherwise keep the step), "hairer" is the classical rule of dopri5.f (without its Lund stabilisation term) and deSolve (RMS error norm, accept at 1, rescale after every step, no growth right after a rejection). The dead band in the gsl rule makes it ride the acceptance edge: on Lorenz to t = 100 it rejects 794 attempts where the hairer rule rejects 491 and does the same 4700 attempts and 28200 evaluations deSolve’s ode45 does. The default stays gsl because changing it changes every adaptive step sequence plant takes; ode_solve(controller = ) and OdeControl$set_controller() opt in now, and the family can flip the default with plant re-baselined.

A finite-difference Jacobian’s small-component floor is 1e-5, not 1 and not the absolute tolerance. fd_jacobian() perturbs component j by rel_step * max(|y_j|, y_floor). With the floor at 1 a component of size 1e-5 was perturbed by a tenth of itself: on Robertson’s kinetics RODAS took 19000 steps where 436 do (435 with the analytic Jacobian; radau takes 144, being fifth order with Newton). With the floor at the absolute tolerance, the first fix tried, a tolerance of 1e-10 makes the perturbation 1e-16 and the difference vanishes in the subtraction (a unit derivative came back as 0, measured). The default is now rodas.f’s 1e-5 everywhere (jac_fd_floor). A component sitting on the upper edge of its domain, whose upward perturbation the right-hand side refuses, is perturbed downwards instead: the Jacobian does not depend on the step size, so retrying smaller could never have helped. ode_solve() also bounds the step by the span of times rather than the control’s plant-sized default of 10, which on the same problem was costing a thousand steps between t = 10 and t = 10⁴.

A minor bump: new capability to pin against, headers changed additively, no consumer’s floor moves (nothing in the family uses RODAS or the new stepper yet).

odelia 0.5.1

The spline reads as fast as 0.4.0’s again, with the same numbers. On a graded knot grid, 0.5.0 found a query’s span by binary search, where 0.4.0 guessed from the mean spacing and stepped from there (#21). plant’s adaptive light field is graded and read in height order, so the search made its FF16 runs 17% slower than on 0.4.0. hermite_spline now uses the guess-and-step lookup again, which returns exactly the span the search did: 2.55 million reads on random graded grids, knots and their neighbouring doubles included, are bit-identical. The front end’s unchecked operator() also skips the initialisation check, as 0.4.0’s did.

Against plant on 0.4.0, interleaved: FF16 full lifetime 0.106 to 0.102 s, K93 0.036 to 0.035 s, and TF24 6 years 2.35 to 2.31 s.

odelia 0.5.0

The spline’s backend becomes a cubic Hermite that can take a slope at each knot; the front end does not change, and neither do its numbers. A model that integrates over crowns knows the derivative at each knot as well as the value, and the global natural spline had no way to accept it. spline.hpp now holds hermite_spline<S>, a local cubic Hermite built from a value and a slope per knot. interpolator.hpp keeps every call the family makes — init(x, y), eval with its out-of-domain refusal, deriv, min/max, set_extrapolate, r_eval — on hermite_interpolator<S> (was basic_interpolator<S>), which derives from the backend so a caller that has slopes reaches init(x, y, m), value_and_slope, set_nodes/set_data and refine.

A caller with values alone gets the natural cubic spline’s own knot slopes (natural_slopes), and a Hermite read through those slopes is the natural spline on every span, so init(x, y) reads the curve it read before: test-drivers.R and test-spline.R are unchanged, and phylloptim 0.8.1 compiled against this release makes the same number of solver evaluations at the same speed (82.3 per solve, 3.0–3.1 µs, interleaved against 0.4.0), with outputs equal to rounding (worst 2e-10 absolute on its golden grid, from evaluating the same cubic in a different order). With supplied slopes the Hermite reproduces a cubic exactly and a read reaches four knots rather than all of them, which is what keeps an adjoint through it O(1). monotone_slopes (Fritsch–Carlson) is there for a caller that needs an interpolant that never leaves its data’s range.

value_with_slope<T> pairs a value with its slope so the two cannot be handed over separately and paired wrongly; util::to_passive strips every derivative layer from a scalar.

A minor bump: breaking only for a consumer that named basic_interpolator<S> or spline::basic_spline<S> directly. Split out of the reverse-mode change (#59).

odelia 0.4.0

Step rejection now works when the integration is pinned to fixed times (plant#642). advance_fixed() steps exactly to a caller-supplied set of times — how a replay reproduces a trajectory recorded earlier. It called the stepper bare, so both of the ways #55 gave a system to refuse a state were unreachable from it: a util::DomainError from a stage killed the solve, and a declared ode_state_valid() was never consulted at all.

The endpoints of a pinned step are the caller’s and cannot be moved, so #55’s answer — take a smaller step — becomes take several smaller ones to the same endpoint. On a refusal the sub-step shrinks through control.reject_step(), the same rule and the same step_size_min floor as the adaptive path, and the walk continues to the original time. That time is still hit exactly, so the times a caller records are unchanged. A system that raises no objection is stepped exactly as before: one RKCK step per interval, six stage evaluations, asserted in tests/standalone/.

An unreachable domain still fails, and now says where it gave up and why, rather than surfacing as whatever the offending stage happened to throw.

This was not a rare corner. plant’s mutant replay pins the stepper to a resident’s recorded times, and its TF24 model reports an empty carbon pool this way as a matter of routine — ~480 rejections in a resident run that goes on to complete normally — so a replay was near-certain to meet one and die. Invasion-fitness analysis was impossible for that model, not merely slow.

One caveat for systems that cache per-stage data through cache(system, rk_step): the stage indices restart at 0 on each sub-step, so a subdivided interval leaves the system holding the last sub-step’s stages rather than stages spanning the whole interval. A consumer recording such a cache for later replay gets a coarser record of a subdivided step than of a plain one.

A minor bump, so downstreams can pin against the capability (odelia (>= 0.4.0)). Systems that never refuse a state are unaffected.

odelia 0.3.1

Removes util::to_string_g(), added in 0.3.0 an hour earlier as a duplicate of util::format_double(). Both formatted a double for an error message with six significant figures — "%g" and "%.6g" are the same format string, since %g’s default precision is 6 — and they were byte-identical on every value tried. The one call site, in the invariant-rejection failure message, now uses format_double; its output is unchanged.

The duplicate arose because #55 was developed on a branch stacked below the 0.2.2 work that introduced format_double. Worth recording as a hazard rather than a review miss: two helpers with different names in different parts of the same file merge with no textual conflict, so neither the rebase nor the diff had anything to show.

A patch bump, although the header core lost a public symbol — the situation that earned 0.2.0 a minor one. The difference is what a version number can usefully say. 0.2.0’s removals had been reachable across released versions, so the bump warned of a break a consumer could actually hit. to_string_g never left this repository: 0.3.0 was never tagged — v0.2.1 remains the only tag and the only release — and both consumers are still on odelia (>= 0.2.x) with @v0.2.1 remotes. Nothing can have compiled against it, so a minor bump would announce an incompatibility that has no possible victim, and spend the signal for nothing.

Downstream floors are deliberately not raised. Neither plant nor phylloptim uses ode_state_valid() yet, and per the precedent set for leaf in 0.2.1 — checked rather than aligned for symmetry — raising a floor forces an upgrade for a change the consumer does not use. They should move to odelia (>= 0.3.0) when plant#609 or plant#599 actually adopts the domain check; 0.3.0 is the version that introduced it, and this release does not change it.

The 0.3.0 entry below has been corrected accordingly; it advertised a function that no longer exists.

odelia 0.3.0

Invariant-aware step rejection (#55). A system may now declare the domain its state lives in, and the adaptive stepper will refuse to commit a step that leaves it. Two ways to say so, both opt-in:

an optional bool ode_state_valid(const state_type&) const on the system, checked after each completed step;
util::stop_domain(msg), which throws the new util::DomainError, from anywhere in a stage.

Either one turns the step into a rejection — shrink and retry — rather than a committed out-of-domain state or, in the throwing case, a dead solve. Previously a throw from a stage ended the whole integration even though the pre-step state was still on the stack one frame up.

Only DomainError is caught. util::stop() and everything else still propagate, which is the point: absorbing them would turn a programming error into step-shrinking until “Cannot achieve the desired accuracy”, a diagnostic that points at the solver instead of at the bug.

This matters for bounded quantities that a finite RK step can overshoot even when the exact flow cannot — a carbon pool at zero, soil water at saturation, a probability at one. It is a discretisation guard, not a way to fix a model whose exact flow leaves its domain: that case shrinks to the minimum step and raises, now naming the reason and the location rather than blaming accuracy.

Version bumped so downstreams can pin against the capability (odelia (>= 0.3.0)); systems declaring neither hook are unaffected, and a Lorenz trajectory over 4127 steps is bit-identical across the change.

Numbers in that message are rendered with util::format_double(), so a step size at its floor reads 1e-08 rather than the "0.000000" std::to_string would give.

odelia 0.2.2

An out-of-domain interpolator lookup now says which point, how far out, and what the domain was. Interpolator::eval() had all three in hand — u, min() and max() — and reported none of them:

Extrapolation disabled and evaluation point outside of interpolated domain.

That sentence is the same whichever spline threw it, so a consumer holding several of them learns nothing about which one, and nothing about whether the point missed the near end or the far end. It cost real time downstream: localising traitecoevo/plant#576 meant instrumenting four call sites by hand to discover which spline was being asked and at what value, and the answer — the lower end, not past the far end as everyone had assumed — inverted the fix. Now:

Extrapolation disabled and evaluation point outside of interpolated domain:
u = -0.0023 lies 0.0023 beyond the lower end of [0, 6.8918].

Which spline, and which caller, is the one thing this layer cannot know; consumers that build several should catch and add it. A patch bump so downstreams can pin against the message.

Behaviour is otherwise unchanged. In particular the guard is still written u < min() || u > max() rather than the negation of an in-range test, because every comparison against NaN is false and a non-finite u must keep falling through to the spline — plant relies on that.

odelia 0.2.1

A patch bump for one reason: #46 has no version number, and a downstream needs one. d8235d1 (“Let a system compute its rates when the solver reads them”, #46) landed after 3bdfcf7 bumped the version to 0.2.0, and the version has not moved since — so 0.2.0 names two different header sets, one with #46 and one without, and no >= requirement can tell them apart.

That is not academic. traitecoevo/plant’s develop does not compile against the released 0.2.0:

odelia/ode_interface.hpp:212:3: error: 'this' argument to member function 'ode_rates'
  has type 'const plant::Patch<plant::FF16_Strategy, plant::FF16_Environment>',
  but function is not marked const

plant #585 made Patch::ode_rates non-const; r_ode_rates(const T& obj) and ode_solver_internal.hpp:155 both call it on a const&. #46 is the fix. Eight errors, four templated <Strategy, Environment> pairs × two call sites, and they surface inside these headers, which points nowhere near the cause. Confirmed pre-existing by syntax-checking plant’s origin/develop unmodified.

This is the same job 0.2.0 was bumped for, and the 0.2.0 entry below says so in as many words: it exists “to give downstream packages something to pin against … so a build against an older odelia fails at dependency resolution with a clear message rather than at compile time”. #46 needed the same courtesy and did not get it.

No header changes. Only DESCRIPTION. Downstream floors after this:

plant -> odelia (>= 0.2.1), because it links the ODE solver and needs #46.
leaf stays at odelia (>= 0.2.0). Checked rather than aligned for symmetry: leaf includes exactly one odelia header, odelia/interpolator.hpp, and never touches ode_rates or the solver. Raising its floor would force an upgrade for a fix in a header it does not include.

Closes #48.

odelia 0.2.0

A minor-version bump rather than a patch, because the header core lost public symbols: util::index, util::index_vector() and the base_1_to_0 / base_0_to_1 helpers are gone, and util::stop / util::warning no longer call into Rcpp. Nothing in the family used them — plant has its own plant::util equivalents — but a consumer that did would fail to compile, which is exactly what a version number is for. It also gives downstream packages something to pin against: leaf now requires odelia (>= 0.2.0), so a build against an older odelia fails at dependency resolution with a clear message rather than at compile time with RcppCommon.h: No such file or directory.

The header-only solver core is now free of R (#43). ode_util.hpp no longer includes RcppCommon.h, so everything reachable through interpolator.hpp, ode_control.hpp and ode_solver.hpp compiles and runs as plain C++ with no R installation — see tests/standalone/, which integrates the Lorenz system on a runner that has no R on it. util::stop() now throws std::runtime_error instead of calling Rcpp::stop(); Rcpp converts that into an ordinary R error with the same message at the package boundary, so R-level behaviour is unchanged apart from the condition’s class vector, which gains std::runtime_error in place of Rcpp::exception. util::warning() writes to std::cerr rather than raising an R warning; it had no callers. The unused util::index struct, its undefined Rcpp::as/wrap specializations, util::index_vector() and the base_1_to_0/base_0_to_1 helpers are removed — plant has its own. R remains where it belongs, in src/ and in solver_interface.hpp / rcpp_interface_helpers.hpp.

odelia 0.1.0

Odelia is a new package, arising out of https://github.com/traitecoevo/plant/. In that project, Rich FitzJohn built a custom ODE solver, using a Runge-Kutta 4-5 method, in C++. I’m spinning that code out into a package, as I want to use it elsewhere.

New implicit, adaptive-step RODAS4(3) Rosenbrock stepper for stiff systems, selectable via method = "rodas" when constructing a solver (#35). It reuses the existing adaptive step-size controller and obtains an exact Jacobian by forward-mode automatic differentiation; systems opt in by providing a template<class U> System<U> rebind() method. The explicit RKCK 4(5) method (method = "rkck") remains the default.
odelia now loads its shared library with global symbol visibility in .onLoad, so packages that LinkingTo: odelia and instantiate Solver can resolve the compiled XAD runtime symbols at load time without per-package linker hacks (#26).