dynamicslab
[BUG] Catch and re-raise simulation errors when ODE blows up
When the discovered equation blows up numerically, scipy.integrate.solve_ivp emits a ValueError. I would like to catch this error, but ValueError is also raised in SINDy.simulate() because it's a pretty general exception. I don't want to catch those ValueErrors
I believe SINDy.simulate() should catch the ValueError and re-raise it as another class, e.g.
SimulationError = type("SimulationError", (UserWarning,), {})
Also - usually this failure is preceded by a long period of the integrator struggling through LSODA warnings (which are just UserWarning). UserWarning is pretty generic too; should we a add an arg to simulate() that raises the LSODA warnings to Exceptions?
Reproducing code example:
import warnings
import pysindy as ps
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
seed = 20
rng = np.random.default_rng(seed)
t = np.linspace(0,1,100)
x_true = np.stack((np.sin(t), np.cos(t)), axis=-1)
x_train = rng.normal(size=x_true.shape)
lib = ps.PolynomialLibrary(degree=3)
opt = ps.STLSQ(threshold=.001, alpha=0)
model = ps.SINDy(feature_library=lib, optimizer=opt)
model.fit(x=x_train, t=t)
x_sim = model.simulate((.5, .5), t)
plt.plot(x_sim)
Error message:
File [~/github/gen-experiments/env/lib/python3.10/site-packages/pysindy/pysindy.py:898](https://file+.vscode-resource.vscode-cdn.net/home/xenophon/github/gen-experiments/~/github/gen-experiments/env/lib/python3.10/site-packages/pysindy/pysindy.py:898), in SINDy.simulate(self, x0, t, u, integrator, stop_condition, interpolator, integrator_kws, interpolator_kws)
[894](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/pysindy/pysindy.py?line=893) # Need to hard-code below, because odeint and solve_ivp
[895](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/pysindy/pysindy.py?line=894) # have different syntax and integration options.
[896](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/pysindy/pysindy.py?line=895) if integrator == "solve_ivp":
[897](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/pysindy/pysindy.py?line=896) return (
--> [898](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/pysindy/pysindy.py?line=897) (solve_ivp(rhs, (t[0], t[-1]), x0, t_eval=t, **integrator_kws)).y
[899](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/pysindy/pysindy.py?line=898) ).T
[900](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/pysindy/pysindy.py?line=899) elif integrator == "odeint":
[901](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/pysindy/pysindy.py?line=900) if integrator_kws.get("method") == "LSODA":
File [~/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/ivp.py:591](https://file+.vscode-resource.vscode-cdn.net/home/xenophon/github/gen-experiments/~/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/ivp.py:591), in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options)
[589](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/ivp.py?line=588) status = None
[590](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/ivp.py?line=589) while status is None:
--> [591](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/ivp.py?line=590) message = solver.step()
[593](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/ivp.py?line=592) if solver.status == 'finished':
[594](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/ivp.py?line=593) status = 0
File [~/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/base.py:181](https://file+.vscode-resource.vscode-cdn.net/home/xenophon/github/gen-experiments/~/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/base.py:181), in OdeSolver.step(self)
[179](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/base.py?line=178) else:
[180](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/base.py?line=179) t = self.t
--> [181](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/base.py?line=180) success, message = self._step_impl()
[183](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/base.py?line=182) if not success:
[184](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/base.py?line=183) self.status = 'failed'
File [~/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/lsoda.py:152](https://file+.vscode-resource.vscode-cdn.net/home/xenophon/github/gen-experiments/~/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/lsoda.py:152), in LSODA._step_impl(self)
[150](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/lsoda.py?line=149) itask = integrator.call_args[2]
[151](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/lsoda.py?line=150) integrator.call_args[2] = 5
--> [152](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/lsoda.py?line=151) solver._y, solver.t = integrator.run(
[153](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/lsoda.py?line=152) solver.f, solver.jac or (lambda: None), solver._y, solver.t,
[154](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/lsoda.py?line=153) self.t_bound, solver.f_params, solver.jac_params)
[155](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/lsoda.py?line=154) integrator.call_args[2] = itask
[157](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/scipy/integrate/_ivp/lsoda.py?line=156) if solver.successful():
...
[159](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/sklearn/utils/validation.py?line=158) "#estimators-that-handle-nan-values"
[160](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/sklearn/utils/validation.py?line=159) )
--> [161](file:///home/xenophon/github/gen-experiments/env/lib/python3.10/site-packages/sklearn/utils/validation.py?line=160) raise ValueError(msg_err)
ValueError: Input contains infinity or a value too large for dtype('float64').
Context:
1.7.6.dev55+g5f224ab https://github.com/dynamicslab/pysindy/blob/5f224ab35508dbd0e8ed564a3a5d522b6dc5f19a/pysindy/pysindy.py#L896
I think it's a good idea to catch the ValueError in simulate and raise a less generic error. As for the LSODA warnings, I found in past experience that it's actually quite difficult to suppress or catch all of them. This is because they are actually printed to the system stderr from a FORTRAN routine that is invoked rather than being raised as warnings that python can actually notice. You can do some elaborate hacks to get around it: https://stackoverflow.com/questions/31681946/disable-warnings-originating-from-scipy. But I would advise just leaving the LSODA warnings as they are -- its only worth the effort if you really need to get rid of them.
Yes Zach - I also found the scipy issue on that topic. I don't want to silence them - I want to add a warning filter to elevate them to an error, so that that trial quits. Generally LSODA warnings mean that the trial failed, but that the integrator will churn a long time. I just want to catch that exception and move to the next trial.
Yes, its a big headache to catch the warnings...my understanding is that you'll have to get a good understanding of the stdout and stderr streams that exist in the process. I think that jupyter and python create some additional stdout and stderr streams for the notebook, and you can use dup2 to redirect these streams. (All this may be architecture/bash dependent) You can probably redirect the system stderr stream and raise a python error/warning if LSODA prints to it somehow, but yeah not easy
ps: this is really all bash nonsense that python wraps around. I suggest trying to understand bash from resources like this: https://man7.org/linux/man-pages/man2/dup.2.html
Dealing with the warning message isn't really where I'm going. My point is that even though the fortran code prints the warning, scipy still detects it and raises a UserWarning quietly. A warnings filter that ignores the warning can't prevent the fortran print-out, but a filter that elevates the UserWarning to an Exception will kill the integration (after fortran prints the warning once).
If the exception is raised, my experimental code can catch it and then move to the next trial.
Ah, gotcha, I didn't know that scipy raised the UserWarning. Well just note that the LSODA warnings are not always an indication of instability or a problem, so sometimes they don't mean that the integration has failed...it is very problem dependent whether they will appear or not. I've seen times where it is just a stiff event that needs care to integrate past and a warning pops up.
One extra note: sometimes I want to show visually that the integrator blew up. solve_ivp doesn't output things if it fails/blows up, but odeint does. This is really nice for e.g. showing that the trapping SINDy produces stable trajectories but regular SINDy does not. So I often use odeint to compare these optimizers and show divergent integration.
LSODA warnings are not always an indication of instability or a problem, so sometimes they don't mean that the integration has failed...it is very problem dependent whether they will appear or not. I've seen times where it is just a stiff event that needs care to integrate past and a warning pops up.
@znicolaou that's a good point! It would be impossible (at least in single-threaded python) to catching and convert the warning, but continuing the integration unless user elevates it to an exception. Let me step back from that suggestion then.
sometimes I want to show visually that the integrator blew up. solve_ivp doesn't output things if it fails/blows up, but odeint does.
@akaptano I've seen solve_ivp raise a ValueError when it blows up enough to get to np.inf/np.nan, but nothing before. What does odeint do?
odeint just returns the trajectory anyways! So you can plot the finite values and see the blow up