[Question]How model.predict() works for implicit ODEs function?

[niuyushuo]

Hi,

I am using sindy-pi feature in PySINDy to fit implicit ODEs function. Just have a quick question, how does the model.predict() function work for implicit ODE functions? We only give x to this function, but there are x_t terms that exist on the right hand side of the fitted model. If I just check the documentation, it seems like the results are calculated by loop x which is 'result = [self.model.predict([xi]) for xi in x]'. Then how to deal with x_t on the right hand side of the fitted model? Thanks a lot.

[Jacob-Stevens-Haas]

PI is my least familiar part of the repo, but I see why it might not. But it looks like it creates coef_ in the correct shape, so what happens if you give it a try?

[niuyushuo]

PI is my least familiar part of the repo, but I see why it might not. But it looks like it creates coef_ in the correct shape, so what happens if you give it a try?

Thanks for your response. After getting the fitted implicit ODEs model, model.predict() works well. To demonstrate the question, here are two models from my results:

x0_t = -0.031 1 + -0.104 x2 + 0.053 x5 + 0.132 x6 + 0.081 x0x1 + 0.024 x0x0 + -0.843 x5_t + 0.728 x0x0_t
x7_t = 0.015 1 + -0.044 x0 + 0.054 x5 + -0.156 x7 + 0.002 x0x1 + 0.030 x0x0

Follow your suggestion, save the model coefficient by: coefs = model.coefficients() If I let the coefficients corresponding to x7_t (which has no x_t terms on the right hand side) loop the dataset, will get the same results compared to the resutls from model.predict(). It is consistent with the documentation which is 'result = [self.model.predict([xi]) for xi in x]'. Then the question is how model.predict() loop data if there are x_t terms on the right hand side of implicit ODEs model? Thanks.

[Jacob-Stevens-Haas]

I'm sorry, I don't understand everything you've done. How many models did SINDyPI fit in your example? And what do you mean when you say you "let the coefficients corresponding to x7_t loop the dataset?"

[niuyushuo]

I'm sorry, I don't understand everything you've done. How many models did SINDyPI fit in your example? And what do you mean when you say you "let the coefficients corresponding to x7_t loop the dataset?"

Hi thanks for your response. For the first question, my dataset has 8 features. The library functions are defined below:

library_functions = [
    lambda x: x,
    lambda x, y: x * y,
    lambda x: x ** 2,
]
x_dot_library_functions = [lambda x: x]

library_function_names = [
    lambda x: x,
    lambda x, y: x + y,
    lambda x: x + x,
] 

So I got a total of 405 models (model0 -model404). For the second question, I picked up model 45 and model 53 which are

x0_t = -0.031 1 + -0.104 x2 + 0.053 x5 + 0.132 x6 + 0.081 x0x1 + 0.024 x0x0 + -0.843 x5_t + 0.728 x0x0_t
x7_t = 0.015 1 + -0.044 x0 + 0.054 x5 + -0.156 x7 + 0.002 x0x1 + 0.030 x0x0

The reason I picked up those two as the first model has x_t terms on the right hand side and the second one does not. I am using those two models to try to understand how model. predict () works. As described in the documentation, the model.predict() works as 'result = [self.model.predict([xi]) for xi in x]'. https://github.com/dynamicslab/pysindy/blob/master/pysindy/pysindy.py, at line 274 If I get the coefficients from x7_t, and loop the dataset. The results are consistent with the results model.predict(dataset)[:,53]. Then the question is just like the title, how model.predict() works if the fitted model (just like model x0_t here) has x_t terms on the right hand side? We only give x (to be input) to this function instead of x_t. I test model.predict(dataset) by myself, and it works without any issues. Thanks.

[Jacob-Stevens-Haas]

Which library are you using? PDELibrary, or CustomLibrary?

[niuyushuo]

It is PDELibrary. Here is my code:

sindy_library = ps.PDELibrary(
    library_functions=library_functions,
    temporal_grid=t,
    function_names=library_function_names,
    include_bias=True,
    implicit_terms=True,
    derivative_order=1
)


sindy_opt = ps.SINDyPI(      ### lambda = threshold^2 / (2 * nu)
    threshold=1e-2,
    tol=1e-8,
    thresholder="l1",
    max_iter=20000,
)

model = ps.SINDy(
    optimizer=sindy_opt,
    feature_library=sindy_library,
    differentiation_method=ps.SINDyDerivative(kind="kalman", alpha=0.6)
)

model.fit(X_sm_poly, t=t)
model.print()

sindy_library.get_feature_names()

[Jacob-Stevens-Haas]

My guess is that predict is giving you bad outputs. I don't think you can trust what you get out. My guess is that SINDyPI is reordering the features internally between each model it's fitting.

Also, you can format your code in Github comments. Heads up, I did this for a few of your comments.

[niuyushuo]

OK, thanks for your reply. I just wonder, for the multiple features implicit ODE model, is there any function/way in pysindy to evaluate each fitted model? If I follow the example 9_sindypi_with_sympy.ipynb. The method from 'Check all the model fits produce sensible models' to 'integrate and plot the new equations for x_dot' part, only works for single feature ODE model. I tested by myself, this method could not work for multiple features. Thanks.

[Jacob-Stevens-Haas]

I'm sorry, you would need to troubleshoot further. That code involves heavy use of sympy, which I've never used. Specifically, can you see where pysindy provides something unexpected?

[mathias-marques]

Hello, I am currently working on identifying a multiple-feature implicit Ordinary Differential Equation (ODE) model, specifically targeting the Rosenzweig-MacArthur predator-prey model. It appears that employing SINDy with a SINDyPI optimizer has yielded some promising models that seem to capture the dynamics well. However, when attempting to simulate the system for various initial conditions, I encounter consistent failures, resulting in an "out of bounds" error. I've searched extensively but haven't been able to locate any examples demonstrating resolutions for a multi-feature implicit Ordinary Differential Equation (ODE) model with sindy.

Thanks a lot !

[Jacob-Stevens-Haas]

Are you using model.predict()? If not, it's a new issue. Regardless, I'd need to know what is meant by "resolutions", what is the exact error you get, and do you have a minimal working example?

[mathias-marques]

No i am trying to use simulate. Here is my minimal example sorry it might be a little bit to big but it is functional.

integrator_keywords = {}
integrator_keywords["rtol"] = 1e-12
integrator_keywords["method"] = "LSODA"
integrator_keywords["atol"] = 1e-12
# Parameter of the model
r=1
K=10
c=1
h=2
b=2
m=1
params_rosmac=(r, K, c, h, b, m)

#equation
def equa_rosmac(t,etat, r, K, c, h, b, m): 
    x,y = etat          # on recupere l'etat
    etat_dot = [r*x*(1-x/K) - c*x/(h+x)*y, #dx
                b*x/(h+x)*y - m*y]  # dy
    return etat_dot

# Time 
dt=0.001
T=40
t = np.arange(0, T + dt, dt)
t_span = (t[0], t[-1])

# Resolution EDO to create data for SINDy
sol = solve_ivp(equa_rosmac,t_span, [1, 2.5],t_eval=t, args=params_rosmac,dense_output=True,**integrator_keywords)
print(sol.y.T.shape)
print(t.shape)
x_train=sol.y.T

# SINDy IMPLEMENTATION
# Initialize custom SINDy library so that we can have x_dot inside it. 
library_functions = [
    lambda x: x,
    lambda x, y: x * y,
    lambda x: x ** 2,
    lambda x, y, z: x * y * z,
    lambda x, y: x * y ** 2,
    lambda x: x ** 3,
    lambda x, y, z, w: x * y * z * w,
    lambda x, y, z: x * y * z ** 2,
]
x_dot_library_functions = [lambda x: x]

# library function names includes both 
# the x_library_functions and x_dot_library_functions names
library_function_names = [
    lambda x: x,
    lambda x, y: x + y,
    lambda x: x + x,
    lambda x, y, z: x + y + z,
    lambda x, y: x + y + y,
    lambda x: x + x + x,
    lambda x, y, z, w: x + y + z + w,
    lambda x, y, z: x + y + z + z,
    lambda x: x,
]
# Need to pass time base to the library so can build the x_dot library from x
sindy_library = ps.SINDyPILibrary(
    library_functions=library_functions,
    x_dot_library_functions=x_dot_library_functions,
    t=t,
    function_names=library_function_names,
    include_bias=True,
)



sindy_opt = ps.optimizers.sindy_pi.SINDyPI()

model_rosmac = ps.SINDy(
    optimizer=sindy_opt,
    feature_library=sindy_library,
    differentiation_method=ps.FiniteDifference(drop_endpoints=True),
    feature_names=["x","y"],
)
print(x_train.shape)
model_rosmac.fit(x_train, t=t,)
model_rosmac.print()

#Time for simulations
dt=0.001
T=50
tpred = np.arange(0, T + dt, dt)
t_span_pred = (t[0], t[-1])
model_rosmac.simulate([2,4],t)

#simulation
dt=0.001
T=50
tpred = np.arange(0, T + dt, dt)
t_span_pred = (t[0], t[-1])
model_rosmac.simulate([2,4],t)

And it gives me this error :

IndexError                                Traceback (most recent call last)
Cell In[4], [line 6](vscode-notebook-cell:?execution_count=4&line=6)
      [4](vscode-notebook-cell:?execution_count=4&line=4) tpred = np.arange(0, T + dt, dt)
      [5](vscode-notebook-cell:?execution_count=4&line=5) t_span_pred = (t[0], t[-1])
----> [6](vscode-notebook-cell:?execution_count=4&line=6) model_rosmac.simulate([2,4],t)

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

File [c:\Users\marquesm\AppData\Local\anaconda3\envs\RL\Lib\site-packages\scipy\integrate\_ivp\ivp.py:602](file:///C:/Users/marquesm/AppData/Local/anaconda3/envs/RL/Lib/site-packages/scipy/integrate/_ivp/ivp.py:602), in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options)
    [600](file:///C:/Users/marquesm/AppData/Local/anaconda3/envs/RL/Lib/site-packages/scipy/integrate/_ivp/ivp.py:600) status = None
    [601](file:///C:/Users/marquesm/AppData/Local/anaconda3/envs/RL/Lib/site-packages/scipy/integrate/_ivp/ivp.py:601) while status is None:
--> [602](file:///C:/Users/marquesm/AppData/Local/anaconda3/envs/RL/Lib/site-packages/scipy/integrate/_ivp/ivp.py:602)     message = solver.step()
    [604](file:///C:/Users/marquesm/AppData/Local/anaconda3/envs/RL/Lib/site-packages/scipy/integrate/_ivp/ivp.py:604)     if solver.status == 'finished':
    [605](file:///C:/Users/marquesm/AppData/Local/anaconda3/envs/RL/Lib/site-packages/scipy/integrate/_ivp/ivp.py:605)         status = 0

File [c:\Users\marquesm\AppData\Local\anaconda3\envs\RL\Lib\site-packages\scipy\integrate\_ivp\base.py:197](file:///C:/Users/marquesm/AppData/Local/anaconda3/envs/RL/Lib/site-packages/scipy/integrate/_ivp/base.py:197), in OdeSolver.step(self)
...
    [216](file:///C:/Users/marquesm/AppData/Local/anaconda3/envs/RL/Lib/site-packages/pysindy/differentiation/finite_difference.py:216)     ),
    [217](file:///C:/Users/marquesm/AppData/Local/anaconda3/envs/RL/Lib/site-packages/pysindy/differentiation/finite_difference.py:217)     np.roll(np.arange(len(x.shape)), self.axis),
    [218](file:///C:/Users/marquesm/AppData/Local/anaconda3/envs/RL/Lib/site-packages/pysindy/differentiation/finite_difference.py:218) )

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

thank you for your quick reply and your help

[Jacob-Stevens-Haas]

Yeah, I think SINDyPI does not provide the ability to simulate. You would need to extract the coef_ that correspond to the model you want, rearrange terms symbolically into a CustomLibrary, and re-assign SINDy.feature_library and SINDy.optimizer.coef_ in order to use either SINDy.predict or SINDy.simulate

[mathias-marques]

Okay thank you for your help I will work on that !