Separate DiscreteSINDy as a new class

docs/conf.py

@@ -31,6 +31,7 @@
     "sphinx.ext.mathjax",
     "sphinx.ext.intersphinx",
     "IPython.sphinxext.ipython_console_highlighting",
+    "matplotlib.sphinxext.plot_directive",
 ]
 
 nb_execution_mode = "off"

examples/1_feature_overview/example.py

@@ -114,7 +114,7 @@ def ignore_specific_warnings():
 x_dot_test_predicted = model.predict(x_test)
 
 # Compute derivatives with a finite difference method, for comparison
-x_dot_test_computed = model.differentiation_method._differentiate(x_test, t=dt)
+x_dot_test_computed = model.differentiation_method(x_test, t=dt)
 
 fig, axs = plt.subplots(x_test.shape[1], 1, sharex=True, figsize=(7, 9))
 for i in range(x_test.shape[1]):
@@ -149,30 +149,6 @@ def ignore_specific_warnings():
 
 fig.show()
 
-# %% [markdown]
-# ## Discrete time dynamical system (map)
-
-# %%
-
-
-def f(x):
-    return 3.6 * x * (1 - x)
-
-
-if __name__ != "testing":
-    n_steps = 1000
-else:
-    n_steps = 10
-eps = 0.001  # Noise level
-x_train_map = np.zeros((n_steps))
-x_train_map[0] = 0.5
-for i in range(1, n_steps):
-    x_train_map[i] = f(x_train_map[i - 1]) + eps * np.random.randn()
-model = ps.DiscreteSINDy()
-model.fit(x_train_map, t=1)
-
-model.print()
-
 # %% [markdown]
 # ## Optimization options
 # In this section we provide examples of different parameters accepted by the built-in sparse regression optimizers `STLSQ`, `SR3`, `ConstrainedSR3`, `MIOSR`, `SSR`, and `FROLS`. The `Trapping` optimizer is not straightforward to use; please check out Example 8 for some examples. We also show how to use a scikit-learn sparse regressor with PySINDy.
@@ -782,7 +758,7 @@ def f(x):
 x_dot_test_predicted = model.predict(x_test)
 
 # Compute derivatives with a finite difference method, for comparison
-x_dot_test_computed = model.differentiation_method._differentiate(x_test, t=dt)
+x_dot_test_computed = model.differentiation_method(x_test, t=dt)
 
 fig, axs = plt.subplots(x_test.shape[1], 1, sharex=True, figsize=(7, 9))
 for i in range(x_test.shape[1]):
@@ -905,7 +881,7 @@ def u_fun(t):
 x_dot_test_predicted = model.predict(x_test, u=u_test)
 
 # Compute derivatives with a finite difference method, for comparison
-x_dot_test_computed = model.differentiation_method._differentiate(x_test, t=dt)
+x_dot_test_computed = model.differentiation_method(x_test, t=dt)
 
 fig, axs = plt.subplots(x_test.shape[1], 1, sharex=True, figsize=(7, 9))
 for i in range(x_test.shape[1]):
@@ -1053,62 +1029,6 @@ def u_fun(t):
     model.fit(x_train, t=t)
     model.print()
 
-# %% [markdown]
-# ## SINDy with control parameters (SINDyCP)
-# The control input in PySINDy can be used to discover equations parameterized by control parameters in conjunction with the `ParameterizedLibrary`. We demonstrate on the logistic map
-# $$ x_{n+1} = r x_n(1-x_n)$$
-# which depends on a single parameter $r$.
-
-# %%
-# Iterate the map and drop the initial 500-step transient. The behavior is chaotic for r>3.6.
-if __name__ != "testing":
-    num = 1000
-    N = 1000
-    N_drop = 500
-else:
-    num = 20
-    N = 20
-    N_drop = 10
-r0 = 3.5
-rs = r0 + np.arange(num) / num * (4 - r0)
-xss = []
-for r in rs:
-    xs = []
-    x = 0.5
-    for n in range(N + N_drop):
-        if n >= N_drop:
-            xs = xs + [x]
-        x = r * x * (1 - x)
-    xss = xss + [xs]
-
-plt.figure(figsize=(4, 4), dpi=100)
-for ind in range(num):
-    plt.plot(np.ones(N) * rs[ind], xss[ind], ",", alpha=0.1, c="black", rasterized=True)
-plt.xlabel("$r$")
-plt.ylabel("$x_n$")
-plt.show()
-
-# %% [markdown]
-# We construct a `parameter_library` and a `feature_library` to act on the input data `x` and the control input `u` independently. The `ParameterizedLibrary` is composed of products of the two libraries output features. This enables fine control over the library features, which is especially useful in the case of PDEs like those arising in pattern formation modeling. See this [notebook](https://github.com/dynamicslab/pysindy/blob/master/examples/17_parameterized_pattern_formation/17_parameterized_pattern_formation.ipynb) for examples.
-
-# %%
-# use four parameter values as training data
-rs_train = [3.6, 3.7, 3.8, 3.9]
-xs_train = [np.array(xss[np.where(np.array(rs) == r)[0][0]]) for r in rs_train]
-
-feature_lib = ps.PolynomialLibrary(degree=3, include_bias=True)
-parameter_lib = ps.PolynomialLibrary(degree=1, include_bias=True)
-lib = ps.ParameterizedLibrary(
-    feature_library=feature_lib,
-    parameter_library=parameter_lib,
-    num_features=1,
-    num_parameters=1,
-)
-opt = ps.STLSQ(threshold=1e-1, normalize_columns=False)
-model = ps.DiscreteSINDy(feature_library=lib, optimizer=opt)
-model.fit(xs_train, u=rs_train, t=1, feature_names=["x", "r"])
-model.print()
-
 # %% [markdown]
 # ## PDEFIND Feature Overview
 # PySINDy now supports SINDy for PDE identification (PDE-FIND) (Rudy, Samuel H., Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. "Data-driven discovery of partial differential equations." Science Advances 3, no. 4 (2017): e1602614.). We illustrate a basic example on Burgers' equation:

pyproject.toml

@@ -60,6 +60,7 @@ docs = [
     "sphinx==8.2.3",
     "pyyaml",
     "sphinxcontrib-apidoc",
+    "matplotlib"
 ]
 miosr = [
     "gurobipy>=9.5.1,!=10.0.0"

pysindy/_core.py

@@ -1,4 +1,3 @@
-import sys
 import warnings
 from abc import ABC
 from abc import abstractmethod
@@ -48,11 +47,9 @@ class _BaseSINDy(BaseEstimator, ABC):
 
     feature_library: BaseFeatureLibrary
     optimizer: _BaseOptimizer
-    discrete_time: bool
     model: Pipeline
     # Hacks to remove later
     feature_names: Optional[list[str]]
-    discrete_time: bool = False
     n_control_features_: int = 0
 
     @abstractmethod
@@ -83,9 +80,7 @@ def _fit_shape(self):
 
     def predict(self, x, u=None):
         """
-        Predict the time derivatives if it is a SINDy model.
-        Predict the next state of the system if it is a DiscreteSINDy model.
-
+        Predict the right hand side of the dynamical system
 
         Parameters
         ----------
@@ -100,8 +95,8 @@ def predict(self, x, u=None):
 
         Returns
         -------
-        x_next: array-like or list of array-like, shape (n_samples, n_input_features)
-            Predicted next state of the system
+        result: array-like or list of array-like, shape (n_samples, n_input_features)
+            Predicted right hand side of the dynamical system
         """
         if not _check_multiple_trajectories(x, None, u):
             x, _, _, u = _adapt_to_multiple_trajectories(x, None, None, u)
@@ -147,7 +142,7 @@ def coefficients(self):
         check_is_fitted(self)
         return self.optimizer.coef_
 
-    def equations(self, precision: int = 3, discrete_time=False) -> list[str]:
+    def equations(self, precision: int = 3) -> list[str]:
         """
         Get the right hand sides of the SINDy model equations.
 
@@ -164,10 +159,7 @@ def equations(self, precision: int = 3, discrete_time=False) -> list[str]:
             input feature.
         """
         check_is_fitted(self, "model")
-        if discrete_time:
-            sys_coord_names = [name + "[k]" for name in self.feature_names]
-        else:
-            sys_coord_names = self.feature_names
+        sys_coord_names = self.feature_names
         feat_names = self.feature_library.get_feature_names(sys_coord_names)
 
         def term(c, name):
@@ -227,25 +219,25 @@ class SINDy(_BaseSINDy):
 
     Parameters
     ----------
-    optimizer : optimizer object, optional
+    optimizer
         Optimization method used to fit the SINDy model. This must be a class
         extending :class:`pysindy.optimizers.BaseOptimizer`.
-        The default is :class:`STLSQ`.
+        The default is :class:`pysindy.optimizers.STLSQ`.
 
-    feature_library : feature library object, optional
+    feature_library
         Feature library object used to specify candidate right-hand side features.
         This must be a class extending
         :class:`pysindy.feature_library.base.BaseFeatureLibrary`.
-        The default option is :class:`PolynomialLibrary`.
+        The default option is :class:`pysindy.feature_library.PolynomialLibrary`.
 
-    differentiation_method : differentiation object, optional
+    differentiation_method
         Method for differentiating the data. This must be a class extending
-        :class:`pysindy.differentiation_methods.base.BaseDifferentiation` class.
+        :class:`pysindy.differentiation.base.BaseDifferentiation` class.
         The default option is centered difference.
 
     Attributes
     ----------
-    model : ``sklearn.multioutput.MultiOutputRegressor`` object
+    model : ``sklearn.multioutput.MultiOutputRegressor``
         The fitted SINDy model.
 
     n_input_features_ : int
@@ -745,6 +737,45 @@ class DiscreteSINDy(_BaseSINDy):
        [0.79648209],
        [0.58382694],
        [0.87479097]])
+
+    >>> import numpy as np
+    >>> num = 1000
+    >>> N = 1000
+    >>> N_drop = 500
+    >>> r0 = 3.5
+    >>> rs = r0 + np.arange(num) / num * (4 - r0)
+    >>> xss = []
+    >>> for r in rs:
+    >>>     xs = []
+    >>>     x = 0.5
+    >>>     for n in range(N + N_drop):
+    >>>         if n >= N_drop:
+    >>>             xs = xs + [x]
+    >>>         x = r * x * (1 - x)
+    >>>     xss = xss + [xs]
+    .. plot::
+        >>> import matplotlib.pyplot as plt
+        >>> plt.figure(figsize=(4, 4), dpi=100)
+        >>> for r, xs in zip(r_values, xss):
+        >>>     plt.plot([r]*len(xs), xs, ",", alpha=0.1, c="black", rasterized=True)
+        >>>     plt.xlabel("$r$")
+        >>>     plt.ylabel("$x_n$")
+        >>>     plt.show()
+    >>> rs_train = [3.6, 3.7, 3.8, 3.9]
+    >>> xs_train = [np.array(xss[np.where(np.array(rs) == r)[0][0]]) for r in rs_train]
+    >>> feature_lib = ps.PolynomialLibrary(degree=3, include_bias=True)
+    >>> parameter_lib = ps.PolynomialLibrary(degree=1, include_bias=True)
+    >>> lib = ps.ParameterizedLibrary(
+    >>>     feature_library=feature_lib,
+    >>>     parameter_library=parameter_lib,
+    >>>     num_features=1,
+    >>>     num_parameters=1,
+    >>> )
+    >>> opt = ps.STLSQ(threshold=1e-1, normalize_columns=False)
+    >>> model = ps.DiscreteSINDy(feature_library=lib, optimizer=opt)
+    >>> model.fit(xs_train, u=rs_train, t=1, feature_names=["x", "r"])
+    >>> model.print()
+    (x)[k+1] = 1.000 r[k] x[k] + -1.000 r[k] x[k]^2
     """
 
     def __init__(
@@ -847,6 +878,43 @@ def fit(
 
         return self
 
+    def equations(self, precision: int = 3) -> list[str]:
+        """
+        Get the right hand sides of the DiscreteSINDy model equations.
+
+        Parameters
+        ----------
+        precision: int, optional (default 3)
+            Number of decimal points to include for each coefficient in the
+            equation.
+
+        Returns
+        -------
+        equations: list of strings
+            List of strings representing the DiscreteSINDy model equations for each
+            input feature.
+        """
+        check_is_fitted(self, "model")
+        sys_coord_names = [name + "[k]" for name in self.feature_names]
+        feat_names = self.feature_library.get_feature_names(sys_coord_names)
+
+        def term(c, name):
+            rounded_coef = np.round(c, precision)
+            if rounded_coef == 0:
+                return ""
+            else:
+                return f"{c:.{precision}f} {name}"
+
+        equations = []
+        for coef_row in self.optimizer.coef_:
+            components = [term(c, i) for c, i in zip(coef_row, feat_names)]
+            eq = " + ".join(filter(bool, components))
+            if not eq:
+                eq = f"{0:.{precision}f}"
+            equations.append(eq)
+
+        return equations
+
     def print(self, precision=3, **kwargs):
         """Print the DiscreteSINDy model equations.
 
@@ -857,7 +925,7 @@ def print(self, precision=3, **kwargs):
 
         **kwargs: Additional keyword arguments passed to the builtin print function
         """
-        eqns = self.equations(precision, discrete_time=True)
+        eqns = self.equations(precision)
         feature_names = self.feature_names
         for i, eqn in enumerate(eqns):
             names = f"({feature_names[i]})[k+1]"
@@ -1016,32 +1084,14 @@ def _check_multiple_trajectories(x, x_dot, u) -> bool:
 
     """
     SequenceOrNone = Union[Sequence, None]
-    if sys.version_info.minor < 10:
-        mixed_trajectories = (
-            isinstance(x, Sequence)
-            and (
-                not isinstance(x_dot, Sequence)
-                and x_dot is not None
-                or not isinstance(u, Sequence)
-                and u is not None
-            )
-            or isinstance(x_dot, Sequence)
-            and not isinstance(x, Sequence)
-            or isinstance(u, Sequence)
-            and not isinstance(x, Sequence)
-        )
-    else:
-        mixed_trajectories = (
-            isinstance(x, Sequence)
-            and (
-                not isinstance(x_dot, SequenceOrNone)
-                or not isinstance(u, SequenceOrNone)
-            )
-            or isinstance(x_dot, Sequence)
-            and not isinstance(x, Sequence)
-            or isinstance(u, Sequence)
-            and not isinstance(x, Sequence)
-        )
+    mixed_trajectories = (
+        isinstance(x, Sequence)
+        and (not isinstance(x_dot, SequenceOrNone) or not isinstance(u, SequenceOrNone))
+        or isinstance(x_dot, Sequence)
+        and not isinstance(x, Sequence)
+        or isinstance(u, Sequence)
+        and not isinstance(x, Sequence)
+    )
     if mixed_trajectories:
         raise TypeError(
             "If x, x_dot, or u are a Sequence of trajectories, each must be a Sequence"