On rudy_algorithm2() in 10_PDEFIND_examples.ipynb

[georgemilosh]

So I tried using rudy_algorithm2() function that is defined in 10_PDEFIND_examples.ipynb although never explicitly used. Are you sure that it can work with PDEs rather than ODEs? For example, I applied it to Kuramoto Sivashinsky example given in the same notebook and got a mismatch in size. Could you please provide an example of its application to PDEs? Also the same from cross-validation. The main problem is when I compute model.score I get errors if the u_test and u_train don't match in size

[Jacob-Stevens-Haas]

You'd have to ask Sam whether it's meant to apply to PDEs as well as ODEs, but given the title, I'd assume so. How is it discussed in the paper?

If you're getting a size/shape mismatch, can you post the error and a minimal working example of how to generate it?

[georgemilosh]

In the supplamentary it seems to be discussed briefly. Nevertheless, since it is provided in the 10_PDEFIND_examples.ipynb I would suggest that some commentary on this should be included because I was not successful in applying the routine rudy_algorithm2() to a PDE there. Here is a minimal example:

import numpy as np
from scipy.io import loadmat
from sklearn.metrics import mean_squared_error
import pysindy as ps

# Ignore matplotlib deprecation warnings
import warnings
warnings.filterwarnings("ignore", category=UserWarning)

# Seed the random number generators for reproducibility
np.random.seed(100)

# Algorithm to scan over threshold values during Ridge Regression, and select
# highest performing model on the test set


def rudy_algorithm2(
    x_train,
    x_test,
    t,
    pde_lib,
    dtol,
    alpha=1e-5,
    tol_iter=25,
    normalize_columns=True,
    optimizer_max_iter=20,
    optimization="STLSQ",
):

    # Do an initial least-squares fit to get an initial guess of the coefficients
    optimizer = ps.STLSQ(
        threshold=0,
        alpha=0,
        max_iter=optimizer_max_iter,
        normalize_columns=normalize_columns,
        ridge_kw={"tol": 1e-10},
    )

    # Compute initial model
    model = ps.SINDy(feature_library=pde_lib, optimizer=optimizer)
    model.fit(x_train, t=t)

    # Set the L0 penalty based on the condition number of Theta
    l0_penalty = 1e-3 * np.linalg.cond(optimizer.Theta)
    coef_best = optimizer.coef_

    # Compute MSE on the testing x_dot data (takes x_test and computes x_dot_test)
    error_best = model.score(
        x_test, metric=mean_squared_error, squared=False
    ) + l0_penalty * np.count_nonzero(coef_best)

    coef_history_ = np.zeros((coef_best.shape[0], 
                              coef_best.shape[1], 
                              1 + tol_iter))
    error_history_ = np.zeros(1 + tol_iter)
    coef_history_[:, :, 0] = coef_best
    error_history_[0] = error_best
    tol = dtol

    # Loop over threshold values, note needs some coding 
    # if not using STLSQ optimizer.
    for i in range(tol_iter):
        if optimization == "STLSQ":
            optimizer = ps.STLSQ(
                threshold=tol,
                alpha=alpha,
                max_iter=optimizer_max_iter,
                normalize_columns=normalize_columns,
                ridge_kw={"tol": 1e-10},
            )
        model = ps.SINDy(feature_library=pde_lib, optimizer=optimizer)
        model.fit(x_train, t=t)
        coef_new = optimizer.coef_
        coef_history_[:, :, i + 1] = coef_new
        error_new = model.score(
            x_test, metric=mean_squared_error, squared=False
        ) + l0_penalty * np.count_nonzero(coef_new)
        error_history_[i + 1] = error_new
        
        # If error improves, set the new best coefficients
        if error_new <= error_best:
            error_best = error_new
            coef_best = coef_new
            tol += dtol
        else:
            tol = max(0, tol - 2 * dtol)
            dtol = 2 * dtol / (tol_iter - i)
            tol += dtol
    return coef_best, error_best, coef_history_, error_history_

# Load data from .mat file
data = loadmat('data/kuramoto_sivishinky.mat')
t = np.ravel(data['tt'])
x = np.ravel(data['x'])
u = data['uu']
dt = t[1] - t[0]
dx = x[1] - x[0]

u = u.reshape(len(x), len(t), 1)

train = range(0, int(len(t) * 0.6))
test = [i for i in np.arange(len(t)) if i not in train]
u_train = u[:, train, :]
u_test = u[:, test, :]
t_train = t[train]
t_test = t[test]

# Define PDE library that is quadratic in u, and 
# fourth-order in spatial derivatives of u.
library_functions = [lambda x: x, lambda x: x * x]
library_function_names = [lambda x: x, lambda x: x + x]
pde_lib = ps.PDELibrary(
    library_functions=library_functions,
    function_names=library_function_names,
    derivative_order=4,
    spatial_grid=x,
    include_bias=True,
    is_uniform=True,
    periodic=True
)

rudy_algorithm2(u_train,
    u_test,
    t_test,
    pde_lib,
    1e-3,
    alpha=1e-5,
    tol_iter=25,
    normalize_columns=True,
    optimizer_max_iter=20,
    optimization="STLSQ",
)

[georgemilosh]

Oh and the error:

(sindy-repo) georgem@amazone:~/pysindy/examples$ python rudy_function.py 
Traceback (most recent call last):
  File "rudy_function.py", line 121, in <module>
    rudy_algorithm2(u_train,
  File "rudy_function.py", line 40, in rudy_algorithm2
    model.fit(x_train, t=t)
  File "/users/cpa/georgem/miniforge3/envs/sindy-repo/lib/python3.8/site-packages/pysindy/pysindy.py", line 233, in fit
    x, x_dot, u = _comprehend_and_validate_inputs(
  File "/users/cpa/georgem/miniforge3/envs/sindy-repo/lib/python3.8/site-packages/pysindy/pysindy.py", line 833, in _comprehend_and_validate_inputs
    x = [comprehend_and_validate(xi, ti) for xi, ti in _zip_like_sequence(x, t)]
  File "/users/cpa/georgem/miniforge3/envs/sindy-repo/lib/python3.8/site-packages/pysindy/pysindy.py", line 833, in <listcomp>
    x = [comprehend_and_validate(xi, ti) for xi, ti in _zip_like_sequence(x, t)]
  File "/users/cpa/georgem/miniforge3/envs/sindy-repo/lib/python3.8/site-packages/pysindy/pysindy.py", line 831, in comprehend_and_validate
    return validate_no_reshape(arr, t)
  File "/users/cpa/georgem/miniforge3/envs/sindy-repo/lib/python3.8/site-packages/pysindy/utils/base.py", line 86, in validate_no_reshape
    raise ValueError("Length of t should match x.shape[-2].")
ValueError: Length of t should match x.shape[-2].

[Jacob-Stevens-Haas]

So broadly speaking, re-using the algorithms in the example notebooks aren't part of the pysindy public API. That said, debugging the code, the algorithm expects you to pass t_train, not t_test. Still raises an error, though: AttributeError: 'STLSQ' object has no attribute 'Theta'.

[Jacob-Stevens-Haas]

It looks like that attribute was removed in late 2021/early 2023, and that Theta was not used for anything after it was assigned, nor was it documented. It used to be assigned in fitting to the most recent x. I'm probably not going to figure this out, though, if you'd like to have a look:

git log -GTheta -- pysindy/optimizers/stlsq.py

^-- This will find all commits to the stlsq module where the diff includes the regex "Theta"

git show --decorate dc27b07 -- pysindy/optimizers/stlsq.py

^-- This will show the diff from a single commit on the stlsq module (in this case, commit dc27b07)