dynamicslab / pysindy

A package for the sparse identification of nonlinear dynamical systems from data
https://pysindy.readthedocs.io/en/latest/
Other
1.46k stars 324 forks source link

[BUG] Can we remove support for PDELibrary with no spatial derivatives or grid? #192

Open Jacob-Stevens-Haas opened 2 years ago

Jacob-Stevens-Haas commented 2 years ago

Issue

A PDELibrary with no spatial derivatives is effectively just a CustomLibrary. The situation is only under test in the pytest fixture data_ode_library() - also the only time we create a PDELibrary without a spatial grid. I'm not familiar enough with all the notebooks to know if that's intentional support.

As far as tests go, data_custom_library() is passed as an alternate parametrized fixture every time data_ode_library() is passed. It has the same library terms and behaves the same, but also with three zero features.

Why it matters

This matters because PDELibrary (and in a few other code blocks under if isinstance(feature_library, PDELibrary)) has a lot of boilerplate code dealing with conditions where len(spatial_grid.ndim) < 1. This cruft makes refactoring difficult, e.g. removing redundant code with CustomLibrary by inheritance or composition (e.g. `PDELibrary(library_functions = another_library)).

I'm not trying to do all that, but in order to solve the notebook errors in #185, I've been making axes explicit in all (ok, most) places. I've been simplifying a lot of loops, nesting, and conditionals, and right now I'm working on moving the reshaping step (--> (n_samples, n_features) AFTER the feature library is applied. Currently, the PDELibrary.transform(), WeakPDELibrary.transform(), and utils.remove_random_rows() all include tricky steps to un-reshape the data and re-reshape the data since much of was written to only deal with 2-D arrays (n_time, n_coordinates -> n_samples, n_features). Would love to be able to remove that code (~50 lines) wholesale, but tricky to dance around len(spatial_grid.ndim) < 1

akaptano commented 2 years ago

@Jacob-Stevens-Haas I agree with what you said. I believe this is a relic of older versions of the library (when PDELibrary and WeakPDELibrary were the same class). The no grid option was available for the PDELibrary because we wanted it to be possible to use the weak form while still doing regular SINDy with ODEs (so we should retain this option for the WeakPDELibrary, although maybe a renaming is in order).

akaptano commented 2 years ago

If you have a branch where you have gone through and fixed this, please open a PR. Otherwise, let me know and I will just make the changes.

Jacob-Stevens-Haas commented 2 years ago

:+1: Thanks Alan. A lot of the grid issues (incl the no-grid option) really boils down to keeping the data in its natural shape through to model.fit() -> model.transform(), which is blocked by how we do ensembling in #194. So I'm going to focus on moving ensembling to its own step in the pipeline post-feature_library.transform(). At that place I'll have approached feature_library.transform from both ends, and the final step will be deleting the reshaping bits from feature_library.tranform(). Once that's done, this issue should just be down to making sure PDELibrary initializes a uniform grid when none is passed in the same way as WeakPDELibrary, or explicitly disallowing derivative_order=0, or both.

akaptano commented 2 years ago

Remind me, is this now addressed in #185, and if not, perhaps we should?

Jacob-Stevens-Haas commented 2 years ago

I think this may actually conflict with the move to use the PDELibrary in place of the SINDyPILibrary. I think that was a part of #195 that was merged into #185, but @znicolaou knows more than me about that.

We should probably simplify initialization of the spatial/temporal grid, however - it's currently 50 lines. I can take a look at it, but I'm not too familiar with the changes that were merged in from #194.

If I understand correctly, the role of PI is to use the time derivatives as additional input coordinates to the feature library, and then mask one feature at a time when optimizing. Not 100% sure where or which part of that are being handled by the feature library.

znicolaou commented 2 years ago

Yeah, there is some tension between the new sindypi approach (which did add an enhancement by enabling sindypi for pdes) and this issue.

There may be an approach to incorporate the most important sindypi functionality at the optimizer stage only, so that the libraries don't need to include an implicit_terms option. If we maintain y=x_dot in the Pipeline, the sindypi optimizier could first append y as a new column to x, then fit all possible submatrices where drop one column from x to use as the "new y" (i.e. the left hand side in the fit).

I'm not sure that this would be quite as general as allowing all the higher-order time and mixed space-time derivatives as columns to x. For example, fitting a second order ode would still require reexpressing it as two first-order odes, and the case with mixed space-time derivatives is even more complicated.

As an aside, I'm also not sure whey the optimizer in the sindypi library is restricted to only the SR3 type now. But to summarize, I'm sorry to say that I don't have confidence about what the best approach here is. For now, I think I'd suggest leaving the ode functionality in the weakpdelibrary and pdelibrary, since that solution is working adequately now.

akaptano commented 2 years ago

The SR3 (really cvxpy) use was an old hack — the issue is that sindypi needs to exclude the part of the library that would give the identity model. One way to do this would be to chop a column for each model fit, but it was easier at the time to just use a constraint (via ConstrainedSR3) to disallow the identity model. Another downside of this approach is then only the L1 norm can be used because cvxpy needs a convex regularizer. We should probably fix this at some point since if we could give each variable it’s own sub library then any of the optimizers could be used like in other libraries.

znicolaou commented 2 years ago

Ah, understood about the cvxpy constraint now, thanks! Since I familiarized myself with the sindypi optimizer some already, I'd like to take a stab at overhauling it some time. But that probably can't happen for a few weeks still, since I've got other stuff to work on.

akaptano commented 2 years ago

Yes, my understanding is SINDyPI with L0 should work much better than the current examples we have too :)