Merge pull request #85 from dynamicslab/expanded-derivatives

Expanded derivatives

Author: briandesilva

docs/conf.py

@@ -1,4 +1,3 @@
-import datetime
 import importlib
 import pathlib
 
@@ -8,7 +7,7 @@
 
 # no need to edit below this line
 
-copyright = f"{datetime.datetime.now().year}, {author}"
+copyright = f"2020, {author}"
 
 module = importlib.import_module(project)
 version = release = getattr(module, "__version__")
@@ -24,6 +23,7 @@
     "sphinx.ext.napoleon",
     "sphinx.ext.mathjax",
     "sphinx_nbexamples",
+    "sphinx.ext.intersphinx",
 ]
 
 apidoc_module_dir = f"../{project}"
@@ -68,6 +68,10 @@ def setup(app):
     pattern=".+.ipynb",
 )
 
+intersphinx_mapping = {
+    "derivative": ("https://derivative.readthedocs.io/en/latest/", None)
+}
+
 # -- Extensions to the  Napoleon GoogleDocstring class ---------------------
 # michaelgoerz.net/notes/extending-sphinx-napoleon-docstring-sections.html
 from sphinx.ext.napoleon.docstring import GoogleDocstring  # noqa: E402

examples/1_feature_overview.ipynb

@@ -20,8 +20,8 @@
    "execution_count": 1,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2020-07-14T21:07:35.177225Z",
-     "start_time": "2020-07-14T21:07:33.573951Z"
+     "end_time": "2020-08-22T22:47:55.509529Z",
+     "start_time": "2020-08-22T22:47:54.323781Z"
     }
    },
    "outputs": [],
@@ -44,8 +44,8 @@
    "execution_count": 2,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2020-07-14T21:07:35.285747Z",
-     "start_time": "2020-07-14T21:07:35.184412Z"
+     "end_time": "2020-08-22T22:47:56.708081Z",
+     "start_time": "2020-08-22T22:47:56.618948Z"
     }
    },
    "outputs": [],
@@ -95,8 +95,8 @@
    "execution_count": 3,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2020-07-14T21:07:44.101945Z",
-     "start_time": "2020-07-14T21:07:35.293352Z"
+     "end_time": "2020-08-22T22:50:47.967603Z",
+     "start_time": "2020-08-22T22:50:40.520544Z"
     }
    },
    "outputs": [
@@ -135,6 +135,62 @@
     "search.best_estimator_.print()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Some extra care must be taken when working with differentiation methods from the `derivative` package (i.e. those accessed via the `SINDyDerivative` class). See the example below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2020-08-22T22:59:18.558877Z",
+     "start_time": "2020-08-22T22:58:57.908732Z"
+    }
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Best parameters: {'differentiation_method__kwargs': {'kind': 'spline', 's': 0.01}, 'optimizer__threshold': 0.1}\n",
+      "x0' = -10.000 x0 + 10.000 x1\n",
+      "x1' = 28.003 x0 + -1.001 x1 + -1.000 x0 x2\n",
+      "x2' = -2.667 x2 + 1.000 x0 x1\n"
+     ]
+    }
+   ],
+   "source": [
+    "model = ps.SINDy(\n",
+    "    t_default=dt,\n",
+    "    differentiation_method=ps.SINDyDerivative(kind='spline', s=1e-2)\n",
+    ")\n",
+    "\n",
+    "param_grid = {\n",
+    "    \"optimizer__threshold\": [0.001, 0.01, 0.1],\n",
+    "    \"differentiation_method__kwargs\": [\n",
+    "        {'kind': 'spline', 's': 1e-2},\n",
+    "        {'kind': 'spline', 's': 1e-1},\n",
+    "        {'kind': 'finite_difference', 'k': 1},\n",
+    "        {'kind': 'finite_difference', 'k': 2},\n",
+    "    ]\n",
+    "}\n",
+    "\n",
+    "# This part is identical to what we did before\n",
+    "search = GridSearchCV(\n",
+    "    model,\n",
+    "    param_grid,\n",
+    "    cv=TimeSeriesSplit(n_splits=5)\n",
+    ")\n",
+    "search.fit(x_train)\n",
+    "\n",
+    "print(\"Best parameters:\", search.best_params_)\n",
+    "search.best_estimator_.print()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -145,11 +201,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 5,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2020-07-14T21:07:44.133790Z",
-     "start_time": "2020-07-14T21:07:44.116536Z"
+     "end_time": "2020-08-22T22:59:50.283609Z",
+     "start_time": "2020-08-22T22:59:50.261397Z"
     }
    },
    "outputs": [],
@@ -203,11 +259,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 6,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2020-07-14T21:07:50.098150Z",
-     "start_time": "2020-07-14T21:07:44.137586Z"
+     "end_time": "2020-08-22T22:59:56.316952Z",
+     "start_time": "2020-08-22T22:59:52.131651Z"
     }
    },
    "outputs": [
@@ -256,11 +312,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 7,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2020-07-14T21:07:50.552501Z",
-     "start_time": "2020-07-14T21:07:50.107361Z"
+     "end_time": "2020-08-22T22:59:59.313750Z",
+     "start_time": "2020-08-22T22:59:58.992166Z"
     }
    },
    "outputs": [
@@ -284,11 +340,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 8,
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2020-07-14T21:07:50.660764Z",
-     "start_time": "2020-07-14T21:07:50.577603Z"
+     "end_time": "2020-08-22T23:00:00.375038Z",
+     "start_time": "2020-08-22T23:00:00.333557Z"
     }
    },
    "outputs": [

examples/4_scikit_learn_compatibility.ipynb

@@ -14,6 +14,7 @@
 from .pysindy import SINDy
 from .differentiation import BaseDifferentiation
 from .differentiation import FiniteDifference
+from .differentiation import SINDyDerivative
 from .differentiation import SmoothedFiniteDifference
 from .feature_library import ConcatLibrary
 from .feature_library import CustomLibrary

examples/5_differentation.ipynb

@@ -1,5 +1,12 @@
 from .base import BaseDifferentiation
 from .finite_difference import FiniteDifference
+from .sindy_derivative import SINDyDerivative
 from .smoothed_finite_difference import SmoothedFiniteDifference
 
-__all__ = ["BaseDifferentiation", "FiniteDifference", "SmoothedFiniteDifference"]
+
+__all__ = [
+    "BaseDifferentiation",
+    "FiniteDifference",
+    "SINDyDerivative",
+    "SmoothedFiniteDifference",
+]

pysindy/__init__.py

@@ -13,7 +13,7 @@ class BaseDifferentiation(BaseEstimator):
     Base class for differentiation methods.
 
     Simply forces differentiation methods to implement a
-    _differentiate function.
+    ``_differentiate`` function.
     """
 
     def __init__(self):

pysindy/differentiation/__init__.py

@@ -0,0 +1,78 @@
+"""
+Wrapper classes for differentiation methods from the :doc:`derivative:index` package.
+
+Some default values used here may differ from those used in :doc:`derivative:index`.
+"""
+from derivative import dxdt
+from numpy import arange
+from sklearn.base import BaseEstimator
+
+from pysindy.utils.base import validate_input
+
+
+class SINDyDerivative(BaseEstimator):
+    """
+    Wrapper class for differentiation classes from the :doc:`derivative:index` package.
+    This class is meant to provide all the same functionality as the
+    `dxdt <https://derivative.readthedocs.io/en/latest/api.html\
+        #derivative.differentiation.dxdt>`_ method.
+
+    This class also has ``_differentiate`` and ``__call__`` methods which are
+    used by PySINDy.
+
+    Parameters
+    ----------
+    derivative_kws: dictionary, optional
+        Keyword arguments to be passed to the
+        `dxdt <https://derivative.readthedocs.io/en/latest/api.html\
+        #derivative.differentiation.dxdt>`_
+        method.
+
+    Notes
+    -----
+    See the `derivative documentation <https://derivative.readthedocs.io/en/latest/>`_
+    for acceptable keywords.
+    """
+
+    def __init__(self, **kwargs):
+        self.kwargs = kwargs
+
+    def set_params(self, **params):
+        """
+        Set the parameters of this estimator.
+        Modification of the pysindy method to allow unknown kwargs. This allows using
+        the full range of derivative parameters that are not defined as member variables
+        in sklearn grid search.
+
+        Returns
+        -------
+        self
+        """
+        if not params:
+            # Simple optimization to gain speed (inspect is slow)
+            return self
+        else:
+            self.kwargs.update(params)
+
+        return self
+
+    def get_params(self, deep=True):
+        """Get parameters."""
+        params = super().get_params(deep)
+
+        if isinstance(self.kwargs, dict):
+            params.update(self.kwargs)
+
+        return params
+
+    def _differentiate(self, x, t=1):
+        if isinstance(t, (int, float)):
+            if t < 0:
+                raise ValueError("t must be a positive constant or an array")
+            t = arange(x.shape[0]) * t
+
+        return dxdt(x, t, axis=0, **self.kwargs)
+
+    def __call__(self, x, t=1):
+        x = validate_input(x, t=t)
+        return self._differentiate(x, t)

pysindy/differentiation/base.py

@@ -215,7 +215,12 @@ def transform(self, X):
             to_stack.append(X)
             for deg in range(2, self.degree + 1):
                 Xp_next = _csr_polynomial_expansion(
-                    X.data, X.indices, X.indptr, X.shape[1], self.interaction_only, deg,
+                    X.data,
+                    X.indices,
+                    X.indptr,
+                    X.shape[1],
+                    self.interaction_only,
+                    deg,
                 )
                 if Xp_next is None:
                     break

pysindy/differentiation/sindy_derivative.py

@@ -158,16 +158,14 @@ def disable_trimming(self):
         self.trimming_fraction = None
 
     def _update_full_coef(self, cho, x_transpose_y, coef_sparse):
-        """Update the unregularized weight vector
-        """
+        """Update the unregularized weight vector"""
         b = x_transpose_y + coef_sparse / self.nu
         coef_full = cho_solve(cho, b)
         self.iters += 1
         return coef_full
 
     def _update_sparse_coef(self, coef_full):
-        """Update the regularized weight vector
-        """
+        """Update the regularized weight vector"""
         coef_sparse = self.prox(coef_full, self.threshold)
         self.history_.append(coef_sparse.T)
         return coef_sparse
@@ -184,8 +182,7 @@ def _trimming_grad(self, x, y, coef_full, trimming_array):
         return 0.5 * np.sum(R2, axis=1)
 
     def _convergence_criterion(self):
-        """Calculate the convergence criterion for the optimization
-        """
+        """Calculate the convergence criterion for the optimization"""
         this_coef = self.history_[-1]
         if len(self.history_) > 1:
             last_coef = self.history_[-2]

pysindy/feature_library/polynomial_library.py

@@ -104,25 +104,22 @@ def __init__(
         self.ridge_kw = ridge_kw
 
     def _sparse_coefficients(self, dim, ind, coef, threshold):
-        """Perform thresholding of the weight vector(s)
-        """
+        """Perform thresholding of the weight vector(s)"""
         c = np.zeros(dim)
         c[ind] = coef
         big_ind = np.abs(c) >= threshold
         c[~big_ind] = 0
         return c, big_ind
 
     def _regress(self, x, y):
-        """Perform the ridge regression
-        """
+        """Perform the ridge regression"""
         kw = self.ridge_kw or {}
         coef = ridge_regression(x, y, self.alpha, **kw)
         self.iters += 1
         return coef
 
     def _no_change(self):
-        """Check if the coefficient mask has changed after thresholding
-        """
+        """Check if the coefficient mask has changed after thresholding"""
         this_coef = self.history_[-1].flatten()
         if len(self.history_) > 1:
             last_coef = self.history_[-2].flatten()