Regselect: accept solver without explicitly setting regularization beforehand (#78)

* fix chicken-egg problem for regselect via _fit_solver()

* ensure model abstracts used for regselect

* DerivativeEstimator.mask(), restructure tests

* use ddt_estimator.mask() in fit_regselect_continuous() for error calc

* generalize mask() to any dimensional array

* fix mask() usage in fit_regselect_continuous()

* regression tests, converted from tutorial notebooks

* add pytest, pytest-cov to dependencies for developers

* ensure correct alignment of train/test data in fit_reselect_continuous()

* version 0.5.14 -> 0.5.15

Author: shanemcq18

docs/source/opinf/changelog.md

@@ -5,10 +5,18 @@
 New versions may introduce substantial new features or API adjustments.
 :::
 
+## Version 0.5.15
+
+- Improvement to `fit_regselect_*()` so that the regularization does not have to be initialized before fitting the model. This fixes a longstanding chicken/egg problem and makes using `fit_regselect_*()` much less cumbersome.
+- Ensured error computation correctly aligns training and predicted states in  `fit_regselect_continuous()`. Previously, this could have been misaligned when using a time derivative estimator that truncates states.
+- Time derivative estimators in `opinf.ddt` now have a `mask()` method that map states to the estimation grid.
+- Added regression tests based on the tutorial notebooks.
+- Added `pytest` and `pytest-cov` to the dependencies for developers.
+
 ## Version 0.5.14
 
 - Catch any errors in `fit_regselect*()` that occur when the model uses `refit()`.
-- Tikhonov-type least-squares solvers do not require the regularizer in the constructor but will raise an `AttributeError` in `solve()` (and other methods) if the regularizer is not set. This makes using `fit_regselect_*()` much less cumbersome.
+- Tikhonov-type least-squares solvers do not require the regularizer in the constructor but will raise an `AttributeError` in `solve()` (and other methods) if the regularizer is not set.
 - `PODBasis.fit(Q)` raises a warning when using the `"method-of-snapshots"`/`"eigh"` strategy if $n < k$ for $\mathbf{Q}\in\mathbb{R}^{n \times k}.$ In this case, calculating the $n \times k$ SVD is likely more efficient than the $k \times k$ eigenvalue problem.
 - Added Python 3.13 to list of tests.
 

pyproject.toml

@@ -45,6 +45,8 @@ dev = [
     "notebook",
     "pandas",
     "pre-commit>=3.7.1",
+    "pytest",
+    "pytest-cov",
     "tox>=4",
 ]
 

src/opinf/__init__.py

@@ -7,7 +7,7 @@
     https://github.com/Willcox-Research-Group/rom-operator-inference-Python3
 """
 
-__version__ = "0.5.14"
+__version__ = "0.5.15"
 
 from . import (
     basis,

src/opinf/ddt/_base.py

@@ -104,7 +104,7 @@ def _check_dimensions(self, states, inputs, check_against_time=True):
             raise errors.DimensionalityError("states must be two-dimensional")
         if check_against_time and states.shape[-1] != self.time_domain.size:
             raise errors.DimensionalityError(
-                "states not aligned with time_domain"
+                "states and time_domain not aligned"
             )
         if inputs is not None:
             if inputs.ndim == 1:
@@ -135,6 +135,41 @@ def estimate(self, states, inputs=None):
         _inputs : (m, k') ndarray or None
             Inputs corresponding to ``_states``, if applicable.
             **Only returned** if ``inputs`` is provided.
+
+        Raises
+        ------
+        opinf.errors.DimensionalityError
+            If the ``states`` or ``inputs`` are not aligned with the
+            :attr:`time_domain`.
+        """
+        raise NotImplementedError  # pragma: no cover
+
+    @abc.abstractmethod
+    def mask(self, arr):
+        """Map an array from the training time domain to the domain of the
+        estimated time derivatives.
+
+        This method is used in post-hoc regularization selection routines.
+
+        Parameters
+        ----------
+        arr : (..., k) ndarray
+            Array (states, inputs, etc.) aligned with the training time domain.
+
+        Returns
+        -------
+        _arr : (..., k') ndarray
+            Array mapped to the domain of the estimated time derivatives.
+
+        Examples
+        --------
+        >>> Q, dQ = estimator.esimate(states)
+        >>> Q2 = estimator.mask(states)
+        >>> np.all(Q2 == Q)
+        True
+        >>> Q3 = estimator.mask(other_states_on_same_time_grid)
+        >>> Q3.shape == Q.shape
+        True
         """
         raise NotImplementedError  # pragma: no cover
 

src/opinf/ddt/_finite_difference.py

@@ -774,9 +774,8 @@ class UniformFiniteDifferencer(DerivativeEstimatorTemplate):
         Time domain corresponding to the snapshot data.
         This class requires uniformly spaced time domains, see
         :class:`NonuniformFiniteDifferencer` for non-uniform domains.
-    scheme : str or callable
-        Finite difference scheme to use.
-        **Options:**
+    scheme : str
+        Finite difference scheme to use. **Options:**
 
         * ``'fwd1'``: first-order forward differences, see :func:`fwd1`.
         * ``'fwd2'``: second-order forward differences, see :func:`fwd2`.
@@ -796,17 +795,6 @@ class UniformFiniteDifferencer(DerivativeEstimatorTemplate):
         * ``'ord2'``: second-order differences, see :func:`ord2`.
         * ``'ord4'``: fourth-order differences, see :func:`ord4`.
         * ``'ord6'``: sixth-order differences, see :func:`ord6`.
-
-        If ``scheme`` is a callable function, its signature must match the
-        following syntax.
-
-        .. code-block:: python
-
-           _states, ddts = scheme(states, dt)
-           _states, ddts, _inputs = scheme(states, dt, inputs)
-
-        Here ``dt`` is a positive float, the uniform time step.
-        Each output should have the same number of columns.
     """
 
     _schemes = types.MappingProxyType(
@@ -832,7 +820,7 @@ class UniformFiniteDifferencer(DerivativeEstimatorTemplate):
         }
     )
 
-    def __init__(self, time_domain, scheme="ord4"):
+    def __init__(self, time_domain, scheme: str = "ord4"):
         """Store the time domain and set the finite difference scheme."""
         DerivativeEstimatorTemplate.__init__(self, time_domain)
 
@@ -842,13 +830,9 @@ def __init__(self, time_domain, scheme="ord4"):
             raise ValueError("time domain must be uniformly spaced")
 
         # Set the finite difference scheme.
-        if not callable(scheme):
-            if scheme not in self._schemes:
-                raise ValueError(
-                    f"invalid finite difference scheme '{scheme}'"
-                )
-            scheme = self._schemes[scheme]
-        self.__scheme = scheme
+        if scheme not in self._schemes:
+            raise ValueError(f"invalid finite difference scheme '{scheme}'")
+        self.__scheme = self._schemes[scheme]
 
     # Properties --------------------------------------------------------------
     @property
@@ -899,6 +883,47 @@ def estimate(self, states, inputs=None):
         states, inputs = self._check_dimensions(states, inputs, False)
         return self.scheme(states, self.dt, inputs)
 
+    def mask(self, arr):
+        """Map an array from the training time domain to the domain of the
+        estimated time derivatives.
+
+        This method is used in post-hoc regularization selection routines.
+
+        Parameters
+        ----------
+        arr : (..., k) ndarray
+            Array (states, inputs, etc.) aligned with the training time domain.
+
+        Returns
+        -------
+        _arr : (..., k') ndarray
+            Array mapped to the domain of the estimated time derivatives.
+
+        Examples
+        --------
+        >>> Q, dQ = estimator.esimate(states)
+        >>> Q2 = estimator.mask(states)
+        >>> np.all(Q2 == Q)
+        True
+        >>> Q3 = estimator.mask(other_states_on_same_time_grid)
+        >>> Q3.shape == Q.shape
+        True
+        """
+        schema = self.scheme.__name__
+        mode, order = schema[:3], int(schema[3])
+        if mode == "ord":
+            return arr
+        elif mode == "fwd":
+            cols = slice(0, -order)
+        elif mode == "bwd":
+            cols = slice(order, None)
+        elif mode == "ctr":
+            margin = order // 2
+            cols = slice(margin, -margin)
+        else:  # pragma: no cover
+            raise RuntimeError("invalid scheme name!")
+        return arr[..., cols]
+
 
 class NonuniformFiniteDifferencer(DerivativeEstimatorTemplate):
     """Time derivative estimation with finite differences for state snapshots
@@ -967,6 +992,35 @@ def estimate(self, states, inputs=None):
             return states, ddts, inputs
         return states, ddts
 
+    def mask(self, arr):
+        """Map an array from the training time domain to the domain of the
+        estimated time derivatives. Since this class provides an estimate at
+        every time step, this method simply returns ``arr``.
+
+        This method is used in post-hoc regularization selection routines.
+
+        Parameters
+        ----------
+        arr : (..., k) ndarray
+            Array (states, inputs, etc.) aligned with the training time domain.
+
+        Returns
+        -------
+        _arr : (..., k') ndarray
+            Array mapped to the domain of the estimated time derivatives.
+
+        Examples
+        --------
+        >>> Q, dQ = estimator.esimate(states)
+        >>> Q2 = estimator.mask(states)
+        >>> np.all(Q2 == Q)
+        True
+        >>> Q3 = estimator.mask(other_states_on_same_time_grid)
+        >>> Q3.shape == Q.shape
+        True
+        """
+        return arr
+
 
 # Old API =====================================================================
 def ddt_uniform(states, dt, order=2):

src/opinf/ddt/_interpolation.py

@@ -189,3 +189,55 @@ def estimate(self, states, inputs=None):
         if inputs is None:
             return states, ddts
         return states, ddts, inputs
+
+    def mask(self, arr):
+        """Map an array from the training time domain to the domain of the
+        estimated time derivatives.
+
+
+        This method is used in post-hoc regularization selection routines.
+
+        Parameters
+        ----------
+        arr : (..., k) ndarray
+            Array (states, inputs, etc.) aligned with the training time domain.
+
+        Returns
+        -------
+        _arr : (..., k) ndarray
+            Array mapped to the domain of the estimated time derivatives.
+
+        Notes
+        -----
+        If :attr:`new_time_domain` is not the same as :attr:`time_domain`,
+        this method interpolates the ``arr`` and evaluates the interpolant
+        (not its derivative) over the :attr:`new_time_domain`.
+
+        Examples
+        --------
+        >>> Q, dQ = estimator.esimate(states)
+        >>> Q2 = estimator.mask(states)
+        >>> np.all(Q2 == Q)
+        True
+        >>> Q3 = estimator.mask(other_states_on_same_time_grid)
+        >>> Q3.shape == Q.shape
+        True
+        """
+        if self.new_time_domain is None:
+            return arr
+
+        # Interpolate to the new time domain.
+        statespline = self.InterpolatorClass(
+            self.time_domain,
+            arr,
+            **self.options,
+        )
+        return statespline(self.new_time_domain)
+
+    # Verification ------------------------------------------------------------
+    def verify(self, plot=False, return_errors=False):
+        new_time_domain = self.__t2
+        self.__t2 = None
+        out = super().verify(plot, return_errors)
+        self.__t2 = new_time_domain
+        return out