Use x_dot kwarg in model.fit in PDEs with WeakPDELibrary()

[georgemilosh]

So the problem occurs when I try to compute x_dot myself and use it in PDEs. Here is my code

import numpy as np
import pysindy as ps
from scipy.io import loadmat

data = loadmat("../../data/kuramoto_sivishinky.mat")
time = np.ravel(data["tt"])
x = np.ravel(data["x"])
u = data["uu"]
dt = time[1] - time[0]
u_dot = ps.FiniteDifference(axis=1)._differentiate(u, t=dt)
u = u.reshape(len(x), len(time), 1)
X, T = np.meshgrid(x, time)
XT = np.asarray([X, T]).T
#pde_lib = ps.WeakPDELibrary(spatiotemporal_grid=XT,)
library_functions = [lambda x: x, lambda x: x * x]
library_function_names = [lambda x: x, lambda x: x + x]
pde_lib = ps.WeakPDELibrary(
    library_functions=library_functions,
    function_names=library_function_names,
    derivative_order=4,
    spatiotemporal_grid=XT,
    is_uniform=True,
    K=200,
    include_interaction=True,
    periodic=True,
)

optimizer = ps.STLSQ()
threshold = 2
optimizer = ps.STLSQ(threshold=threshold, alpha=0.05, normalize_columns=True)
model = ps.SINDy(feature_library=pde_lib, optimizer=optimizer)
model.fit(u, x_dot =u_dot)
model.print()

I get

ValueError: operands could not be broadcast together with shapes (200,) (1024,) 

If I try to flatten u_dot I get

ValueError: operands could not be broadcast together with shapes (200,) (257024,) 

I don't get this error if I work with standard kwargs:

import numpy as np
import pysindy as ps
from scipy.io import loadmat

data = loadmat("../../data/kuramoto_sivishinky.mat")
time = np.ravel(data["tt"])
x = np.ravel(data["x"])
u = data["uu"]
dt = time[1] - time[0]
u_dot = ps.FiniteDifference(axis=1)._differentiate(u, t=dt)
u = u.reshape(len(x), len(time), 1)
X, T = np.meshgrid(x, time)
XT = np.asarray([X, T]).T
pde_lib = ps.WeakPDELibrary(spatiotemporal_grid=XT,)

optimizer = ps.STLSQ()
model = ps.SINDy()
model.fit(u, x_dot=u_dot.flatten())
model.print()

Here flatten u helps. So how can I provide my own x_dot info for PDEs?

My package version:
0.1.dev1618+g694c904.d20231219 3.8.18 | packaged by conda-forge | (default, Oct 10 2023, 15:44:36)
[GCC 12.3.0]

[Jacob-Stevens-Haas]

WeakPDELibrary doesn't even use the derivative - even if it did accept x_dot without error, it would be ignored. A good argument for #351.

Does it work for PDELibrary?

[Jacob-Stevens-Haas]

Heads up, LaTeX on github doesn't allow a space between the dollar signs and the text between them. (I edited your post to format the equations)

In general, you need to use ConstrainedSR3 or MIOSR with equality constraints on the terms that are known.

[georgemilosh]

Thanks, I will check them out later. The decision to set time axis as the last (-1) rather than first (0) seems so counter-intuitive to me. Also when I try to apply GridSearchCV to these high dimensional cases (with spatial dimensions rather than time), there are errors because I think GridSearchCV will treat axis=0 as sampling axis, which in case of PDEs corresponds to the spatial dimension. Has anyone tried applying GridSearchCV here to PDEs?

[Jacob-Stevens-Haas]

People probably have in the past. But the time axis isn't last, it's second-to-last, and as I mentioned, it's always accessible with arr.ax_time (or arr.ax_sample once we've calculated features, which is first)

I don't think GridSearchCV assumes anything about axes. It's either what you're telling it, or it's what we're doing internally.

[georgemilosh]

Thanks for the useful explanations. I see you are doing some upgrades of the package which is quite welcome. What I was wondering if it will be possible to create WeakPDELibrary() that is fully spatial. That is to say, it would treat time coordinate as samples and not integrate over them. This would be useful in situations where time is sampled sparsely and only spatial derivatives are of interest. This is again applicable for the regression problem where l.h.s. is something else, not a time derivative.

[Jacob-Stevens-Haas]

If you want to do it yourself, you can't use SINDy class, but you can manually fit/transform the data with a WeakPDELibrary, then send it to any optimizer to solve.