Differentiation step: return derivatives AND smoothed coordinates

[Jacob-Stevens-Haas]

Problem?

I recently added Kalman smoothing to the derivative package, and the next step is figuring out how to also pass the smoothed values back to SINDy to use. This seems like a useful effort, since several of the differentiation methods also smooth the coordinate values, but don't return them. I debated coding this change to demonstrate why the feature would be useful, but in reading through all the pysindy/differentiation/* and the derivative package, it seemed the "right way" to do this isn't straightforward. I'd hope this 🧵 can also be a discussion of refactoring path in order to make this an easier feature to implement, if indeed this is a desirable feature.

Ideal world

Option 1

There's some code smells and design asymmetry in how pysindy uses derivatives. In an ideal world, I could precede this line with

x = self.smoothing_method(x,t)

A mock implementation:

def kalman_smoother_posits(x,t):
    return kalman_smoother(x,t)[0]

def kalman_smoother_velocities(x,t):
    return kalman_smoother(x,t)[1]

@lru_cache(1)
def kalman_smoother(x,t):
    # the implementation

...

x = self.kalman_smoother_posits(x,t)
x_dot = self.kalman_smoother_velocities(x,t)

where the cache is really the only state needed. __init__ would have to detect smoothing_method is None and replace it with the identity function.

Alternatives

Option 2: Classy

For less dramatic change that would replace the cache with standard object-oriented statefulness, we could ask Differentiators/smoothers to add a smoothed_x attribute, if they want. SINDy.fit() would have the lines:

x_dot = self.differentiation_method(x,t)
try:
    x = self.differentiation_method.smoothed_x
except AttributeError:
    pass

Here, the differentiation_method object would be responsible for setting the smoothed_x property. We could obviate the hasattr try-except by modifying BaseDifferentiation.__call__() to set smoothed_x = x. The reason I don't like this option is that PySINDY differentiation objects are stateless objects whose only external API comprises __init__() and __call__(). In other words, they're not really classes. The only nice things about having classes is that you can set the kwargs first and called the function later, but you can do that with functools.partial.

On a related note: SmoothedFiniteDifference objects are a subclass of FiniteDifference. I think this is a case of "composition over inheritance." To illustrate the subclass problem, if I want to apply the smoother to a different differentiation method, or do some other preprocessing, I could end up with combinatorial subclasses. The first method solves this problem, but if going with option 2, I would recommend replacing SmoothedFiniteDifference with SmoothedDifferentiation that handles smoothing and takes another BaseDifferentiation subclass as an argument.

Additional context

• SINDyDerivative doesn't inherit from BaseDifferentiation. I'm not sure why not - it repeats the definition of __call__(), and overrides all the other methods. In addition, because SINDyDerivative calls derivative.dydx() instead of the classes in derivative directly, it can't implement the smoothing in either Option 1 or Option 2. Since dydx() uses a public dictionary of differentiation methods (derivative.methods), we could have SINDyDerivative just use the same dictionary to access the exact class it wants, bypassing dydx().

• one con of option 1: the functional interface might allow nonsensical combinations of smoothing and derivative. E.g. the kalman positions and velocities are returned from the same maximum likelihood expression, so it wouldn't make sense to allow a different smoother. This problem suggests option 3, which is to implement the interface:

x_dot, x = self.differentiation_and_smoothing_method(x,t)

Minimal way of accomplishing this is to modify BaseDifferentiation.__call__() to return self.derivative(x, t), x

Refactoring path??:

Rapid thoughts because I've gotta wrap this up:

1. Make SINDyDerivative wrap the exact class it wants and not rely on dydx.
2. Remove FiniteDifference as the superclass of SmoothedFiniteDifference.
3. Implement Option 2. Stop here if that's all we want.
4. Make BaseDifferentiator do self.differentiate = functools.partial(self.differentiate, **kwargs)
5. implement cacheing and smoothing function of Option 1, but keep inside classes for now.
6. Expose class differentiate methods as module level functions & deprecate the classes. Or not if you want to keep the classes.

Additional, Additional context

• Currently, pysindy.fit() and score() have identical, 55-line blocks of code for calculating derivatives. SINDy objects also have a differentiate() function that is currently unused in the package except in tests, and handles fewer cases than the blocks in fit() and score(). Can I refactor those 55 lines into differentiate() and repurpose that function?

[akaptano]

@briandesilva @znicolaou andgoldschmidt Thoughts here? Zack and I have discussed more generally that pysindy.py needs some reworking, there is a lot of repeated code and the shape checking/verifying is done all over the place.

[znicolaou]

Yes, I've seen that pysindy.py could benefit from some reorganizing, but I wouldn't know what the optimal approach is without spending some time. I think I contributed a lot to is the inconsistent use of self.differentiation_method, SINDy.differentiate. In fact, in the multiple trajectories cases that I hacked together around lines 349 and 832, the class is ignoring the self.differentiation_method and using FiniteDifference._differentiate!

Part of the problem originated when we started allowing shaped input for the PDELibrary and WeakPDELibrary cases. Prior to that, the validate_input and validate_control_input (in utils/base.py) could be applied to the input immediately, but since these functions reshape the input, they forget about the dimension of the system, which would make differentiating along a specific axis impossible for the library classes. We never properly updated the input validation, which I think should be done.

For me, anyone who is willing to spend the time to try to clean things up is welcome to give it a shot!

[znicolaou]

PS: In regards to the problem that @Jacob-Stevens-Haas brought up, would it be possible to apply smoothing to the input in validate_input if a smoothing method is used?

I think we should include a new function for all the preprocessing (validate inputs, optionally smooth, calculate time derivatives, reshape the data and deal with multiple_trajectories, etc), so that SINDy.fit is less bloated.

[Jacob-Stevens-Haas]

PS: In regards to the problem that @Jacob-Stevens-Haas brought up, would it be possible to apply smoothing to the input in validate_input if a smoothing method is used?

I think we should include a new function for all the preprocessing (validate inputs, optionally smooth, calculate time derivatives, reshape the data and deal with multiple_trajectories, etc), so that SINDy.fit is less bloated.

In the interest of flat is better than nested, I tend to prefer this style:

def fit():
    step_1()
    step_2()
    step_3()
    ...

rather than this style

def fit():
    step_1()

def step_1():
    # stuff
    step_2()

def step_2():
    # stuff
    step_3()
...

because the former is easier to test. In this case, that argues against a comprehensive preprocessing function.

That said, I think I've mentally moved towards Option 3, which avoids illegal combinations of smoothing and differentiation by introducing a new abstraction comprised of both a smoothing step (that can be the identity function) and a differentiation step. Not quite sure yet how to keep backwards compatibility with this method though.

---

You also hit on a good point about fit(): it has some heavily nested conditionals to manage the WeakPDELibrary and PDELibrary. Here's what I think would work. It relies on (a) adapting all situations to multiple trajectories to simplify those conditionals, (b) indirection over validate_input and differentiate/convert_udot_integral. Then all cases could be covered by:

SINDy.fit()

...

if not isinstance(x, List):
    x = [x]
    t = [t]
for x_i, t_i in zip(x, t):
    self.feature_library.validate_input(x, t)
    # replace x_dot with more generic "lhs"
    lhs = self.feature_library.calc_trajectory(self.differentiation_method, x, t)
...

WeakPDELibrary.__init__():

self.validate_input = lambda x: x
self.calc_trajectory = lambda differentiation_method, x, t: self.convert_udot_integral(x)

BaseFeatureLibrary.__init__():

self.validate_input = validate_input
self.calc_trajectory = lambda differentiation_method, x, t: differentiation_method.differentiate(x, t)

In general, this should be able to handle all of the weak_libraries-style flags. The method is also compatible with all of the smoothing/differentiation options. e.g. for option 3, calc_trajectory would take a combined smoothing/differentiation object and return the smoothed x along with the lhs.

[Jacob-Stevens-Haas]

I'm gonna build a PR with the indirection stuff from my last comment just to see how much can be cleaned up, but not 100% sure I've fully mapped out all the conditionals.

[Jacob-Stevens-Haas]

Done in #244