Obtain regenerable results initializing ps.WeakPDELibrary

[Thirie]

Hi all!

I am using the ps.WeakPDELibrary in combination with the SR3 optimizer for my model, but i am obtaining different results using the model.score() everytime i run it. It varies a lot, with some really good fits, alternated by poorly fit models. I think this has to do with the "If you initialize a new PDE library, a new set of subdomains is randomly chosen" (as stated in the Weak/Integral SINDy Feature Overview example). However, I want to enhance reproducability, thus wanting to eliminate the random factor in the set of subdomains. I've tried setting np.random.seed and random.seed, but have yet to see replicable results. Is there an option to reduce the randomness in the creation of the library?

This is the code I use:
opt = ps.SR3(
threshold=0.005,
thresholder="l1",
max_iter=1000,
normalize_columns=False,
tol=1e-1)

ode_lib = ps.WeakPDELibrary(

    library_functions=library_functions,
    spatiotemporal_grid= spatio_temporal,
    include_bias=True,
    is_uniform=True,
    K=optimal_K,
    H_xt=(optimal_H_t,optimal_H_x)
    
)

model_opt = ps.SINDy(feature_library=ode_lib, optimizer=opt)
model_opt.fit(normalized_transformed_data, t=time_list, multiple_trajectories=True)
print("R2 score", model_opt.score(normalized_transformed_data, t=time_list, multiple_trajectories=True))

Kind regards,
Thijs

[Jacob-Stevens-Haas]

Unfortunately, predict/score are not implemented for WeakPDELibrary. See example 12, in the markdown that discusses "for weak form, model.predict returns a prediction of the weak form of x_dot, rather than a prediction of x_dot."

You need to manually assign the coefficients discovered from the weak form into another model

[Thirie]

Thanks for the clarification Jacob! However, when i assign the coefficients of the weak model to the non-weak model, i obtain a purely flat line using .predict. Any ideas as to why? The fit of my weak model also seems really poor (R2 score of 0.004), which i don't quite understand since i cannot manage to improve it with hyperparameter tuning. My data is EEG data.
I hope you can help me out!

Kind regards,
Thijs

The code for context:

time_vector = time_list[0]
num_electrodes = comp_needed
space_vector = np.arange(num_electrodes)
time_grid, space_grid = np.meshgrid(time_vector, space_vector, indexing='ij')  
spatio_temporal = np.zeros((len(time_vector), len(space_vector), 2))  
spatio_temporal[:,:,0], spatio_temporal[:,:,1] = np.meshgrid(time_vector, space_vector, indexing='ij')

#%% Library functions

library_functions = [lambda x: x, lambda x: x * x, lambda x, y: x * y]
library_function_names = [lambda x: x, lambda x: x + x, lambda x, y: x + y]

#%% Current best WEAK model

ode_lib = ps.WeakPDELibrary(

        library_functions=library_functions,
        function_names=library_function_names,
        spatiotemporal_grid=time_vector,
        include_bias=True,
        is_uniform=True,
        K=6,#Performed het beste!
        
    )
opt = ps.SR3(
    threshold=0.005,
    thresholder="l1",
    max_iter=1000,
    normalize_columns=True,
    tol=1e-1,
)

model = ps.SINDy(feature_library=ode_lib, optimizer=opt)
model.fit(normalized_transformed_data, t=time_list, multiple_trajectories=True)

#%% Generate normal model for plotting

ode_lib2 = ps.PolynomialLibrary(degree=2,include_bias=True)
optimizer = ps.SR3(
    threshold=0.0005, thresholder="l1", max_iter=10000, normalize_columns=True, tol=1e-1
)
original_model = ps.SINDy(feature_library=ode_lib2, optimizer=optimizer)
original_model.fit(normalized_transformed_data, t=time_list, multiple_trajectories=True)
print("Normal form model: ")
original_model.print()
print("R2 score", original_model.score(normalized_transformed_data, t=time_list, multiple_trajectories=True))
print("MSE score",original_model.score(normalized_transformed_data, t=time_list, multiple_trajectories=True,metric=mean_squared_error))

#%% Plotting params
smoother_kws={'window_length': 25,
'polyorder': 3}
dt = time_list[1] - time_list[0]
sfd= ps.SmoothedFiniteDifference(smoother=savgol_filter, smoother_kws=smoother_kws)
u_dot_list = []
#sfd._differentiate(normalized_transformed_data, t=dt)
for data in normalized_transformed_data:
    dt=time_list[0][1]-time_list[0][0]
    #print(dt)
    u_dot = sfd._differentiate(data, t=dt)
    u_dot_list.append(u_dot)
u_dot_stack = np.vstack(u_dot_list)

u_pred = original_model.simulate(normalized_transformed_data[0][1000, :], time_list[0][1000:6000], integrator="odeint")
u_dot_pred = original_model.predict(normalized_transformed_data[0])
#%%
# Plot the ground truth, weak form, and non-weak derivative predictions
feature_names = pca_comp_names
optimizer.coef_ = opt.coef_
u_dot_weak = original_model.predict(normalized_transformed_data[0])

print("Weak form model: ")
print("R2 score Weak model", original_model.score(normalized_transformed_data, t=time_list, multiple_trajectories=True))
print("MSE score Weak model",original_model.score(normalized_transformed_data, t=time_list, multiple_trajectories=True,metric=mean_squared_error))

#%% Actual plotting weak model

# Plot the ground truth, weak form, and non-weak trajectory predictions
u_weak = original_model.simulate(normalized_transformed_data[0][1000, :], time_list[0][1000:6000], integrator="odeint")
plt.figure(figsize=(16, 4))
for i in range(3):
    plt.subplot(1, 3, i + 1)
    plt.plot(time_list[0][1000:6000], normalized_transformed_data[0][1000:6000, i], "c", label=r"$q$ Original")
    plt.plot(time_list[0][1000:6000], u_pred[:, i], "r", label=r"$q$ prediction")
    plt.plot(time_list[0][1000:6000], u_weak[:, i], "b", label=r"$q$ weak form prediction", alpha=0.5)
    #plt.plot(time_list[0][1000:1100], u_clean[:, i], "k", label=r"$q$ without added noise")
    plt.grid(True)
    plt.ylabel(feature_names[i], fontsize=14)
    plt.xlabel("t", fontsize=14)
    plt.ylim(0,1)
    if i == 2:
        plt.legend()

[Jacob-Stevens-Haas]

It's impossible for me to tell why the fit is so poor, as I don't know what EEG data looks like. Could you share a bit of info (shape, plots)?

There is one weird thing here: K=6 is really low. That means you're fitting a regression against six data points.

[Thirie]

I understand. The input data is a list of arrays with shape (3,150000,8), indicating 3 persons, 150000 time points with sample freq 500Hz and 8 electrodes. Since the data is from patients, I'm not comfortable with uploading it to Github. However, I would be glad to transfer an example to you personally if that is OK for you. Would it be possible for you to send me a mail (my email is on my profile). I'll send you some data afterwards with more explanation!

About the K=6, I agree it is very low. I've run some different values and the best R2 scores were obtained with low K values. The scores were calculated by swapping the coefficients of a 'normal' model to circumvent the inability of weakPDElibrary to yield a score.

[Jacob-Stevens-Haas]

Sorry Thirie, I'm not going to hold onto sensitive data. But generally speaking, these kinds of research questions can be tough and depend upon the data. So I can only answer generally

• You shared a lot of code that prints things and plots things, but didn't share the code output or plots. Can you?
• In the above code, you fit model to a WeakPDELibrary, but the evaluation code all uses original_model, which is fit on a PolynomialLibrary. Did you mean to use model instead of original_model for some of the above, or transfer the coefficients?
• How did you determine the R2 score to identify the best K value?
• You've reused an optimizer. Scikit-learn estimators should only be used once; they forget their prior state.