capi_return is NULL result from model.simulate

[itsjacobhere]

After fitting a sindy model I'm trying to simulate the model forward on a set of initial conditions:

X0 = array([ 4.1145, 29.1101, 22.3   ])

t = array([  0.        ,  29.48601293,  58.97202586,  88.45803879,
        117.94405173, 147.43006466, 176.91607759, 206.40209052,
        235.88810345, 265.37411638, ... ])

sim = model.simulate(X0, t=t)

However I've been getting a capi_return is NULL Call-back cb_f_in_lsoda__user__routines failed. error, I'll paste the entire error below, it seems likely to be an issue with the scipy solv_ivp integrator, though I also get an error using odeint and haven't been able to figure this out:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[37], line 1
----> 1 sim = model.simulate(X0, t=t)

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/pysindy.py:896, in SINDy.simulate(self, x0, t, u, integrator, stop_condition, interpolator, integrator_kws, interpolator_kws)
    892 # Need to hard-code below, because odeint and solve_ivp
    893 # have different syntax and integration options.
    894 if integrator == "solve_ivp":
    895     return (
--> 896         (solve_ivp(rhs, (t[0], t[-1]), x0, t_eval=t, **integrator_kws)).y
    897     ).T
    898 elif integrator == "odeint":
    899     if integrator_kws.get("method") == "LSODA":

File ~/anaconda3/envs/dodo/lib/python3.9/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 status = None
    590 while status is None:
--> 591     message = solver.step()
    593     if solver.status == 'finished':
    594         status = 0

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py:181, in OdeSolver.step(self)
    179 else:
    180     t = self.t
--> 181     success, message = self._step_impl()
    183     if not success:
    184         self.status = 'failed'

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/scipy/integrate/_ivp/lsoda.py:152, in LSODA._step_impl(self)
    150 itask = integrator.call_args[2]
    151 integrator.call_args[2] = 5
--> 152 solver._y, solver.t = integrator.run(
    153     solver.f, solver.jac or (lambda: None), solver._y, solver.t,
    154     self.t_bound, solver.f_params, solver.jac_params)
    155 integrator.call_args[2] = itask
    157 if solver.successful():

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/scipy/integrate/_ode.py:1344, in lsoda.run(self, f, jac, y0, t0, t1, f_params, jac_params)
   1341     self.acquire_new_handle()
   1342 args = [f, y0, t0, t1] + self.call_args[:-1] + \
   1343        [jac, self.call_args[-1], f_params, 0, jac_params]
-> 1344 y1, t, istate = self.runner(*args)
   1345 self.istate = istate
   1346 if istate < 0:

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py:138, in OdeSolver.__init__.<locals>.fun(t, y)
    136 def fun(t, y):
    137     self.nfev += 1
--> 138     return self.fun_single(t, y)

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py:20, in check_arguments.<locals>.fun_wrapped(t, y)
     19 def fun_wrapped(t, y):
---> 20     return np.asarray(fun(t, y), dtype=dtype)

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/pysindy.py:860, in SINDy.simulate.<locals>.rhs(t, x)
    859 def rhs(t, x):
--> 860     return self.predict(x[np.newaxis, :])[0]

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/pysindy.py:477, in SINDy.predict(self, x, u, multiple_trajectories)
    475     u = validate_control_variables(x, u)
    476     x = [np.concatenate((xi, ui), axis=xi.ax_coord) for xi, ui in zip(x, u)]
--> 477 result = [self.model.predict([xi]) for xi in x]
    478 result = [
    479     self.feature_library.reshape_samples_to_spatial_grid(pred)
    480     for pred in result
    481 ]
    483 # Kept for backwards compatibility.

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/pysindy.py:477, in <listcomp>(.0)
    475     u = validate_control_variables(x, u)
    476     x = [np.concatenate((xi, ui), axis=xi.ax_coord) for xi, ui in zip(x, u)]
--> 477 result = [self.model.predict([xi]) for xi in x]
    478 result = [
    479     self.feature_library.reshape_samples_to_spatial_grid(pred)
    480     for pred in result
    481 ]
    483 # Kept for backwards compatibility.

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/sklearn/utils/metaestimators.py:113, in _AvailableIfDescriptor.__get__.<locals>.<lambda>(*args, **kwargs)
    110         raise attr_err
    112     # lambda, but not partial, allows help() to work with update_wrapper
--> 113     out = lambda *args, **kwargs: self.fn(obj, *args, **kwargs)  # noqa
    114 else:
    116     def fn(*args, **kwargs):

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/sklearn/pipeline.py:469, in Pipeline.predict(self, X, **predict_params)
    467 Xt = X
    468 for _, name, transform in self._iter(with_final=False):
--> 469     Xt = transform.transform(Xt)
    470 return self.steps[-1][1].predict(Xt, **predict_params)

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/base.py:191, in x_sequence_or_item.<locals>.func(self, x, *args, **kwargs)
    189 if isinstance(x, Sequence):
    190     xs = [AxesArray(xi, comprehend_axes(xi)) for xi in x]
--> 191     result = wrapped_func(self, xs, *args, **kwargs)
    192     if isinstance(result, Sequence):  # e.g. transform() returns x
    193         return [AxesArray(xp, comprehend_axes(xp)) for xp in result]

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/base.py:323, in ConcatLibrary.transform(self, x_full)
    321 xp_full = []
    322 for x in x_full:
--> 323     feature_sets = [lib.transform([x])[0] for lib in self.libraries_]
    324     xp = np.concatenate(feature_sets, axis=feature_sets[0].ax_coord)
    326     xp = AxesArray(xp, comprehend_axes(xp))

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/base.py:323, in <listcomp>(.0)
    321 xp_full = []
    322 for x in x_full:
--> 323     feature_sets = [lib.transform([x])[0] for lib in self.libraries_]
    324     xp = np.concatenate(feature_sets, axis=feature_sets[0].ax_coord)
    326     xp = AxesArray(xp, comprehend_axes(xp))

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/base.py:191, in x_sequence_or_item.<locals>.func(self, x, *args, **kwargs)
    189 if isinstance(x, Sequence):
    190     xs = [AxesArray(xi, comprehend_axes(xi)) for xi in x]
--> 191     result = wrapped_func(self, xs, *args, **kwargs)
    192     if isinstance(result, Sequence):  # e.g. transform() returns x
    193         return [AxesArray(xp, comprehend_axes(xp)) for xp in result]

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/base.py:323, in ConcatLibrary.transform(self, x_full)
    321 xp_full = []
    322 for x in x_full:
--> 323     feature_sets = [lib.transform([x])[0] for lib in self.libraries_]
    324     xp = np.concatenate(feature_sets, axis=feature_sets[0].ax_coord)
    326     xp = AxesArray(xp, comprehend_axes(xp))

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/base.py:323, in <listcomp>(.0)
    321 xp_full = []
    322 for x in x_full:
--> 323     feature_sets = [lib.transform([x])[0] for lib in self.libraries_]
    324     xp = np.concatenate(feature_sets, axis=feature_sets[0].ax_coord)
    326     xp = AxesArray(xp, comprehend_axes(xp))

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/base.py:191, in x_sequence_or_item.<locals>.func(self, x, *args, **kwargs)
    189 if isinstance(x, Sequence):
    190     xs = [AxesArray(xi, comprehend_axes(xi)) for xi in x]
--> 191     result = wrapped_func(self, xs, *args, **kwargs)
    192     if isinstance(result, Sequence):  # e.g. transform() returns x
    193         return [AxesArray(xp, comprehend_axes(xp)) for xp in result]

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/base.py:323, in ConcatLibrary.transform(self, x_full)
    321 xp_full = []
    322 for x in x_full:
--> 323     feature_sets = [lib.transform([x])[0] for lib in self.libraries_]
    324     xp = np.concatenate(feature_sets, axis=feature_sets[0].ax_coord)
    326     xp = AxesArray(xp, comprehend_axes(xp))

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/base.py:323, in <listcomp>(.0)
    321 xp_full = []
    322 for x in x_full:
--> 323     feature_sets = [lib.transform([x])[0] for lib in self.libraries_]
    324     xp = np.concatenate(feature_sets, axis=feature_sets[0].ax_coord)
    326     xp = AxesArray(xp, comprehend_axes(xp))

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/base.py:191, in x_sequence_or_item.<locals>.func(self, x, *args, **kwargs)
    189 if isinstance(x, Sequence):
    190     xs = [AxesArray(xi, comprehend_axes(xi)) for xi in x]
--> 191     result = wrapped_func(self, xs, *args, **kwargs)
    192     if isinstance(result, Sequence):  # e.g. transform() returns x
    193         return [AxesArray(xp, comprehend_axes(xp)) for xp in result]

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/feature_library/pde_library.py:425, in PDELibrary.transform(self, x_full)
    422         s[axis] = slice(self.spatiotemporal_grid.shape[axis])
    423         s[-1] = axis
--> 425         derivs = self.differentiation_method(
    426             d=multiindex[axis],
    427             axis=axis,
    428             **self.diff_kwargs,
    429         )._differentiate(derivs, self.spatiotemporal_grid[tuple(s)])
    430 library_derivatives[
    431     ..., library_idx : library_idx + n_features
    432 ] = derivs
    433 library_idx += n_features

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/differentiation/finite_difference.py:252, in FiniteDifference._differentiate(self, x, t)
    250 else:
    251     coeffs = self._coefficients(t)
--> 252     interior = self._accumulate(coeffs, x)
    253 s[self.axis] = slice((self.n_stencil - 1) // 2, -(self.n_stencil - 1) // 2)
    254 x_dot[tuple(s)] = interior

File ~/anaconda3/envs/dodo/lib/python3.9/site-packages/pysindy/differentiation/finite_difference.py:214, in FiniteDifference._accumulate(self, coeffs, x)
    207 # a new axis is introduced after self.axis for the stencil indices
    208 # To contract with the coefficients, roll by -self.axis to put it first
    209 # Then roll back by self.axis to return the order
    210 trans = np.roll(np.arange(len(x.shape) + 1), -self.axis)
    211 return np.transpose(
    212     np.einsum(
    213         "ij...,ij->j...",
--> 214         np.transpose(x[tuple(s)], axes=trans),
    215         np.transpose(coeffs),
    216     ),
    217     np.roll(np.arange(len(x.shape)), self.axis),
    218 )

IndexError: index 1 is out of bounds for axis 0 with size 1

I'm also wondering if anyone knows how to 'predict' on a subsample of a dataset which is a different size than the sample the model is fitted on? For example fitting on time series with size 1000, but scoring on a subsample or window of only 100 datapoints.

[Jacob-Stevens-Haas]

Since I don't know how the model was built, I can't reproduce the issue. Since the traceback goes back into our code, I think it may actually be our error. Try passing X0 = array([[ 4.1145, 29.1101, 22.3 ]]) (i.e. a 2-D array)

To predict on a subsample which is a different size in the time axis, just pass smaller arrays for x and t in model.score or model.predict. I believe that will work.

[Jacob-Stevens-Haas]

The choice of the PDE library is incorrect - you don't have a spatial data, and

pde_lib = ps.PDELibrary(derivative_order=2, ...)

Indicates you want spatial derivatives.

Also,

lastly, when I try to score or predict on a different size than the time axis I get IndexError...

FWIW, predict doesn't take a t argument. If you want to score with all of your timesteps the same but leave the number of timesteps unspecified, you can pass a float as the t argument.

[itsjacobhere]

Could you explain why PDE library is the wrong choice? Would you suggest modelling as a simple ODE? Can you not model a PDE without spatial data here? Why not? The black-scholes is one such PDE without spatial variables, so I would imagine its feasible to discover a similar model. Thanks

[Jacob-Stevens-Haas]

IIUC, I believe the nuance here is that SINDy doesn't model differential equations, it discovers the models. If the model you're discovering is a PDE, that means it needs a grid of data to estimate spatial (or other non-spatial variable) derivatives from data. In the Black-Scholes system, that would mean it needs a price grid for which you've observed data.

There are often clever ways to manipulate equations in PySINDy to do what you want: formulating $\partial V/\partial S = \frac{\partial V/\partial t}{\partial S/\partial t}$ in the Black-Scholes equation, for instance, might allow using SINDy-PI to allow time derivatives in the denominator.

The other reason is that nobody has been motivated to do the research for how to do SINDy on those types of equations, so nobody in our research group knows how it would work or what to engineer.

[itsjacobhere]

Very interesting, thanks for your response! I'll start to think about ways to reformulate the problem, that makes sense about requiring a grid of data for the black-scholes

[Jacob-Stevens-Haas]

Glad I could help! Can I ask what your interest in PySINDy is? Research, industry, grad student, etc?