Incorporate make_constraints into trapping_sr3.py

[Jacob-Stevens-Haas]

Problem

TrappingSR3 requires the user to provide the constraints on quadratic terms in order to work according to the results in the paper, "Promoting Global Stability in Data-Driven Models of Quadratic Nonlinear Dynamics". The constraints are equation 4/5 in the paper. The example notebook 8 has the function make_constraints(), but It looks like it lines up the constraints wrongly

Example:

I took the code from cell 2 of here and adapted it here, but I'm working to simplify the indexing variables involved so it's easier to read. I don't quite get how the indexes of $Q_{ijk}$ line up with the constraint matrix indexes for a PolynomialLibrary. Since example 8 isn't tested, I want to make sure that this is working before I bring it into the package. To wit:

import pysindy as ps
from pysindy.feature_library.polynomial_library import n_poly_features
import numpy as np

inputs = ["x", "y"]
lib = ps.PolynomialLibrary(2, include_bias=False)
lib.fit(np.zeros((3,2)))
print(lib.get_feature_names(inputs))
constraint_rhs, constraint_lhs = _make_constraints(2)
print(constraint_lhs)

produces:

['x', 'y', 'x^2', 'x y', 'y^2']
[[0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]
 [0. 0. 0. 0. 0. 1. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1. 0. 0. 1. 0. 0.]]

Additional context

I also don't quite understand how equation 4 is derived in the paper. @akaptano , it looks like equation 4 represents $$1/2\frac{d( a^2)}{dt} = a\dot a =0.$$ Since the section discusses conservation of energy, I was thought it would be $$1/2\frac{d(\dot a^2)}{dt} = \dot a \ddot a = 0$$?

[Jacob-Stevens-Haas]

As an addition, the constraints made by example 8's make_constraints() appear to be in neither in "row" or "target" format, since reorder_constraints(constraint_lhs, 2, "row") and reorder_constraints(constraint_lhs, 2, "target") both change it.

[akaptano]

@Jacob-Stevens-Haas there is a small error in the constraints -- but I think we fixed this on @MPeng5 's trapping branch of the pysindy code and have been meaning to start a pull request to get that and some new stuff added. However, we think this was for dimension > 3 systems and it left the trapping example 8 virtually unchanged (and working a bit better fortunately!).

For the example above, there is no error I believe. It looks incorrect because you are not using the custom ordered polynomial that we use in the trapping method. This was a stupid thing that I implemented the constraints in this way, but basically it is expecting a Polynomial Library (a custom library of polynomial terms) that looks like [x, y, xy, x^2, y^2]. Then the constraints are correct:

For instance, the first row is saying the x^2 term should vanish in the 1st equation for x_dot (assuming I'm reading that off correctly), second constraint says the y^2 term should vanish in the equation for y_dot, third constraint says xy term in y_dot must match the y^2 term in x_dot, and fourth constraint says xy term in x_dot must match the x^2 term in y_dot. These are all correct with the reordered library!

For the energy, it is just E = 0.5 * u^2 and plug in Equation (2) and use the orthogonality of the spatial modes and integrate over the volume to get K = 0.5 a^2 as the total time dependent part of the energy in the fluid. Now we are interested dK/dt = a * a_dot and we do not assume dK/dt = 0. We assume a-priori that the dynamics have vanishing cubic contributions to dK/dt, i.e. dK/dt = C_ia_i + L_{ij}a_ia_j + Q_{ijk}a_ia_ja_k, and only the last term that is cubic vanishes.

[Jacob-Stevens-Haas]

Ohhh, that makes sense if the terms are flattened in fortran order/interspersed, like if $\dot x = c_{11} f_1 + c_{12} f_2$ and $\dot y = c_{21} f_1 + c_{22} f_2 $, the constraint order goes $[c_{11}, c_{21}, c_{12}, c_{22}]$!

For the energy, it is just E = 0.5 * u^2 and plug in Equation (2) and use the orthogonality of the spatial modes and integrate over the volume to get K = 0.5 a^2 as the total time dependent part of the energy in the fluid. Now we are interested dK/dt = a * a_dot and we do not assume dK/dt = 0. We assume a-priori that the dynamics have vanishing cubic contributions to dK/dt, i.e. dK/dt = C_ia_i + L_{ij}a_ia_j + Q_{ijk}a_ia_ja_k, and only the last term that is cubic vanishes.

I can follow that derivation, but I guess I'm unsure why for PDEs we use $E = .5 u^2$ and not $E = .5 \dot u^2$. I admit my fluids knowledge isn't great and the latter is just what was in kinematics so many years ago.

I obviously have been trying to use trapping recently, and I made a PR to clean up the code and reduce the number of lines so it was clearer what was the methods did. I don't think there are significant conflicts in functionality with @MPeng5's branch, but there are some cases where I refactored a function to be clearer while y'all expanded the functionality in a way that's algorithmically compatible but would show up as a git conflict. Example: yous, mine, master. It's tricky to tell because the "compare branches" doesn't show the git conflicts, and because there might be more issues if @MPeng5 or you have commits more recent than the August 17 one on git. For instance, it doesn't show that the make_constraint() function was added to pysindy.

I'm not at all wedded to my code, but I would like to be able to develop and work on trapping asynchronously. My spidey sense is that this will be a difficult PR and we should start working through any structural changes sooner rather than later. Could you guys at least push your most recent commits so that I can see what you're working on?

[akaptano]

Sorry it took me a while to get to this. Remember that "u" here is the fluid velocity, not the fluid position, so there is no need to have u_dot there. The kinetic energy in the fluid is just $E = 0.5\rho u^2$ and we are taking the incompressible case with no density time evolution so can ignore $\rho$ safely. Agreed this is a tricky change and I encourage @MPeng5 to take a look at this soon :)

[MPeng5]

@Jacob-Stevens-Haas I did push my latest version a few days ago but I did not make structural change at all -- just played with those examples you have already seen a bit more for paper revisions. But can you share some of the major conflicts you have now? We can actually meet on campus if you are available this or next week to talk about this so that you can show me directly.
edited:
the conflict in the example you mentioned for checking symmetries of all those tensors is because a new tensor is introduced so we changed the algorithm a little bit. And it is exactly what you found -- the old version only works when we have those constraints. As what I replied in the other thread, we were trying to get rid of the make_constraints() so we need to rewrite this part.

[Jacob-Stevens-Haas]

Remember that "u" here is the fluid velocity

🤦‍♂️ Of course. That makes sense.

Hey @MPeng5 , thanks for following up! I'd love to talk about this but I'm not on campus (in London until April). We'll have to Zoom, if that's OK? What's your email?

the conflict in the example you mentioned for checking symmetries of all those tensors is because a new tensor is introduced so we changed the algorithm a little bit

Yeah, it seems simple enough. New tensor, PT, that appears to use the same checks as PQ.

But can you share some of the major conflicts you have now?

Firstly, my edits to make_constraints() are in trapping-constraints, and make_constraints() isn't part of trapping_sr3.py in either my refactor in trapping-helper or your work in trapping-extended. So #437 should be fine to review and merge independently.

Broadly, the conflicts between trapping-helper and trapping-extended are where I found a way to condense the code and make it more readable, whereas you added functionality in the old style (like _check_P_matrix). So my recommended way forwards is to review/merge #435, pull it into trapping-extended, then create a PR for trapping-extended.

• all __init__() validation is now done in __post_init_guard(), which is also called by set_params(), so that we follow the scikit-learn API
• TrappingSR3 is made a subclass of ConstrainedSR3. This removes a lot of cvxpy code from TrappingSR3, since it can just inherit those methods.
• The main implication for conflicts with trapping-extended is that conditional calls to cvxpy are extracted into their own functions, and conditional calls to other solution methods are also extracted for consistency. In other words, there's the main iteration loop -> a function to update either m or coef -> preliminary steps, plus a function to solve an optimization problem (e.g. calling _update_coef_cvxpy(), inherited from superclass).
• It also looks like you removed the relaxed_optim option, and thus the solve_direct methods?
• Some variable names are changed to be more readable, most significantly m -> trap_ctr, r -> n_tgts, N -> n_features
• _set_Ptensors(), like _check_P_matrix(), was condensed and made more readable.
• Other things, like sorting the docstring to match the args might be git conflicts but aren't real conflicts.

It might also help to merge master into trapping-extended first to get caught up on the non-trapping changes to SINDy first.

[Jacob-Stevens-Haas]

@MPeng5 why don't you also open a PR for trapping-extended so it's possible to attach questions to specific lines?

[MPeng5]

@Jacob-Stevens-Haas Of course! Totally forgot Nathan brought a bunch of people to London and you are one of them. Let's try to meet next week on zoom and my email is maipeng@uw.edu.