[BUG] Finite difference shape mismatch

NolanBrb commented 1 month ago

When trying to run the basic examples found at https://pysindy.readthedocs.io/en/latest/examples/3_original_paper/example.html , I run into multiple errors. It seems that there are some issues with the finite difference code, with numpy functions called with arrays of the wrong shapes.

A first error occurs when np.linalg.solve(matrices, b) is called to compute finite differences coefficients. I fixed it by changing b for its transpose b.T. However, I run into a second error in the np.einsum() function.

Reproducing code example:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.cm import rainbow
import numpy as np
from scipy.integrate import solve_ivp
from scipy.io import loadmat
from pysindy.utils import linear_damped_SHO
from pysindy.utils import cubic_damped_SHO
from pysindy.utils import linear_3D
from pysindy.utils import hopf
from pysindy.utils import lorenz

import pysindy as ps

# ignore user warnings
import warnings
warnings.filterwarnings("ignore", category=UserWarning)

np.random.seed(1000)  # Seed for reproducibility

# Integrator keywords for solve_ivp
integrator_keywords = {}
integrator_keywords['rtol'] = 1e-12
integrator_keywords['method'] = 'LSODA'
integrator_keywords['atol'] = 1e-12

# Generate training data

dt = 0.01
t_train = np.arange(0, 25, dt)
t_train_span = (t_train[0], t_train[-1])
x0_train = [2, 0]
x_train = solve_ivp(linear_damped_SHO, t_train_span,
                    x0_train, t_eval=t_train, **integrator_keywords).y.T

# Fit the model

poly_order = 5
threshold = 0.05

model = ps.SINDy(
    optimizer=ps.STLSQ(threshold=threshold),
    feature_library=ps.PolynomialLibrary(degree=poly_order),
)
model.fit(x_train, t=dt)
model.print()

Error message #1

ValueError Traceback (most recent call last) Cell In[1], line 45 39 threshold = 0.05 41 model = ps.SINDy( 42 optimizer=ps.STLSQ(threshold=threshold), 43 feature_library=ps.PolynomialLibrary(degree=poly_order), 44 ) ---> 45 model.fit(x_train, t=dt) 46 model.print()

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:246, in SINDy.fit(self, x, t, x_dot, u) 240 u = validate_control_variables( 241 x, 242 u, 243 trim_last_point=(self.discrete_time and x_dot is None), 244 ) 245 self.n_controlfeatures = u[0].shape[u[0].ax_coord] --> 246 x, x_dot = self._process_trajectories(x, t, x_dot) 248 # Append control variables 249 if u is not None:

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:491, in SINDy._process_trajectories(self, x, t, x_dot) 488 x = [xi[:-1] for xi in x] 489 else: 490 x, x_dot = zip( --> 491 *[ 492 self.feature_library.calc_trajectory( 493 self.differentiation_method, xi, ti 494 ) 495 for xi, ti in _zip_like_sequence(x, t) 496 ] 497 ) 498 return x, x_dot

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:492, in (.0) 488 x = [xi[:-1] for xi in x] 489 else: 490 x, x_dot = zip( 491 *[ --> 492 self.feature_library.calc_trajectory( 493 self.differentiation_method, xi, ti 494 ) 495 for xi, ti in _zip_like_sequence(x, t) 496 ] 497 ) 498 return x, x_dot

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\feature_library\base.py:66, in BaseFeatureLibrary.calc_trajectory(self, diff_method, x, t) 65 def calc_trajectory(self, diff_method, x, t): ---> 66 x_dot = diff_method(x, t=t) 67 x = AxesArray(diff_method.smoothedx, x.axes) 68 return x, AxesArray(x_dot, x.axes)

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\base.py:53, in BaseDifferentiation.call(self, x, t) 52 def call(self, x, t=1): ---> 53 return self._differentiate(x, t)

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\finite_difference.py:278, in FiniteDifference._differentiate(self, x, t) 275 if not self.drop_endpoints: 276 # Forward differences on boundary 277 if not self.periodic: --> 278 coeffs = self._coefficients_boundary_forward(t) 279 boundary = self._accumulate(coeffs, x) 281 if self.order % 2 == 0:

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\finite_difference.py:148, in FiniteDifference._coefficients_boundary_forward(self, t) 146 b = np.zeros(self.stencil_inds.shape).T 147 b[:, self.d] = factorial(self.d) --> 148 return np.linalg.solve(matrices, b)

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\numpy\linalg_linalg.py:413, in solve(a, b) 410 signature = 'DD->D' if isComplexType(t) else 'dd->d' 411 with errstate(call=_raise_linalgerror_singular, invalid='call', 412 over='ignore', divide='ignore', under='ignore'): --> 413 r = gufunc(a, b, signature=signature) 415 return wrap(r.astype(result_t, copy=False))

ValueError: solve: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (m,m),(m,n)->(m,n) (size 2 is different from 3)

Fix for error #1 : replace all lines with np.linalg.solve(matrices, b) with np.linalg.solve(matrices, b.T)

Error message #2

AFTER replacing all np.linalg.solve(matrices, b) with np.linalg.solve(matrices, b.T), I get the following error

ValueError Traceback (most recent call last) Cell In[13], line 45 39 threshold = 0.05 41 model = ps.SINDy( 42 optimizer=ps.STLSQ(threshold=threshold), 43 feature_library=ps.PolynomialLibrary(degree=poly_order), 44 ) ---> 45 model.fit(x_train, t=dt) 46 model.print()

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:246, in SINDy.fit(self, x, t, x_dot, u) 240 u = validate_control_variables( 241 x, 242 u, 243 trim_last_point=(self.discrete_time and x_dot is None), 244 ) 245 self.n_controlfeatures = u[0].shape[u[0].ax_coord] --> 246 x, x_dot = self._process_trajectories(x, t, x_dot) 248 # Append control variables 249 if u is not None:

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:491, in SINDy._process_trajectories(self, x, t, x_dot) 488 x = [xi[:-1] for xi in x] 489 else: 490 x, x_dot = zip( --> 491 *[ 492 self.feature_library.calc_trajectory( 493 self.differentiation_method, xi, ti 494 ) 495 for xi, ti in _zip_like_sequence(x, t) 496 ] 497 ) 498 return x, x_dot

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:492, in (.0) 488 x = [xi[:-1] for xi in x] 489 else: 490 x, x_dot = zip( 491 *[ --> 492 self.feature_library.calc_trajectory( 493 self.differentiation_method, xi, ti 494 ) 495 for xi, ti in _zip_like_sequence(x, t) 496 ] 497 ) 498 return x, x_dot

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\feature_library\base.py:66, in BaseFeatureLibrary.calc_trajectory(self, diff_method, x, t) 65 def calc_trajectory(self, diff_method, x, t): ---> 66 x_dot = diff_method(x, t=t) 67 x = AxesArray(diff_method.smoothedx, x.axes) 68 return x, AxesArray(x_dot, x.axes)

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\base.py:53, in BaseDifferentiation.call(self, x, t) 52 def call(self, x, t=1): ---> 53 return self._differentiate(x, t)

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\finite_difference.py:279, in FiniteDifference._differentiate(self, x, t) 277 if not self.periodic: 278 coeffs = self._coefficients_boundary_forward(t) --> 279 boundary = self._accumulate(coeffs, x) 281 if self.order % 2 == 0: 282 right_len = (self.n_stencil - 1) // 2

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\finite_difference.py:227, in FiniteDifference._accumulate(self, coeffs, x) 222 x = AxesArray(x, {"ax_unk": list(range(x.ndim))}) 223 x_expanded = AxesArray( 224 np.transpose(x[tuple(s)], axes=trans), x.insert_axis(0, "ax_offset") 225 ) 226 return np.transpose( --> 227 np.einsum( 228 "ij...,ij->j...", 229 x_expanded, 230 np.transpose(coeffs), 231 ), 232 np.roll(np.arange(x.ndim), self.axis), 233 )

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\numpy_core\einsumfunc.py:1429, in einsum(out, optimize, *operands, *kwargs) 1427 if specified_out: 1428 kwargs['out'] = out -> 1429 return c_einsum(operands, **kwargs) 1431 # Check the kwargs to avoid a more cryptic error later, without having to 1432 # repeat default values here 1433 valid_einsum_kwargs = ['dtype', 'order', 'casting']

ValueError: operand has more dimensions than subscripts given in einstein sum, but no '...' ellipsis provided to broadcast the extra dimensions.

PySINDy/Python version information:

1.7.6.dev465+g69e51e8 3.11.2 | packaged by conda-forge | (main, Mar 31 2023, 17:45:57) [MSC v.1934 64 bit (AMD64)]

chaitanya94 commented 1 week ago

I had the same problem. I tried both the stable version 1.7.5 available from pip as well as compiling from source. The bug is present in both.

An easy way to reproduce the error is:

import numpy as np
import pysindy as ps

t = np.linspace(0, 1, 100)
x = 3 * np.exp(-2 * t)
y = 0.5 * np.exp(t)
X = np.stack((x, y), axis=-1)

model = ps.SINDy(feature_names=["x", "y"])
model.fit(X, t=t)

Solution For pysindy version 1.7.5 the issue is resolved when I downgrade my numpy to version 1.26.4.

Jacob-Stevens-Haas commented 1 week ago

Thanks @chaitanya94! @NolanBrb, can you report your numpy version?

dynamicslab / pysindy