Merge pull request #244 from dynamicslab/kalman

Save DifferentiationMethod smoothed values

examples/5_differentiation/example.py

@@ -1,10 +1,9 @@
-#!/usr/bin/env python
-# coding: utf-8
-# # Differentiators in PySINDy
+# %% [markdown]
+# Differentiators in PySINDy
 #
 # This notebook explores the differentiation methods available in PySINDy. Most of the methods are powered by the [derivative](https://pypi.org/project/derivative/) package. While this notebook explores these methods on temporal data, these apply equally well to the computation of spatial derivatives for SINDy for PDE identification (see example Jupyter notebooks 10 and 12, on PDEs and weak forms).
 # [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/dynamicslab/pysindy/v1.7?filepath=examples/5_differentiation.ipynb)
-# In[1]:
+# %%
 import warnings
 
 import matplotlib.pyplot as plt
@@ -23,11 +22,11 @@
 integrator_keywords["atol"] = 1e-12
 
 from utils import (
-    plot_kws,
-    pal,
     compare_methods,
     print_equations,
     compare_coefficient_plots,
+    plot_sho,
+    plot_lorenz,
 )
 
 if __name__ != "testing":
@@ -50,6 +49,7 @@
     )
 
 
+# %% [markdown]
 # In the cell below we define all the available differentiators. Note that the different options in `SINDyDerivative` all originate from `derivative`.
 #
 # * `FiniteDifference` - First order (forward difference) or second order (centered difference) finite difference methods with the ability to drop endpoints. Does *not* assume a uniform time step. Appropriate for smooth data.
@@ -60,7 +60,7 @@
 # * `trend_filtered` - Use total squared variations to fit the data (computes a global derivative that is a piecewise combination of polynomials of a chosen order). Set `order=0` to obtain the total-variational derivative. Appropriate for noisy data
 # * `spectral` - Compute the spectral derivative of the data via Fourier Transform. Appropriate for very smooth (i.e. analytic) data. There is an in-house PySINDy version for speed but this is also included in the derivative package.
 
-# In[2]:
+# %%
 
 
 diffs = [
@@ -75,21 +75,23 @@
     ("Trend Filtered", ps.SINDyDerivative(kind="trend_filtered", order=0, alpha=1e-2)),
     ("Spectral", ps.SINDyDerivative(kind="spectral")),
     ("Spectral, PySINDy version", ps.SpectralDerivative()),
+    ("Kalman", ps.SINDyDerivative(kind="kalman", alpha=0.05)),
 ]
 
 
+# %% [markdown]
 # ## Compare differentiation methods directly
 # First we'll use the methods to numerically approximate derivatives to measurement data directly, without bringing SINDy into the picture. We'll compare the different methods' accuracies when working with clean data ("approx" in the plots) and data with a small amount of white noise ("noisy").
 
-# In[3]:
+# %%
 
 
 noise_level = 0.01
 
 
 # ### Sine
 
-# In[4]:
+# %%
 
 
 # True data
@@ -100,69 +102,61 @@
 
 # ### Absolute value
 
-# In[5]:
+# %%
 
 
 # Shrink window for Savitzky Golay method
 diffs[3] = (
     "Savitzky Golay",
     ps.SINDyDerivative(kind="savitzky_golay", left=0.1, right=0.1, order=3),
 )
+diffs[8] = ("Kalman", ps.SINDyDerivative(kind="kalman", alpha=0.01))
 
 x, y, y_dot, y_noisy = gen_data_step(noise_level)
 
 axs = compare_methods(diffs, x, y, y_noisy, y_dot)
 plt.show()
 
 
+# %% [markdown]
 # ## Compare differentiators when used in PySINDy
 # We got some idea of the performance of the differentiation options applied to raw data. Next we'll look at how they work as a single component of the SINDy algorithm.
-
+#
 # ### Linear oscillator
 # $$ \frac{d}{dt} \begin{bmatrix}x \\ y\end{bmatrix} = \begin{bmatrix} -0.1 & 2 \\ -2 & -0.1 \end{bmatrix} \begin{bmatrix}x \\ y\end{bmatrix} $$
+#
+# +
+#
 
-# In[6]:
+# %%
 
 
 noise_level = 0.1
 
-
-# In[7]:
+# %%
 
 
 # Generate training data
 
 dt, t_train, x_train, x_train_noisy = gen_data_sho(noise_level, integrator_keywords)
 
-
-# In[8]:
-
-
-fig, ax = plt.subplots(1, 1, figsize=(5, 5))
-
-ax.plot(x_train[:, 0], x_train[:, 1], ".", label="Clean", color=pal[0], **plot_kws)
-ax.plot(
-    x_train_noisy[:, 0],
-    x_train_noisy[:, 1],
-    ".",
-    label="Noisy",
-    color=pal[1],
-    **plot_kws
-)
-
-ax.set(title="Training data", xlabel="$x_0$", ylabel="$x_1$")
-ax.legend()
-plt.show()
-
-
-# In[9]:
-
+# %%
+figure = plt.figure(figsize=[5, 5])
+plot_sho(x_train, x_train_noisy)
 
 # Allow Trend Filtered method to work with linear functions
 diffs[5] = (
     "Trend Filtered",
     ps.SINDyDerivative(kind="trend_filtered", order=1, alpha=1e-2),
 )
+diffs[8] = ("Kalman", ps.SINDyDerivative(kind="kalman", alpha=0.5))
+diffs.append(("Smooth FD, reuse old x", ps.SmoothedFiniteDifference(save_smooth=False)))
+diffs.append(
+    (
+        "Kalman, reuse old x",
+        ps.SINDyDerivative(kind="kalman", alpha=0.5, save_smooth=False),
+    )
+)
 
 equations_clean = {}
 equations_noisy = {}
@@ -187,14 +181,12 @@
     equations_noisy[name] = model.equations()
     coefficients_noisy[name] = model.coefficients()
 
-
-# In[10]:
+# %%
 
 
 print_equations(equations_clean, equations_noisy)
 
-
-# In[11]:
+# %%
 
 
 feature_names = model.get_feature_names()
@@ -206,50 +198,51 @@
 )
 plt.show()
 
-
+# %% [markdown]
+#
+# We can take a look at the smoothed values of x that some differentiation
+# methods implicitly calculate:
+
+# %%
+fig = plt.figure(figsize=[12, 5])
+fig.suptitle("Training Data Coordinates")
+plt.subplot(1, 2, 1)
+ax = plot_sho(x_train, x_train_noisy, diffs[2][1].smoothed_x_)
+ax.set_title("Savitzky-Golay filtered for Smoothed FD method")
+plt.subplot(1, 2, 2)
+ax = plot_sho(x_train, x_train_noisy, diffs[8][1].smoothed_x_)
+ax.set_title("Kalman smoothed")
+
+# %% [markdown]
 # ### Lorenz system
 #
 # $$ \begin{aligned} \dot x &= 10(y-x)\\ \dot y &= x(28 - z) - y \\ \dot z &= xy - \tfrac{8}{3} z, \end{aligned} $$
 #
 
-# In[12]:
+# %%
 
 
 noise_level = 0.5
 
-
-# In[13]:
+# %%
 
 
 # Generate measurement data
 dt, t_train, x_train, x_train_noisy = gen_data_lorenz(noise_level, integrator_keywords)
 
-
-# In[14]:
+# %%
 
 
 fig = plt.figure(figsize=(8, 8))
-ax = fig.add_subplot(1, 1, 1, projection="3d")
-ax.plot(
-    x_train[:, 0], x_train[:, 1], x_train[:, 2], color=pal[0], label="Clean", **plot_kws
-)
-
-ax.plot(
-    x_train_noisy[:, 0],
-    x_train_noisy[:, 1],
-    x_train_noisy[:, 2],
-    ".",
-    color=pal[1],
-    label="Noisy",
-    alpha=0.3,
+plot_lorenz(x_train, x_train_noisy)
+fig.show()
+
+# %%
+diffs[8] = ("Kalman", ps.SINDyDerivative(kind="kalman", alpha=0.0015))
+diffs[10] = (
+    "Kalman, reuse old x",
+    ps.SINDyDerivative(kind="kalman", alpha=0.0015, save_smooth=False),
 )
-ax.set(title="Training data", xlabel="$x$", ylabel="$y$", zlabel="$z$")
-ax.legend()
-plt.show()
-
-
-# In[15]:
-
 
 equations_clean = {}
 equations_noisy = {}
@@ -275,14 +268,12 @@
     equations_noisy[name] = model.equations()
     coefficients_noisy[name] = model.coefficients()
 
-
-# In[16]:
+# %%
 
 
 print_equations(equations_clean, equations_noisy)
 
-
-# In[17]:
+# %%
 
 
 feature_names = model.get_feature_names()
@@ -294,10 +285,20 @@
 )
 plt.show()
 
+# %%
 
-# In[18]:
 
+# %%
+fig = plt.figure(figsize=(16.5, 8))
+fig.suptitle("Training Data Coordinates")
+ax = fig.add_subplot(1, 2, 1, projection="3d")
+ax = plot_lorenz(x_train, x_train_noisy, diffs[2][1].smoothed_x_, ax=ax)
+ax.set_title("Savitzky-Golay filtered for Smoothed FD method")
+ax = fig.add_subplot(1, 2, 2, projection="3d")
+ax = plot_lorenz(x_train, x_train_noisy, diffs[8][1].smoothed_x_, ax=ax)
+ax.set_title("Kalman smoothed")
 
+# %%
 import timeit
 
 N_spectral = np.logspace(1, 8, n_spectral, dtype=int)
@@ -320,8 +321,7 @@
     stop = timeit.default_timer()
     spectral_times[i, 1] = stop - start
 
-
-# In[19]:
+# %%
 
 
 plt.figure()
@@ -333,8 +333,7 @@
 plt.legend()
 plt.show()
 
-
-# In[20]:
+# %%
 
 
 # Check error improves as order increases
@@ -352,6 +351,3 @@
 plt.ylabel("Derivative error")
 plt.xlabel("Finite difference order")
 plt.show()
-
-
-# In[ ]:

examples/5_differentiation/utils.py

@@ -116,3 +116,65 @@ def signed_sqrt(x):
         axs[1, 0].set_ylabel("Noisy", labelpad=10)
 
         fig.tight_layout()
+
+
+def plot_sho(x_train, x_train_noisy, x_smoothed=None):
+    ax = plt.gca()
+    ax.plot(x_train[:, 0], x_train[:, 1], ".", label="Clean", color=pal[0], **plot_kws)
+    ax.plot(
+        x_train_noisy[:, 0],
+        x_train_noisy[:, 1],
+        ".",
+        label="Noisy",
+        color=pal[1],
+        **plot_kws,
+    )
+    if x_smoothed is not None:
+        ax.plot(
+            x_smoothed[:, 0],
+            x_smoothed[:, 1],
+            ".",
+            label="Smoothed",
+            color=pal[2],
+            **plot_kws,
+        )
+
+    ax.set(title="Training data", xlabel="$x_0$", ylabel="$x_1$")
+    ax.legend()
+    return ax
+
+
+def plot_lorenz(x_train, x_train_noisy, x_smoothed=None, ax=None):
+    if ax is None:
+        ax = plt.axes(projection="3d")
+    ax.plot(
+        x_train[:, 0],
+        x_train[:, 1],
+        x_train[:, 2],
+        color=pal[0],
+        label="Clean",
+        **plot_kws,
+    )
+
+    ax.plot(
+        x_train_noisy[:, 0],
+        x_train_noisy[:, 1],
+        x_train_noisy[:, 2],
+        ".",
+        color=pal[1],
+        label="Noisy",
+        alpha=0.3,
+    )
+    if x_smoothed is not None:
+        ax.plot(
+            x_smoothed[:, 0],
+            x_smoothed[:, 1],
+            x_smoothed[:, 2],
+            ".",
+            color=pal[2],
+            label="Smoothed",
+            alpha=0.3,
+        )
+    ax.set(title="Training data", xlabel="$x$", ylabel="$y$", zlabel="$z$")
+    ax.legend()
+    return ax

pyproject.toml

@@ -26,8 +26,8 @@ classifiers = [
 ]
 readme = "README.rst"
 dependencies = [
-    "scikit-learn>=0.23",
-    "derivative",
+    "scikit-learn>=1.1",
+    "derivative>=0.5.4",
 ]
 
 [project.optional-dependencies]