Request for advice about shape of control_features when using PDELibrary

[anirbanm93]

I have a data frame of size (123, 2001, 1). I am trying to fit this data into a SINDy model with a control value, u. u has the same shape as the data frame, the following error comes up:

Traceback (most recent call last):
File "SINDy_trial.py", line 100, in
model.fit(my_d1, u=depthOnenm, t=time_arr)
File "/home/ee_student/.local/lib/python3.8/site-packages/pysindy/pysindy.py", line 299, in fit
u = validate_control_variables(
File "/home/ee_student/.local/lib/python3.8/site-packages/pysindy/utils/base.py", line 75, in validate_control_variables
u_arr = _check_control_shape(x, u, trim_last_point)
File "/home/ee_student/.local/lib/python3.8/site-packages/pysindy/utils/base.py", line 98, in _check_control_shape
raise ValueError(
ValueError: control variables u must have same number of rows as x. u has 246123 rows and x has 123 rows

In the pysindy.py, line 296 shows self.n_control_features_ = u.shape[1]. So, I have reshaped my control_feature u to an array of size (246123, 1). In this case, I am getting the same valuerror:

At line 99 of base.py file, I made the following change,

x1 = x.reshape(x.size // x.shape[-1], x.shape[-1])
if len(x1) != u.shape[0]:
    raise ValueError(
        "control variables u must have same number of rows as x. "
        "u has {} rows and x has {} rows".format(u.shape[0], len(x))
    )

And the following error comes up,

Traceback (most recent call last):
File "SINDy_trial.py", line 101, in
model.print()
File "/home/ee_student/.local/lib/python3.8/site-packages/pysindy/pysindy.py", line 676, in print
eqns = self.equations(precision)
File "/home/ee_student/.local/lib/python3.8/site-packages/pysindy/pysindy.py", line 658, in equations
return equations(
File "/home/ee_student/.local/lib/python3.8/site-packages/pysindy/utils/base.py", line 391, in equations
input_features = pipeline.steps[0][1].get_feature_names(input_features)
File "/home/ee_student/.local/lib/python3.8/site-packages/pysindy/feature_library/pde_library.py", line 212, in get_feature_names
self.function_names[i]([input_features[j] for j in c])
File "/home/ee_student/.local/lib/python3.8/site-packages/pysindy/feature_library/pde_library.py", line 212, in
self.function_names[i]([input_features[j] for j in c])
IndexError: list index out of range

Please let me know how to correctly employ the control input in the PDELibrary, i.e. what should be the size of the control array. I would greatly appreciate any help you could provide.

Regards,
Anirban

[akaptano]

Hi Anirban,

Could you explain why the data is shape (123, 2001, 1)? For instance, is the training data a set of 123 trajectories, each of 2001 time points, and in a 1D state space? Or are you trying to fit a 123D system of ODEs? Or is this a spatiotemporal field that is function of 123 spatial grid points and 2001 temporal grid points? Please also post your PySINDy version, since there are some known bugs in the reshaping in the previous few versions.

Best,
Alan

[anirbanm93]

Thanks for your prompt response. As you said, my data is a spatiotemporal field with 123 spatial grid points and 2001 temporal grid points. I want to employ two different control variables separately,i.e., one is constant over the spatiotemporal grid, and another is varying.

I have reshaped my spatiotemporal field to (Nx, Nt, 1) following the PDEFIND example.

Name: pysindy
Version: 1.6.3
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: /home/anirbanm/.local/lib/python3.8/site-packages
Requires: derivative, scipy, matplotlib, numpy, cvxpy, scikit-learn
Required-by:

Regards,
Anirban

[akaptano]

Excellent, so your data appears to be in the correct shape to fit a PDE model with controls. Can you print how you are initializing the pysindy.PDELibrary that you are using for the model? Note that for spatiotemporal fields you need to be using the PDELibrary class

[anirbanm93]

Initialization of SINDy model:

library_functions = [lambda x: x,
lambda x: x * x,
lambda x: x * x * x,
lambda x, y: x * y,
lambda x, y: x * x * y,
lambda x, y: x * y * y]
library_function_names = [lambda x: x,
lambda x: x + x,
lambda x: x + x + x,
lambda x, y: x + y,
lambda x, y: x + x + y,
lambda x, y: x + y + y]
pde_lib = ps.PDELibrary(
library_functions=library_functions,
function_names=library_function_names,
derivative_order=4,
spatial_grid=space_arr,
include_bias=True,
is_uniform=True,
periodic=True
)
print('STLSQ model: ')
optimizer = ps.STLSQ(threshold=0.1, alpha=1e-5)
model = ps.SINDy(differentiation_method=ps.FiniteDifference(), feature_library=pde_lib, optimizer=optimizer, feature_names=["my"])

[akaptano]

What is the shape of space_arr and can you try this with periodic=False?

[anirbanm93]

space_arr has a size of 123.
The same kind of error I am getting.

STLSQ model:
x size: (123, 2001, 1), u size: (123, 2001, 1)
Traceback (most recent call last):
File "SINDy_trial.py", line 100, in
model.fit(my_d1, u=depthOnenm, t=time_arr)
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy/pysindy.py", line 299, in fit
u = validate_control_variables(
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy/utils/base.py", line 75, in validate_control_variables
u_arr = _check_control_shape(x, u, trim_last_point)
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy/utils/base.py", line 98, in _check_control_shape
raise ValueError(
ValueError: control variables u must have same number of rows as x. u has 246123 rows and x has 123 rows

STLSQ model:
x size: (123, 2001, 1), u size: (246123, 1)
Traceback (most recent call last):
File "SINDy_trial.py", line 100, in
model.fit(my_d1, u=depthOnenm, t=time_arr)
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy/pysindy.py", line 299, in fit
u = validate_control_variables(
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy/utils/base.py", line 75, in validate_control_variables
u_arr = _check_control_shape(x, u, trim_last_point)
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy/utils/base.py", line 98, in _check_control_shape
raise ValueError(
ValueError: control variables u must have same number of rows as x. u has 246123 rows and x has 123 rows

[anirbanm93]

time_arr has a size of 2001.

[anirbanm93]

Hi, any suggestions on how to solve this error.

[akaptano]

Yes, this is a legit bug that we will fix in the coming release. For now, you need to take the lines

        if u is None:
            self.n_control_features_ = 0
        else:
            trim_last_point = self.discrete_time and (x_dot is None)
            u = validate_control_variables(
                x,
                u,
                multiple_trajectories=multiple_trajectories,
                trim_last_point=trim_last_point,
            )
            self.n_control_features_ = u.shape[1]

in pysindy.py and move them to right above the lines

        # Set ensemble variables
        self.ensemble = ensemble
        self.library_ensemble = library_ensemble

This will fix the bug. You should now be able to run the following full example:

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

x = np.random.rand(123, 2001, 1)
u = np.random.rand(123, 2001, 1)
space_arr = np.linspace(1, 123, 123)
print(x.shape, u.shape)

library_functions = [lambda x: x,
lambda x: x * x,
lambda x: x * x * x,
lambda x, y: x * y,
lambda x, y: x * x * y,
lambda x, y: x * y * y]
library_function_names = [lambda x: x,
lambda x: x + x,
lambda x: x + x + x,
lambda x, y: x + y,
lambda x, y: x + x + y,
lambda x, y: x + y + y]
pde_lib = ps.PDELibrary(
library_functions=library_functions,
function_names=library_function_names,
derivative_order=4,
spatial_grid=space_arr,
include_bias=True,
is_uniform=True,
periodic=True
)
print('STLSQ model: ')
optimizer = ps.STLSQ(threshold=0, alpha=1e-5)
model = ps.SINDy(differentiation_method=ps.FiniteDifference(), 
                 feature_library=pde_lib, optimizer=optimizer, 
                 feature_names=["x", "u"])
model.fit(x, u=u)
model.print()

[anirbanm93]

Thanks a lot.

[anirbanm93]

@akaptano I tried to modify the same code for multiple trajectories as follows:

from typing import Sequence
t_arr = np.linspace(1, 2001, 2001)
X = [x, x]
U = [u, u]
if isinstance(X, Sequence):
print("X is a list.")

if isinstance(U, Sequence):
print("U is a list.")

model.fit(X, u=U, t=t_arr, multiple_trajectories=True)
model.print()

I am getting the following error:

X is a list.
U is a list.
Traceback (most recent call last):
File "SINDy_trial.py", line 53, in
model.fit(X, u=U, t=t_arr, multiple_trajectories=True)
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy_issue164/pysindy.py", line 419, in fit
u = validate_control_variables(
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy_issue164/utils/base.py", line 60, in validate_control_variables
raise ValueError("x must be a list when multiple_trajectories is True")
ValueError: x must be a list when multiple_trajectories is True

Moreover, when I am trying to simulate the same data from the model from the initial condition, I am facing the following error:

t_arr = np.linspace(1, 2001, 2001)

def ufunc(time):
ind = int(time//1)
return u[:, ind, 0]

sim = model.simulate(x[:, 0, 0], u=ufunc, t=t_arr, interpolator=CubicSpline)

capi_return is NULL
Call-back cb_f_in_lsoda__user__routines failed.
Traceback (most recent call last):
File "SINDy_trial.py", line 61, in
sim = model.simulate(x[:, 0, 0], u=ufunc, t=t_arr, interpolator=CubicSpline)
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy_issue164/pysindy.py", line 1108, in simulate
(solve_ivp(rhs, (t[0], t[-1]), x0, t_eval=t, **integrator_kws)).y
File "/home/anirbanm/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/ivp.py", line 580, in solve_ivp
message = solver.step()
File "/home/anirbanm/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/base.py", line 181, in step
success, message = self._step_impl()
File "/home/anirbanm/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/lsoda.py", line 152, in _step_impl
solver._y, solver.t = integrator.run(
File "/home/anirbanm/.local/lib/python3.8/site-packages/scipy/integrate/_ode.py", line 1346, in run
y1, t, istate = self.runner(*args)
File "/home/anirbanm/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/base.py", line 138, in fun
return self.fun_single(t, y)
File "/home/anirbanm/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/base.py", line 20, in fun_wrapped
return np.asarray(fun(t, y), dtype=dtype)
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy_issue164/pysindy.py", line 1095, in rhs
return self.predict(x[np.newaxis, :], u_fun(t).reshape(1, -1))[
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy_issue164/pysindy.py", line 635, in predict
return self.model.predict(np.concatenate((x, u), axis=1)).reshape(
File "/home/anirbanm/.local/lib/python3.8/site-packages/sklearn/utils/metaestimators.py", line 113, in
out = lambda *args, **kwargs: self.fn(obj, *args, **kwargs) # noqa
File "/home/anirbanm/.local/lib/python3.8/site-packages/sklearn/pipeline.py", line 469, in predict
Xt = transform.transform(Xt)
File "/home/anirbanm/.local/lib/python3.8/site-packages/pysindy_issue164/feature_library/pde_library.py", line 314, in transform
raise ValueError("x shape does not match training shape")
ValueError: x shape does not match training shape

Regards,
Anirban

[akaptano]

Hi Anirban,

Yes this is a similar issue. To my knowledge you are the first person trying to identify PDEs with control inputs with the new functionality. I will implement a fix in the next release. For now you can simply omit multiple_trajectories and just concatenate all your trajectories together (this is what the code does anyways). So now your data should have shape (123, 2001 * n_trajectories, 1) and similar for u. Sorry for the inconvenience.

[akaptano]

I believe this is resolved in the new release, closing this now.

[anirbanm93]

@akaptano
I am facing a similar issue.
Dataframe shape: (122, 8, 2001, 2)
Training data shape: (122, 8, 1000, 2)
Testing data shape: (122, 8, 1001, 2)

After building the SINDy model, when I am trying to simulate the future behaviour with inputs: initial condition of shape (122, 8, 1, 2) and time-array of size 1001, I get the following error:

Traceback (most recent call last):
File "SINDy_code.py", line 92, in
Xfuture = model.simulate(Xtest[:,:,0,:], t=time_arr[ind_t0:])
File "/home/ee_student/.local/lib/python3.8/site-packages/pysindy_issue164/pysindy.py", line 1108, in simulate
(solve_ivp(rhs, (t[0], t[-1]), x0, t_eval=t, **integrator_kws)).y
File "/home/ee_student/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/ivp.py", line 546, in solve_ivp
solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
File "/home/ee_student/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/lsoda.py", line 113, in init
super().init(fun, t0, y0, t_bound, vectorized)
File "/home/ee_student/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/base.py", line 119, in init
self._fun, self.y = check_arguments(fun, y0, support_complex)
File "/home/ee_student/.local/lib/python3.8/site-packages/scipy/integrate/_ivp/base.py", line 17, in check_arguments
raise ValueError("y0 must be 1-dimensional.")
ValueError: y0 must be 1-dimensional.

pysindy package details:
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: /home/ee_student/.local/lib/python3.8/site-packages
Requires: cmake, cvxpy, derivative, matplotlib, numpy, scikit-learn, scipy

[akaptano]

This one is exactly what the error message says. The y0 in model.simulate (using scipy solve_ivp) should be a single point in your state at time t0... so if your state space is 2, it should be size 2. Assuming that your "extra" dimensions are actually spatial dimensions, it would not even make sense to use model.simulate, because this is an ODE solver not a PDE solver!

[anirbanm93]

@akaptano Thanks, I have missed that. To simulate future behavior for spatio-temporal dynamical system, what one should do? And yes, my data-frame is dependent on x, y, t. 2001 sample points are taken along time-axis.

[akaptano]

General PDE solvers are very complicated! PDE solvers are typically chosen based on the type of PDE (for instance, linear PDEs can be solved with Fourier transforms). For these reasons we do not provide any PDE solvers in PySINDy -- you must build or find PDE solvers appropriate for your identified PDE.