using SINDy-PI to find implicit ode sets

[penguinaugustus]

Hi, I am testing sindy with finding implicit ode sets with optimizer pysindy.SINDyPI.

I would expect in this case one model would have three odes like the example of the yeast glycolisis model in this paper: SINDy-PI: A Robust Algorithm for Parallel Implicit Sparse Identification of Nonlinear Dynamics ( https://arxiv.org/abs/2004.02322 ), however, what i get is still one ode for each model.

Model  0
Model  1
...
Model  27
1 = 1.068 z + -0.294 zz + -0.006 xx_dot + -0.001 xxx_dot + 0.030 zzx_dot + -0.024 xy_dot + 0.165 yy_dot + 0.118 zy_dot + -0.033 yyy_dot + -0.141 zzy_dot + 0.029 xz_dot + 0.125 yz_dot + -0.025 xxz_dot + 0.012 yyz_dot + -0.010 zzz_dot
x = 0.147 y + 0.128 xx + 0.175 zz + 0.041 x_dot + 0.114 y_dot + 0.176 z_dot + 0.037 xx_dot + -0.040 yx_dot + -0.005 xxx_dot + 0.020 yyx_dot + 0.394 zy_dot + -0.021 xxy_dot + 0.013 yyy_dot + 0.025 zzy_dot + -0.020 yz_dot + 0.116 zz_dot + -0.032 yyz_dot + 0.012 zzz_dot
...
zzz_dot = -0.459 1 + -1.646 x + 0.886 z + 0.173 xx + -0.818 x_dot + -3.287 z_dot + -0.071 xx_dot + 2.284 zx_dot + 0.011 xxx_dot + 0.009 yyx_dot + -1.198 zzx_dot + 0.435 xy_dot + -0.106 xxy_dot + -0.049 yyy_dot + 0.564 zzy_dot + -0.586 yz_dot + 4.129 zz_dot + 0.208 xxz_dot + -0.325 yyz_dot

I have read the document for SINDy-PI https://pysindy.readthedocs.io/en/latest/examples/9_sindypi_with_sympy/example.html#SINDy-PI-Feature-Overview . However, i didn't find any relevant example on finding ode sets. Am I missing a trick? Thank you for your help!

[Jacob-Stevens-Haas]

Please post a relevant minimal working example.

[penguinaugustus]

Thank you for your reply! Here is the relevant minimal working example:

import numpy as np
import matplotlib.pyplot as plt
import pysindy as ps
import gurobipy
from scipy.integrate import odeint
from scipy.integrate import solve_ivp

### generate simulation data
a = 360
k = 1.368
b = 1
alpha = 1
beta = 0.6
gamma = 1
delta = 0.8
n = 2
def system(state, t):
    x, y, z = state
    dxdt = a / (k**n + z**n) - b * x
    dydt = alpha * x - beta * y
    dzdt = gamma * y - delta * z
    return [dxdt, dydt, dzdt]
x0 = 1
y0 = 1
z0 = 1
t = np.linspace(0, 10,50)
solution = odeint(system, [x0, y0, z0], t)
X = np.stack((solution[:, 0],solution[:, 1],solution[:, 2]), axis=-1)

### create library
# Initialize custom SINDy library
x_library_functions = [
    lambda x: x,
    lambda x: x ** 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: x + x,
    lambda x: x,
]

# Need to pass time base to the library
sindy_library = ps.SINDyPILibrary(
    library_functions=x_library_functions,
    x_dot_library_functions=x_dot_library_functions,
    t=t,
    function_names=library_function_names,
    include_bias=True,
)

sindy_opt = ps.SINDyPI(
    threshold=1e-3,
    tol=1e-8,
    max_iter=20000,
)
### create fit problem
model = ps.SINDy(
    optimizer=sindy_opt,
    feature_library=sindy_library,
    differentiation_method=ps.FiniteDifference(drop_endpoints=True),
    feature_names=[" x"," y"," z"],
)

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

Here i generate data [x,y,z] which looks like this:

They are solved by

[Jacob-Stevens-Haas]

Ah, thanks for the update.

I would expect in this case one model would have three odes...however, what i get is still one ode for each model

I would point you to the docs before the examples. For instance, SINDyPI describes the following argument:

model_subset : np.ndarray, shape(n_models), optional (default None)
   List of indices to compute models for. If list is not provided,
   the default is to compute SINDy-PI models for all possible
   candidate functions.

So SINDyPI fits a regression model for every term in your library. FWIW, the stability of the algorithm can be improved by setting normalize_columns=True. Doing that makes the output more obvious:

 1 = 0.004  x + 0.014  z + 0.001  z x_dot + 0.001  z y_dot + 0.001  z z_dot
 x = 0.006 1 + 0.007  y + 0.011  y_dot + 0.001  z_dot
 y = 0.002 1 + 0.008  z + 0.010  z_dot
 z = 0.027 1 + 0.007  y + -0.008  z_dot
 x x = 0.002  x
 y y = 0.001  x + 0.003  y
 z z = 0.004  z + -0.001  z y_dot
 x_dot = 0.001 1 + 0.005  y_dot + -0.004  z_dot + 0.001  x x_dot
 y_dot = 0.008  x + -0.003  z + 0.001  x_dot
 z_dot = -0.052 1 + 0.011  x + -0.001  x_dot + 0.002  z z_dot
 x x_dot = -0.002 1 + 0.006  x_dot + 0.001  z_dot
 y x_dot = 0.001 1
 z x_dot = 0.019 1 + 0.001  x_dot + 0.002  y_dot + -0.002  z_dot
 x x x_dot = 0.004 1 + -0.012  x_dot + -0.002  z_dot + 0.001  x x_dot
 y y x_dot = -0.002 1 + 0.001  y x_dot
 z z x_dot = 0.016 1
 x y_dot = 0.003  y_dot
 y y_dot = 0.001 1 + 0.001  y_dot
 z y_dot = 0.028 1 + 0.001  y_dot + 0.001  z_dot
 x x y_dot = -0.002 1 + -0.002  x_dot + -0.001  z_dot + 0.001  x y_dot
 y y y_dot = -0.005 1 + 0.001  y y_dot
 z z y_dot = -0.004  y_dot + -0.001  z_dot + 0.002  z y_dot + -0.001  z z_dot
 x z_dot = -0.001 1 + 0.001  z_dot
 y z_dot = 0.001 1 + 0.002  z_dot
 z z_dot = 0.014  z_dot
 x x z_dot = -0.001 1 + -0.001  z_dot + 0.001  x z_dot
 y y z_dot = -0.003  z_dot + 0.001  y z_dot
 z z z_dot = 0.001 1 + 0.001  x_dot + -0.020  z_dot + 0.003  z z_dot

If you don't know a priori which equations are important to fit, It's up to you to examine the equations to see if they contain the system you're looking for. It may require substitution of one equation into another, and I don't believe they are necessarily linearly independent.

Also, in the docstring/docs for SINDyPILibrary, note:

 WARNING: This library is deprecated in PySINDy versions > 1.7. Please
 use the PDE or WeakPDE libraries instead

[penguinaugustus]

Thank you so much for the explanation! Now everything makes sense to me.