Converting a d/dt() model puts nonlinear terms inside the rate constants,
where they cannot be iterated and where they violate an assumption the
event-sensitivity code already documents.
Updated while implementing this in #1213. Two things in the original
write-up turned out to be wrong and are corrected below: consuming the flag
column is not sufficient (it is per-equation, not per-term), and
rxIndLinStrategy()/rxIndLinState() cannot keep working the way this issue
assumed. A hand-written state-dependent k_from_to is also now a parse error
rather than being left alone -- a deliberate change of scope, not a
discovery.
Current behavior
indLin() (R/indLin.R:9-109) calls rxIndLin_() (:25), which returns a
character matrix carrying the state columns, a _rxF column and an indLin
flag column (built at src/expm.cpp:37-48). Lines :28-29 keep only the state
block and _rxF, and discard the flag column.
- So a Michaelis-Menten model converts to
k_central_output = vmax/(v*(km + central/v)) -- a state inside a rate
constant. Verified against current HEAD.
- That contradicts what the event-sensitivity jump code states at
inst/include/rxode2parseHandleEvid.h:1710-1719: "indLin's whole premise is a
CONSTANT (time-invariant) Jacobian over the interval -- the k_from_to rate
constants are parameter-functions only, never state-dependent".
- It is also why
rxSolve(<mm ode model>, method="indLin") is ~70% off against
liblsoda at the default hmax, measured on the 13-point sampling grid from
Sharif et al. 2022.
Change
R/indLin.R:9-109 -- consume the flag column. Corrected: the flag is
per equation, not per term -- .rxIndLinLine() returns one .fullIndLin
for the whole d/dt(state) line (R/rxIndLin.R:167) -- so consuming it cannot
produce the target shape below, which keeps ka*depot linear while routing only
the Michaelis-Menten residual to the forcing. Term-level routing is needed, and
the place to do it is where the per-term decision is already made and where the
flag is currently set: the three branches of .rxIndLinLine() that set
.fullIndLin <- TRUE (R/rxIndLin.R:95, :100, :115) route the whole term
to _rxF instead of emitting a state-dependent coefficient into A.
rxIndLinStrategy() / rxIndLinState() (R/rxIndLin.R:19,30) keep
controlling how a multi-state product is split. Corrected: they cannot.
Dividing any one state out of a multi-state product always leaves a coefficient
that still reads a state -- van der Pol's mu*y^2*dy leaves mu*y*dy or
mu*y^2 whichever is chosen -- so such a term goes to the forcing regardless of
the preference. Both stay exported and documented, with an @details note that
they no longer affect the split.
Target shape for the Michaelis-Menten example (unchanged, and achieved):
matExp()
cmt(depot)
cmt(central)
k_depot_central <- ka
cp <- central/v
indLin(central) <- -vmax*cp/(km + cp)
Correct value after this issue
rxSolve(<ode model>, method="indLin") is accurate at the default hmax.
Hand-written state-dependent k_from_to models are untouched and keep
working. Superseded: a matExp() rate constant that depends on a
compartment is now a parse error, naming the constant, the compartment it
reaches and the indLin() form to use instead. The matrix exponential is only
valid when the rate matrix is constant over the step, so such a model was
silently wrong rather than merely unsupported. Models generated by
rxSensMatExp() are exempt for now: their sensitivity blocks are built out of
rate constants throughout and reference the primal ones by name, so rewriting
that generator is separate work and is not done here.
Test
indLin(mmOde) emits indLin(central) <- ... with no k_central_output at
all (the elimination is entirely in the forcing)
- solving matches
liblsoda at the default settings -- measured relative error
3.9e-06 against the ~70% above
the van der Pol block that exercises rxIndLinStrategy()/rxIndLinState()
still runs unchanged Corrected: that block was rewritten rather than
re-enabled. Its method="lsoda" arm was never lsoda -- rxSolve.default
force-selects method="indLin" for any model carrying an indLin descriptor,
so it was comparing the model against itself and the assertion was vacuous. It
now compares against a plain ODE model, which is what the commented-out
assertions were reaching for, and additionally asserts that the three strategy
variants produce identical models.
Resolution
Implemented in #1213 (branch fix/indlin-iteration), together with #1185.
Converting a
d/dt()model puts nonlinear terms inside the rate constants,where they cannot be iterated and where they violate an assumption the
event-sensitivity code already documents.
Current behavior
indLin()(R/indLin.R:9-109) callsrxIndLin_()(:25), which returns acharacter matrix carrying the state columns, a
_rxFcolumn and anindLinflag column (built at
src/expm.cpp:37-48). Lines:28-29keep only the stateblock and
_rxF, and discard the flag column.k_central_output = vmax/(v*(km + central/v))-- a state inside a rateconstant. Verified against current HEAD.
inst/include/rxode2parseHandleEvid.h:1710-1719: "indLin's whole premise is aCONSTANT (time-invariant) Jacobian over the interval -- the k_from_to rate
constants are parameter-functions only, never state-dependent".
rxSolve(<mm ode model>, method="indLin")is ~70% off againstliblsodaat the defaulthmax, measured on the 13-point sampling grid fromSharif et al. 2022.
Change
Corrected: the flag isR/indLin.R:9-109-- consume the flag column.per equation, not per term --
.rxIndLinLine()returns one.fullIndLinfor the whole
d/dt(state)line (R/rxIndLin.R:167) -- so consuming it cannotproduce the target shape below, which keeps
ka*depotlinear while routing onlythe Michaelis-Menten residual to the forcing. Term-level routing is needed, and
the place to do it is where the per-term decision is already made and where the
flag is currently set: the three branches of
.rxIndLinLine()that set.fullIndLin <- TRUE(R/rxIndLin.R:95,:100,:115) route the whole termto
_rxFinstead of emitting a state-dependent coefficient intoA.Corrected: they cannot.rxIndLinStrategy()/rxIndLinState()(R/rxIndLin.R:19,30) keepcontrolling how a multi-state product is split.
Dividing any one state out of a multi-state product always leaves a coefficient
that still reads a state -- van der Pol's
mu*y^2*dyleavesmu*y*dyormu*y^2whichever is chosen -- so such a term goes to the forcing regardless ofthe preference. Both stay exported and documented, with an
@detailsnote thatthey no longer affect the split.
Target shape for the Michaelis-Menten example (unchanged, and achieved):
Correct value after this issue
rxSolve(<ode model>, method="indLin")is accurate at the defaulthmax.Hand-written state-dependentSuperseded: ak_from_tomodels are untouched and keepworking.
matExp()rate constant that depends on acompartment is now a parse error, naming the constant, the compartment it
reaches and the
indLin()form to use instead. The matrix exponential is onlyvalid when the rate matrix is constant over the step, so such a model was
silently wrong rather than merely unsupported. Models generated by
rxSensMatExp()are exempt for now: their sensitivity blocks are built out ofrate constants throughout and reference the primal ones by name, so rewriting
that generator is separate work and is not done here.
Test
indLin(mmOde)emitsindLin(central) <- ...with nok_central_outputatall (the elimination is entirely in the forcing)
liblsodaat the default settings -- measured relative error3.9e-06 against the ~70% above
the van der Pol block that exercisesCorrected: that block was rewritten rather thanrxIndLinStrategy()/rxIndLinState()still runs unchanged
re-enabled. Its
method="lsoda"arm was never lsoda --rxSolve.defaultforce-selects
method="indLin"for any model carrying an indLin descriptor,so it was comparing the model against itself and the assertion was vacuous. It
now compares against a plain ODE model, which is what the commented-out
assertions were reaching for, and additionally asserts that the three strategy
variants produce identical models.
Resolution
Implemented in #1213 (branch
fix/indlin-iteration), together with #1185.