[BUG] Using SINDy-PI with multiple_trajectories does not work

[BMP-TUD]

Hi everyone,

I wanted to use SINDy-PI with multiple trajectories, similar to the original paper by Kaheman (2020) but in python. However, when I try to fit the list of trajectories in model.fit the function stacks the different trajectories and does not handle them as single trajectories. Also providing the x_dot directly has not worked. I also cannot find an example in your documentation and except model.fit no function has the argument multiple_trajectories.

Thanks in advance for your answer, you can find an example below:

Best,
Bartosz

Reproducing code example:

import pysindy as ps
import numpy as np
import matplotlib.pyplot as plt
from functions import*
from scipy.integrate import solve_ivp

x=sim_calc.x # is a list of trajectories in this case 10 trajectories --> [np.array(), np.array(),...]
x_dot_train=[]
for i_traj in range(len(x)):
     x_dot_train.append(ps.FiniteDifference()._differentiate(x[i_traj], sim_calc.t))
        
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: x ** 8,
lambda x: x ** 9,
lambda x: x ** 10,
lambda x: x ** 11,
]
x_dot_library_functions = [lambda x: x]

library_function_names = [
lambda x: x,
lambda x, y: x + y,
lambda x: x + '²',
lambda x, y, z: x + y + z,
lambda x, y: x + y +'²',
lambda x: x +'³',
lambda x: x +'⁸',
lambda x: x + '⁹',
lambda x: x + '¹⁰',
lambda x: x + '¹¹'
]
    
library_function_names.append(lambda x: x)

sindy_library = ps.SINDyPILibrary(
library_functions=library_functions,
x_dot_library_functions=x_dot_library_functions,
t=sim_calc.t,
function_names=library_function_names,
include_bias=True,
)

sindy_opt = ps.SINDyPI(
    threshold=1e-6,
    tol=1e-8,
    thresholder="l1",
    max_iter=2000,
)

model_pi = ps.SINDy(
    optimizer=sindy_opt,
    feature_library=sindy_library,
)
model_pi.fit(x, x_dot=x_dot_train, multiple_trajectories=multiple_trajectories)
model_pi.print()
        

Error message:

  File "L:\...\spiralwaves_Bartosz\spiral_waves\sindy_pi.py", line 117, in <module>
    model_pi.fit(x, x_dot=x_dot_train, multiple_trajectories=multiple_trajectories)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\pysindy.py", line 563, in fit
    self.model.fit(x, x_dot)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\pipeline.py", line 390, in fit
    Xt = self._fit(X, y, **fit_params_steps)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\pipeline.py", line 348, in _fit
    X, fitted_transformer = fit_transform_one_cached(

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\joblib\memory.py", line 349, in __call__
    return self.func(*args, **kwargs)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\pipeline.py", line 893, in _fit_transform_one
    res = transformer.fit_transform(X, y, **fit_params)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\base.py", line 855, in fit_transform
    return self.fit(X, y, **fit_params).transform(X)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\feature_library\sindy_pi_library.py", line 356, in transform
    x_dot = nan_to_num(self.differentiation_method(x, self.t))

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\differentiation\base.py", line 52, in __call__
    return self._differentiate(x, t)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\differentiation\finite_difference.py", line 254, in _differentiate
    x_dot[tuple(s)] = interior

ValueError: could not broadcast input array from shape (396701,3) into shape (3967028,3)


runcell('Goodwin oscillator', 'L:/.../spiralwaves_Bartosz/spiral_waves/sindy_pi.py')
Traceback (most recent call last):

  File "L:\GBW-0101_DiBS_Backup\spiralwaves_Bartosz\spiral_waves\sindy_pi.py", line 117, in <module>
    model_pi.fit(x, x_dot=x_dot_train, multiple_trajectories=multiple_trajectories)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\pysindy.py", line 563, in fit
    self.model.fit(x, x_dot)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\pipeline.py", line 390, in fit
    Xt = self._fit(X, y, **fit_params_steps)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\pipeline.py", line 348, in _fit
    X, fitted_transformer = fit_transform_one_cached(

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\joblib\memory.py", line 349, in __call__
    return self.func(*args, **kwargs)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\pipeline.py", line 893, in _fit_transform_one
    res = transformer.fit_transform(X, y, **fit_params)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\base.py", line 855, in fit_transform
    return self.fit(X, y, **fit_params).transform(X)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\feature_library\sindy_pi_library.py", line 356, in transform
    x_dot = nan_to_num(self.differentiation_method(x, self.t))

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\differentiation\base.py", line 52, in __call__
    return self._differentiate(x, t)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\differentiation\finite_difference.py", line 254, in _differentiate
    x_dot[tuple(s)] = interior

ValueError: could not broadcast input array from shape (394628,3) into shape (3946298,3)

PySINDy/Python version information:

pysindy: 1.7, sys version: '3.8.12 (default, Oct 12 2021, 03:01:40) [MSC v.1916 64 bit (AMD64)]'

[Jacob-Stevens-Haas]

Full pysindy version, with patch?

Also, please work to minimize your example so it's easier to read and run. You rely on several variables that you haven't defined.

[BMP-TUD]

Hi, I have the full version of pysindy, but what do you mean with patch?
The minimized example is the following:

import pysindy as ps
import numpy as np
from scipy.integrate import solve_ivp

class goodwin3_oscillator():
        def __init__(self,a=1,b=1,c=1,g=0.1,h=0.1,i=0.1,n=10,K=1, arrest=[]):
            self.a=a
            self.b=b
            self.c=c
            self.g=g
            self.h=h
            self.i=i
            self.n=n
            self.K=K
            self.arrest=arrest
        def oscillator(self,t,x):
            u=x[0]
            v=x[1]
            w=x[2]
            eq= [
                self.a*hill_function(self.n, self.K, w)-self.g*u,
                self.b*u-self.h*v,
                self.c*v-self.i*w]
            if self.arrest != []:
                for i in range(len(self.arrest)):
                    eq[self.arrest[i]]=0
            return eq
    
ic_scan=[[0,0.2,2.1],[0.1,0.1,2],[0.05,0.15,2.05] ]# just an example list of different ics

sim_calc=[] #list of multiple trajectories
for ic in ic_scan:
  osci=goodwin3_oscillator(1,1,1,0.1,0.1,0.1,10,1).oscillator
  t=np.arange(0,250,0.01)
  t_span=(t[0],t[-1])
  x=solve_ivp(osci, t_span,ic,t_eval=t).y.T
  sim_calc.append(x)

library_functions = [
        lambda x: x,
        lambda x, y: x * y,
        lambda x: x ** 10,
        lambda x,y: y*x ** 10,
        ]
x_dot_library_functions = [lambda x: x]
        
library_function_names = [
        lambda x: x,
        lambda x, y: x + y,
        lambda x: x + '¹⁰',
        lambda x,y: y + x + '¹⁰'
        ]

sindy_library = ps.SINDyPILibrary(
        library_functions=library_functions,
        x_dot_library_functions=x_dot_library_functions,
        t=sim_calc.t,
        function_names=library_function_names,
        include_bias=True,
        )
        
        sindy_opt = ps.SINDyPI(
            threshold=1e-6,
            tol=1e-8,
            thresholder="l1",
            max_iter=2000,
        )
        
        model_pi = ps.SINDy(
            optimizer=sindy_opt,
            feature_library=sindy_library,
        )
        model_pi.fit(x, multiple_trajectories=True)
        model_pi.print()

The error is still the same:

Traceback (most recent call last):

  File "L:\....\spiralwaves_Bartosz\spiral_waves\sindy_pi.py", line 146, in <module>
    model_pi.fit(x, multiple_trajectories=True)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\pysindy.py", line 563, in fit
    self.model.fit(x, x_dot)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\pipeline.py", line 390, in fit
    Xt = self._fit(X, y, **fit_params_steps)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\pipeline.py", line 348, in _fit
    X, fitted_transformer = fit_transform_one_cached(

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\joblib\memory.py", line 349, in __call__
    return self.func(*args, **kwargs)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\pipeline.py", line 893, in _fit_transform_one
    res = transformer.fit_transform(X, y, **fit_params)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\sklearn\base.py", line 855, in fit_transform
    return self.fit(X, y, **fit_params).transform(X)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\feature_library\sindy_pi_library.py", line 356, in transform
    x_dot = nan_to_num(self.differentiation_method(x, self.t))

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\differentiation\base.py", line 52, in __call__
    return self._differentiate(x, t)

  File "C:\Users\u0149745\Anaconda3\envs\spirals\lib\site-packages\pysindy\differentiation\finite_difference.py", line 254, in _differentiate
    x_dot[tuple(s)] = interior

ValueError: could not broadcast input array from shape (1558,3) into shape (4678,3)

From the ValueError I assume that pysindy is just flattening the three trajectories, thus the length of the arrays does not fit. The main question is: How can I use multiple trajectories with SINDy_PI, what do I have to change to make it work?

Thanks in advance.

[Jacob-Stevens-Haas]

patch version is the third number, e.g. the 5 in 1.7.5. You should be able to see this if you do pip show pysindy

[BMP-TUD]

It is just the 1.7 version, there is no patch version:

Name: pysindy
Version: 1.7
Summary: Sparse Identification of Nonlinear Dynamics
Home-page: https://github.com/dynamicslab/pysindy
Author: Brian de Silva, Kathleen Champion, Markus Quade, Alan Kaptanoglu
Author-email: bdesilva@uw.edu, kpchamp@uw.edu, info@markusqua.de, akaptano@uw.edu
License: MIT
Location: c:\users\u0149745\anaconda3\envs\spirals\lib\site-packages
Requires: scipy, scikit-learn, derivative, cvxpy, numpy, cmake, matplotlib
Required-by: 
Note: you may need to restart the kernel to use updated packages.
WARNING: Ignoring invalid distribution -atplotlib (c:\users\u0149745\anaconda3\envs\spirals\lib\site-packages)

[Jacob-Stevens-Haas]

Some issues with your MWE: indentation errors, hill_function is undefined, sim_calc is a list and has no attribute 't'. Try saving that code in a file like debug.py and running it. Also, if the data itself isn't important, just the data shape, just generate random arrays or an arange, shaped the way you use it. Similarly, if the exact library functions aren't important, just use the identity function.

[Jacob-Stevens-Haas]

It is just the 1.7 version, there is no patch version:

Gotcha... That's an error on our part to release a version without a patch number. Probably 1.7.0, thanks!

Try pip install -U pysindy to upgrade. 1.7.5 is the most recent release

[BMP-TUD]

Some issues with your MWE: indentation errors, hill_function is undefined, sim_calc is a list and has no attribute 't'. Try saving that code in a file like debug.py and running it. Also, if the data itself isn't important, just the data shape, just generate random arrays or an arange, shaped the way you use it. Similarly, if the exact library functions aren't important, just use the identity function.

Got it, sorry about that. The structure of the data is indeed not that important, I think that after what you have written, I just can not use multiple trajectories because it is an older version of the 1.7. (1.7.0) release? I will check the release notes for the current version to ensure that my other code will not be messed up by other changes in the releases.