Integrator state not updated after setValue when using oneStep with reset=False

[matthiaskoenig]

Hi all,

I did some additional digging into the following issue:
#1232
This turns out to be a major problem for us.

We can partially work around it by calling oneStep with reset=False, which continues the integration while preserving the integrator state. This flag is exactly what we need. However, there is a bug: values set via setValue are not propagated into the integrator state between oneStep calls.

Example:

import roadrunner
r: roadrunner.RoadRunner = roadrunner.RoadRunner("spt_liver.xml")
r.setValue("[S_ext]", 1.0)
print({key: r.getValue(key) for key in r.selections})

print("--- [S_ext] change lost ---")
_new_time = r.oneStep(currentTime=0, stepSize=10, reset=False)
print({key: r.getValue(key) for key in r.selections})

This produces:

/home/mkoenig/git/pkdb_models/.venv/bin/python /home/mkoenig/git/pkdb_models/pkdb_models/models/tpm_spt/models/results/models/bug_onestep_reset.py 
{'time': 0.0, '[S_ext]': 1.0, '[P_ext]': 0.0, '[S]': 0.0, '[P]': 0.0, '[T]': 0.0, '[protein]': 0.3333333333333333, 'tsim': 0.0}
--- [S_ext] change lost ---
{'time': 10.0, '[S_ext]': 0.0, '[P_ext]': 0.0, '[S]': 0.0, '[P]': 0.0, '[T]': 0.0, '[protein]': 0.3333333333333333, 'tsim': 10.0}

As shown, the model modification to [S_ext] is lost during integration.
Values changed via setValue are simply not pushed into the integrator state.

A minimal example project is attached here:
[spt_liver.zip](https://github.com/user-attachments/files/23794143/spt_liver.zip)

Source of the issue

The problematic behavior is visible in the implementation of oneStep.
When reset == False, the integration continues without ensuring that any updated model values (set by setValue) are transferred to the integrator:

double RoadRunner::oneStep(const double currentTime, const double stepSize, const bool reset) {
    get_self();
    check_model();
    applySimulateOptions();

    try {
        if (reset) {
            self.integrator->restart(currentTime);
        }
        return self.integrator->integrate(currentTime, stepSize);
    }
    catch (EventListenerException &e) {
        rrLog(Logger::LOG_NOTICE) << e.what();
        return self.model->getTime();
    }
}

A missing branch should ensure that changed model values are synchronized when reset == false, e.g.:

if (!reset) {
    self.integrator->updateChangedValuesFromModel;
}

Alternatively, setValue should directly update the integrator state.

Documentation issue

The reset argument in oneStep is also poorly documented and misleadingly named.
The option does not "reset" the model; instead, it controls whether the integrator is restarted.

The Python docstring currently shows only two arguments and does not describe reset at all:

def oneStep(self, currentTime: "double", stepSize: "double", reset: "bool"=True) -> "double":
    r"""

    RoadRunner.oneStep(startTime, stepSize)

    Carry out a one step integration of the model. The method takes two arguments,
    the current time and the step size to us in the integration. The method returns
    the new time which will be currentTime + StepSize::

        newTime = rr.oneStep (10, 0.5)

    """

This should be updated to correctly document all three parameters and clarify the meaning of the reset/restart flag.

[matthiaskoenig]

Probably requires something along the line of the following for the CVODE integrator to update the states before continuing the integration (see restart function). I.e. if the states were updated due to setValues the updated state must be copied in the CVODE memory

// copy state std::vector into cvode memory
if (mStateVector && mCVODE_Memory) {
    mModel->getStateVector(NV_DATA_S(mStateVector));
}

[matthiaskoenig]

@luciansmith Could you have a look at the issue? I have the feeling this should be pretty simple to fix (see suggestions I made). As I understand the issue it only requires to synchronize the changes between the model object and the integrator.

[luciansmith]

OK, I did indeed look into this, and unfortunately, I think what you are asking for might be impossible.

In order to get the current state vector into CVODE, we need to call two functions: the first to update the state vector, and the second to copy that vector into CVODE. The first you do with

mModel->getStateVector(NV_DATA_S(mStateVector));

and the second is done with:

CVodeReInit(mCVODE_Memory, time, mStateVector);

I don't see any other way to copy a state vector into CVODE without calling 'init' or 'reinit'.

This makes sense conceptually to me: if you look at the documentation for the reinit function, it says that the function is used when you are asking CVODE to solve a new problem with the same shape as the old problem. If you've changed a value, the internal state of CVODE (which presumably involves slopes and second derivatives) needs to also change, so it's a new problem with the same shape as the old problem, and you really can't 'just continue'.

To make sure that you transition smoothly when you don't have to change anything, I think I would recommend doing something like the following:

for x in range(steps):
    reset = False
    if (needToChangeSomething):
        r.S1 = 5
        reset = True
    r.oneStep(t, hstep, reset)
    t = t + hstep

[luciansmith]

One other thing to try: there's another function in roadrunner called 'internalOneStep', which behaves exactly the same as 'oneStep' except that it sets the integrator's 'variable_step_size' setting to 'True'. So you might see if turning on variable step size helps with your original problem (#1232)

[matthiaskoenig]

Here’s an improved version of your message that maintains the content while enhancing clarity and tone:
Hi Lucian,

Thanks for looking into this!

Unfortunately, the models do not integrate correctly with multiple oneStep calls. As discussed in [this issue](#1232), the model can be integrated using a single simulate call, but it does not work properly when using repeated oneStep calls. This behavior makes sense to me because of the internal solver state that tracks whether the step size was reduced in the previous time steps. CVODE's internal settings switch accordingly, which means that integrating via simulate (or a single oneStep call) differs significantly from multiple oneStep calls where the solver is reinitialized at each step.

It would be helpful to have a brief call where I can show you some examples and explain in more detail what’s going on.

Regarding the variable step size, I don't think it will help in this case. It essentially runs the same integration process, but without interpolation at the end. Moreover, it complicates matters, as I would lose control over the timestep, which is critical for synchronizing with other simulation frameworks.

Please see the attached example, where:

• The model can be successfully simulated using simulate.
• The model can be successfully simulated with oneStep when reset=False.
• However, if you change reset=True in line 94, the integration fails.

_new_time = r.oneStep(t_model, dt, reset=False)

In summary, the model fails to integrate when using multiple oneStep calls, especially with reset=True. We've observed this issue with several of our models, and it strongly depends on the model's trajectory. The problem arises because we lose the internal CVODE state during the reinitialization.

It would be great if we could have a short chat to discuss this further.

Best regards,
Matthias

[luciansmith]

I am amused by the "Here’s an improved version of your message" bit ;-) Although at first I thought you were improving my message!

I think I do understand the issue, though I'm also more than happy to have a call about it.

The fundamental problem is (I think!) that if you change an element's value, you can't keep 'the internal CVODE state' any more, because that internal state is fundamentally dependent on the particular values of every element. Change one, and you'd have to change the internal state as well. I could not find a function in the CVODE interface that does this, at any rate: the only way I found to update the state vector inside CVODE is to use the CVodeReInit function, which presumably ditches the internal state and rebuilds a new one (and if you call it all the time, causes some of your simulations to fail).

So as a compromise, I would change your line 94 (which I couldn't find? neither spt_liver.zip nor droplets_iri_liver.zip had it) to something like:

_new_time = r.oneStep(t_model, dt, reset=needToReset)

and have needToReset be False by default, but change it to True when you need to change a model element's value.

(I suppose another option that might work in some cases is to put the logic that changes particular values in as events in your SBML model directly. That's not always going to be possible, of course.)

[matthiaskoenig]

Hi Lucian,

I occasionally clean up my English to make things clearer — I'm not a native speaker, so this helps ;).
Here are my thoughts:

• Yes, one has to use CVodeReInit. I looked into CVODE and what is exposed, and this seems to be the only safe/good way to do it.
• I perform updates after every onestep because this is a coupled simulation. After each timestep, another simulation framework is updated (i.e., it’s part of a hybrid simulation).
• The system can be integrated — the solver is able to handle the full time course; otherwise, simulate would fail. So it must be something about the integration settings that changes in later onestep calls compared to the initial one.
• These settings must be state-dependent. The main suspects that come to mind are tolerances:
  • Tolerance settings are likely computed based on the initial state. I assume there is an absolute tolerance vector that is recalculated on each ReInit based on the current state (and thus differs at every onestep). Ideally, I would compute this vector once at t = 0 and reuse it for all subsequent ReInit calls.
  • Possibly also step-size-related settings (max step size, min step size, initial step size), which might be derived from the state on initialization.

So my question is: which state-dependent integrator settings does RoadRunner/CVODE apply (e.g., absolute tolerance vectors), and how can I reuse the values computed at the initial onestep?

Best,
Matthias

[luciansmith]

Oh, hmm. The absolute tolerances do change based on the initial state, and there is indeed a line in oneStep that resets the tolerances.

You can get around this by setting the absolute tolerances to a vector instead of a single value:

abstol_vec = r.integrator.getAbsoluteToleranceVector()
r.integrator.setValue("absolute_tolerance", list(abstol_vec))

(Having to cast it to a list is annoying; I'll add that to the 'to fix' list...)

Once it's a vector and not a single value, roadrunner will use the vector as-is, and not try to scale it any more.

Unfortunately, I tried this in the iri_droplets_liver example you provided, and it made a difference, but not a helpful one: instead of making it to time point 180 before failing, it made it to time point 174 before failing. But the absolute tolerance vector was the same the whole time!

I'll experiment with other things to see what else we do that might be the culprit.

[matthiaskoenig]

Thanks. I will experiment with this. The other setting which is definitely changed is the initial time step which is estimated base on the state vector. If I could read this out on the initial step I could set this on the later steps instead of providing 0.0 which results in estimation. See below:

[matthiaskoenig]

The setting of the initial absolute tolerance vector on all subsequent one-steps solved the issue ;).
Perfect.

So if this could be made easier/more efficient to set, i.e. addressing #1301 would be great.

It would be amazing to have a flag for onestep to keep the initial tolerances so these don't have to be set in every step. We basically run 10000 of micro-simulations in coupled FEM simulations and would now require to reset for every roadrunner instance on every onestep (also a few thousands per simulator) the absolute tolerance vector. I.e. we have to pass the tolerance vector a few million times per simulations around ;)
Everything which we can avoid passing around especially in the context of MPI on clusters helps.

[luciansmith]

Once the tolerance is a vector, you shouldn't have to keep setting it; once is enough.

Internally, though, it would certainly be more efficient to not reset it on every call to oneStep. I'll have to check if we actually need it at all; it may be already set.

I'm surprised it worked for you, though pleased! I tried the iri_droplets_liver python program you uploaded to the previous issue, and as I mentioned, it changed the behavior, but didn't solve the overall problem. (Though maybe resetting it every time helped for some unknown reason?) Here's what I changed:

# FEBio version (onestep)
print("--- onestep --- ")
r.resetAll()
abstol_vec = r.integrator.getAbsoluteToleranceVector()
r.integrator.setValue("absolute_tolerance", list(abstol_vec))
for key, value in changes.items():
    r.setValue(key, value)

t_model = 0.0
while t_model < tend:
    s = r.internalOneStep(t_model, dt)
    print(s)
    t_model += dt
#df2 = pd.DataFrame(s, columns=s1.colnames) #<-Doesn't work with onestep!
#print(df2.head)

[matthiaskoenig]

I do use the OneStep not the internalOneStep.
Attached the full code which works. Also made 2 other models work with the absolute tolerances so the reuse of the absolute tolerances fixes the issue.

"""Example script for porous media coupling for IRI & Droplet model."""
from pathlib import Path

import pandas as pd
import roadrunner

# -------------------------------------------------------------------------------------
# load model
# -------------------------------------------------------------------------------------
model_path = Path(__file__).parent /  "droplets_iri_liver.xml"
r = roadrunner.RoadRunner(str(model_path))

# -------------------------------------------------------------------------------------
# set changes for given simulator
# -------------------------------------------------------------------------------------
vol_point = 1.0  # [l]
f_fluid = 0.3  # [dimensionless] fluid fraction
f_solid = 1 - f_fluid  # [dimensionless] solid fraction

changes = {
    # setting volumes
    "Vext": vol_point * f_fluid,  # [liter] fluid phase volume
    "Vli": vol_point * f_solid,  # [liter] solid phase volume

    # concentrations of external concentrations in fluid phase
    "[glc_ext]": 5.0,  # [mM]
    "[lac_ext]": 1.0,  # [mM]
    "[o2_ext]": 1.2,  # [mM]
    "[alt_ext]": 0.0,  # [mM]
    "[ast_ext]": 0.0,  # [mM]
    "[ffa_ext]": 0.0,  # [mM]
    "[vldl_ext]": 0.0,  # [mM]

    "[glc]": 5.0,  # [mM]
    "[lac]": 1.0,  # [mM]
    "[o2]": 1.2,  # [mM]
    "[alt]": 1.0,  # [mM]
    "[ast]": 0.7,  # [mM]
    "[pyr]": 0.8,  # [mM]
    "[atp]": 3.0,  # [mM]
    "[adp]": 1.0,  # [mM]
    "[nadh]": 3.9,  # [mM]
    "[nad]": 0.1,  # [mM]
    "[ros]": 0.0,  # [mM]

    # set temperature
    "temperature": 37 + 273.15,  # [Kelvin] warm simulation
    # "temperature": 4 + 273.15,  # [Kelvin] cold simulation
}

# -------------------------------------------------------------------------------------
# Timecourse
# -------------------------------------------------------------------------------------
# run simulation
tend = 600  # [min]
steps = 10
dt = tend/steps  # [min]

# single simulation (single simulation works)
print("--- simulate --- ")

# apply changes
r.resetAll()
for key, value in changes.items():
    r.setValue(key, value)

# simulate
s1 = r.simulate(start=0.0, end=tend, steps=steps) #*60
df1 = pd.DataFrame(s1, columns=s1.colnames)
print(df1.head())

# FEBio version (onestep/simulate)
print("--- onestep --- ")
# apply chances
r.resetAll()
for key, value in changes.items():
    r.setValue(key, value)

# initial tolerances
integrator: roadrunner.Integrator = r.integrator
init_tolerances = integrator.getAbsoluteToleranceVector()

t_model = 0.0
all_data = [
    {key: r.getValue(key) for key in r.selections}
]
while t_model < tend:
    r.integrator.setValue("absolute_tolerance", list(init_tolerances))

    _new_time = r.oneStep(t_model, dt, reset=True)
    all_data.append(
        {key: r.getValue(key) for key in r.selections}
    )

    t_model += dt

df_onestep = pd.DataFrame(all_data)
print(df_onestep.head())

Results in identical results:

--- simulate --- 
    time     [ins]     [ffa]  [vldl_ext]      [dg]  [protein_dgat2]  [protein_ldsyn]       [tg1]       [tg2]  ...     [adp]    [nadh]         [nad]     [ros]         [alt]         [ast]  [necrosis]       [tgd]   tsim
0    0.0  0.000001  2.142857    0.000000  1.714286         1.428571         1.428571    0.000000    0.000000  ...  1.000000  3.900000  1.000000e-01  0.000000  1.000000e+00  7.000000e-01    0.000000    0.000000    0.0
1   60.0  0.000001  0.000449    0.032462  1.733773         1.428571         1.428571  186.116700  184.993032  ...  1.167704  3.993897  6.103139e-03  0.000000  1.000000e+00  7.000000e-01    0.000000  496.160247   60.0
2  120.0  0.000001  0.000449    0.045561  1.741488         1.428571         1.428571   68.630327  122.233641  ...  2.128399  4.000000  5.262470e-12  0.098465  1.000000e+00  7.000000e-01    0.000000  498.800256  120.0
3  180.0  0.000001  0.000449    0.049510  1.743797         1.428571         1.428571   55.132707  101.127744  ...  4.000000  4.000000  3.443320e-15  0.098716  5.130573e-08  3.591401e-08    1.428571  499.047986  180.0
4  240.0  0.000001  0.000449    0.049510  1.743797         1.428571         1.428571   55.132707  101.127744  ...  4.000000  4.000000  3.443320e-15  0.098716  2.983829e-13  2.088680e-13    1.428571  499.047986  240.0

[5 rows x 30 columns]
--- onestep --- 
    time     [ins]     [ffa]  [vldl_ext]      [dg]  [protein_dgat2]  [protein_ldsyn]       [tg1]       [tg2]  ...     [adp]    [nadh]         [nad]     [ros]         [alt]         [ast]  [necrosis]       [tgd]   tsim
0    0.0  0.000001  2.142857    0.000000  1.714286         1.428571         1.428571    0.000000    0.000000  ...  1.000000  3.900000  1.000000e-01  0.000000  1.000000e+00  7.000000e-01    0.000000    0.000000    0.0
1   60.0  0.000001  0.000449    0.032462  1.733773         1.428571         1.428571  186.116700  184.993032  ...  1.167704  3.993897  6.103139e-03  0.000000  1.000000e+00  7.000000e-01    0.000000  496.160247   60.0
2  120.0  0.000001  0.000449    0.045561  1.741488         1.428571         1.428571   68.630335  122.233645  ...  2.128399  4.000000  5.264234e-12  0.098465  1.000000e+00  7.000000e-01    0.000000  498.800256  120.0
3  180.0  0.000001  0.000449    0.049510  1.743797         1.428571         1.428571   55.132750  101.127796  ...  4.000000  4.000000  3.442532e-15  0.098716  5.130657e-08  3.591460e-08    1.428571  499.047986  180.0
4  240.0  0.000001  0.000449    0.049510  1.743797         1.428571         1.428571   55.132750  101.127796  ...  4.000000  4.000000  3.442532e-15  0.098716 -2.045335e-12 -1.431734e-12    1.428571  499.047986  240.0

[luciansmith]

OK, I'm glad you've found something that works! A couple observations:

• If you change the number of steps from 10 to 100, you get the behavior I saw, where oneStep doesn't work even with the absolute tolerance vector. I'm glad it works sometimes, but I think you might be right about the initial time step issue, so I'll keep looking.
• You don't need to set the absolute tolerance vector every time. Just set it once before the oneStep loop, and it'll still work.

[matthiaskoenig]

I think the core issue is that the system is poorly conditioned: the state variables span many orders of magnitude, from around 1×10⁻⁹ up to 10³. The SBML model is formally valid, but integrators struggle when variables differ so drastically in scale—especially in larger models.

Two main approaches come to mind:

1. Rewrite or rescale the equations

In principle, one could rewrite the SBML model so that all state variables lie in similar numerical ranges. In practice, this quickly becomes complicated and error-prone. For example, drug concentrations may be in the pM range (10⁻⁹ mM) while glucose is around 5 mM—a difference of 10⁹ in scale. Handling species in different units requires conversion factors, which some simulators do not support correctly, and forces parameters such as Km to be expressed in different units (e.g., pM for drugs, mM for glucose).
This significantly complicates the modeling process and does not guarantee robustness, especially when combining multiple models or using hierarchical model composition. We have experimented with such rescaling, but it is not sustainable.
In theory, one could imagine an SBML-to-SBML converter that automatically rescales variables to produce a numerically better-conditioned model.

2. Apply preconditioning / state-variable scaling in CVODE

A much more practical solution is to use preconditioning and appropriate state-variable scaling directly within the integrator. CVODE is designed for this: preconditioning can greatly improve robustness and should ideally relieve modelers from having to think about numerical tolerances, step sizes, or internal solver details when writing SBML models—which is, after all, the intent of SBML.

When preconditioning is used, user-supplied functions define the linear operators corresponding to the left and right preconditioner matrices (P_1) and (P_2).
My systems are exactly the kind of models that benefit from preconditioning. You should use a preconditioner when:

• Corrector convergence fails (→ yes)
• Krylov iterations regularly exceed 5–10
• The Jacobian is stiff or ill-conditioned (→ yes)
• The system size is > ~50 ODEs (→ yes)
• There are slow convergences after events or discontinuities
• The system has strong multiscale dynamics (PBPK, hepatic models, metabolic networks) (→ yes)

Without preconditioning:

• Newton iterations fail
• The step size collapses to hmin (→ yes)
• CVODE reports: “corrector convergence test failed repeatedly or with |h| = hmin” (→ yes)

(See figure below.)

Enabling preconditioners in CVODE would allow much more robust simulation of large, stiff, multi-scale systems with widely differing state-variable scales—for example, by treating subsystems with domain-specific preconditioners.