Skip to content

Forward-sensitivity integration exhausts mxstep at a fitted if(t>=...) switch time, killing gradient fits of models with moving discontinuities #82

Description

@wshlavacek

On a model whose rate laws gate on if(t>=p, ...) with p a fitted parameter, forward-sensitivity integration exhausts its step budget at the switch time and the fit dies:

[ERROR][CVode] At t = 49.4914, mxstep steps taken before reaching tout.

t = 49.49 is where one of the model's switch times falls for the candidate being evaluated. The plain state integration survives the same model; the sensitivity integration does not.

Model

Lin-2021/nyc_multiphase (PyBNF jobs collection). Its functions block gates on several fitted times:

U()      = if(t>=t0, 1, 0)
Lambda() = if(t>=sigma && t<tau1, lambda0, if(t>=tau1 && t<200, lambda1, 0))
kQfunc() = if(t>=t0, kQ, 0)      jQfunc(), cIfunc(), phiM() -- same

t0, sigma, tau1 are all fit variables, so the discontinuity moves with every candidate and cannot be registered as a fixed stop time in advance.

Likely mechanism, and how it relates to #55

These are Functional rate laws, so df/dp goes down the CVODES difference-quotient path that #55 exists to remove. A difference quotient perturbs p — which moves the switch time — so the quotient straddles the jump and reports an enormous, essentially meaningless derivative. The integrator then chases it with ever-smaller steps until the budget is gone.

If that reading is right, #55 (and its staged children #66/#67/#68) is the fix rather than a step-count tweak: an analytic df/dp never perturbs the switch time to begin with. #49 is adjacent but distinct — it concerns the correctness of dt*/dp for a switch time, not the integration blowing its budget.

The same pathology is visible without sensitivities, more mildly

Running the plain state solve on this model, CVODES reports step-size collapse and continues:

[CVode] Internal t = 69.8066 and h = 3.5245e-15 are such that t + h = t on the next step.
        The solver will continue anyway.
[CVode] The above warning has been issued mxhnil times and will not be issued again.

Two things worth noting there. The mxhnil cap means the warning count systematically understates how much of this is happening. And "will continue anyway" means the caller gets no signal at all — in a 5,936-simulation fit, 89 solves started and never logged completion, and step counts ran to a median of 1,085 and a max of 1,962 for a small ODE over 170 time units.

This is the same class as #54 (closed), which bounded the unbounded case; it did not make the crossing cheap, and the sensitivity path evidently still exhausts the budget.

Impact

This is the one model of six in a PyBNF fitting benchmark where the in-process bngsim backend is slower than the legacy BNG2.pl subprocess stack (26.6 min against 11.3 min for identical work). BNG2.pl integrates to each sample time — 170 short solves that restart every 1.0 time units — which bounds how far a mid-interval discontinuity can poison the integration; a single continuous t=[0,170] call has no such natural barrier. That may be worth considering independently of the sensitivity failure.

Environment

bngsim 0.11.35 (v0.11.35-77-g32395e4), PyBNF 1.6.0, Python 3.12, macOS 15.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions