diff --git a/.githooks/pre-commit.sh b/.githooks/pre-commit.sh new file mode 100755 index 0000000..3e68739 --- /dev/null +++ b/.githooks/pre-commit.sh @@ -0,0 +1,83 @@ +#!/usr/bin/env bash + +# Git pre-commit hook to check staged Python files for formatting issues with +# yapf. +# +# INSTALLING: Copy this script into `.git/hooks/pre-commit`, and mark it as +# executable. +# +# This requires that yapf is installed and runnable in the environment running +# the pre-commit hook. +# +# When running, this first checks for unstaged changes to staged files, and if +# there are any, it will exit with an error. Files with unstaged changes will be +# printed. +# +# If all staged files have no unstaged changes, it will run yapf against them, +# leaving the formatting changes unstaged. Changed files will be printed. +# +# BUGS: This does not leave staged changes alone when used with the -a flag to +# git commit, due to the fact that git stages ALL unstaged files when that flag +# is used. + +# Find all staged Python files, and exit early if there aren't any. +PYTHON_FILES=() +while IFS=$'\n' read -r line; do PYTHON_FILES+=("$line"); done \ + < <(git diff --name-only --cached --diff-filter=AM | grep --color=never '.py$') +if [ ${#PYTHON_FILES[@]} -eq 0 ]; then + exit 0 +fi + +########## PIP VERSION ############# +# Verify that yapf is installed; if not, warn and exit. +if ! command -v yapf >/dev/null; then + echo 'yapf not on path; can not format. Please install yapf:' + echo ' pip install yapf' + exit 2 +fi +######### END PIP VERSION ########## + +########## PIPENV VERSION ########## +# if ! pipenv run yapf --version 2>/dev/null 2>&1; then +# echo 'yapf not on path; can not format. Please install yapf:' +# echo ' pipenv install yapf' +# exit 2 +# fi +###### END PIPENV VERSION ########## + + +# Check for unstaged changes to files in the index. +CHANGED_FILES=() +while IFS=$'\n' read -r line; do CHANGED_FILES+=("$line"); done \ + < <(git diff --name-only "${PYTHON_FILES[@]}") +if [ ${#CHANGED_FILES[@]} -gt 0 ]; then + echo 'You have unstaged changes to some files in your commit; skipping ' + echo 'auto-format. Please stage, stash, or revert these changes. You may ' + echo 'find `git stash -k` helpful here.' + echo 'Files with unstaged changes:' "${CHANGED_FILES[@]}" + exit 1 +fi + +# Format all staged files, then exit with an error code if any have uncommitted +# changes. +echo 'Formatting staged Python files . . .' + +########## PIP VERSION ############# +yapf -i -r "${PYTHON_FILES[@]}" +######### END PIP VERSION ########## + +########## PIPENV VERSION ########## +# pipenv run yapf -i -r "${PYTHON_FILES[@]}" +###### END PIPENV VERSION ########## + + +CHANGED_FILES=() +while IFS=$'\n' read -r line; do CHANGED_FILES+=("$line"); done \ + < <(git diff --name-only "${PYTHON_FILES[@]}") +if [ ${#CHANGED_FILES[@]} -gt 0 ]; then + echo 'Reformatted staged files. Please review and stage the changes.' + echo 'Files updated: ' "${CHANGED_FILES[@]}" + exit 1 +else + exit 0 +fi diff --git a/.github/workflows/check-formatting.yml b/.github/workflows/check-formatting.yml new file mode 100644 index 0000000..565a10e --- /dev/null +++ b/.github/workflows/check-formatting.yml @@ -0,0 +1,12 @@ +name: Check formatting +on: [push] +jobs: + formatting-check: + name: Formatting Check + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + - name: run YAPF to test if python code is correctly formatted + uses: AlexanderMelde/yapf-action@master + with: + args: --verbose diff --git a/.isort.cfg b/.isort.cfg new file mode 100644 index 0000000..46d7a41 --- /dev/null +++ b/.isort.cfg @@ -0,0 +1,3 @@ +[isort] +multi_line_output=3 +include_trailing_comma = true diff --git a/.mypy.ini b/.mypy.ini index 5de30ce..d4a7c52 100644 --- a/.mypy.ini +++ b/.mypy.ini @@ -1,29 +1,4 @@ [mypy] # https://mypy.readthedocs.io/en/stable/running_mypy.html#missing-imports -[mypy-sklearn.*] -ignore_missing_imports = True - -[mypy-abc] -ignore_missing_imports = True - -[mypy-pandas] -ignore_missing_imports = True - -[mypy-scipy] -ignore_missing_imports = True - -[mypy-scipy.*] -ignore_missing_imports = True - -[mypy-joblib] -ignore_missing_imports = True - -[mypy-picos] -ignore_missing_imports = True - -[mypy-optht] -ignore_missing_imports = True - -[mypy-pytest] ignore_missing_imports = True diff --git a/README.rst b/README.rst index 4c2eaca..2a2bbbc 100644 --- a/README.rst +++ b/README.rst @@ -130,6 +130,16 @@ The documentation can be compiled using $ cd doc $ make html +If you want a hook to check source code formatting before allowing a commit, +you can use + +.. code-block:: sh + + $ cd .git/hooks/ + $ ln -s ../../.githooks/pre-commit.sh . + $ chmod +x ./pre-commit.sh + +You will need ``yapf`` installed for this. Related packages ================ diff --git a/constraints.txt b/constraints.txt index a37a08c..742fc3e 100644 --- a/constraints.txt +++ b/constraints.txt @@ -4,6 +4,7 @@ scikit-learn>=1.0.1 PICOS>=2.2.52 pandas>=1.3.1 optht>=0.2.0 +Deprecated>=1.2.13 pytest matplotlib diff --git a/doc/conf.py b/doc/conf.py index 1fcdf3a..7e7a622 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -12,8 +12,8 @@ # import os import sys -sys.path.insert(0, os.path.abspath('../')) +sys.path.insert(0, os.path.abspath('../')) # -- Project information ----------------------------------------------------- @@ -21,7 +21,6 @@ copyright = '2021, Steven Dahdah' author = 'Steven Dahdah' - # -- General configuration --------------------------------------------------- # Add any Sphinx extension module names here, as strings. They can be @@ -49,7 +48,6 @@ # This pattern also affects html_static_path and html_extra_path. exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] - # -- Options for HTML output ------------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for diff --git a/doc/index.rst b/doc/index.rst index faa6207..a705d7e 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -136,6 +136,16 @@ The documentation can be compiled using $ cd doc $ make html +If you want a hook to check source code formatting before allowing a commit, +you can use + +.. code-block:: sh + + $ cd .git/hooks/ + $ ln -s ../../.githooks/pre-commit.sh . + $ chmod +x ./pre-commit.sh + +You will need ``yapf`` installed for this. Related packages ================ @@ -182,7 +192,10 @@ References .. [bilinear] Daniel Bruder, Xun Fu, and Ram Vasudevan. "Advantages of bilinear Koopman realizations for the modeling and control of systems with unknown dynamics." arXiv:2010.09961v3 [cs.RO] (2020). - https://arxiv.org/abs/2010.09961v3 +.. [local] Giorgos Mamakoukas, Ian Abraham, and Todd D. Murphey. "Learning + Stable Models for Prediction and Control." arXiv:2005.04291v2 [cs.RO] + (2022). https://arxiv.org/abs/2005.04291v2 https://arxiv.org/abs/2010.09961v3 +.. [scorers] https://scikit-learn.org/stable/modules/model_evaluation.html Citation ======== diff --git a/examples/example_pipeline_vdp.py b/examples/example_pipeline_vdp.py new file mode 100644 index 0000000..6c24f31 --- /dev/null +++ b/examples/example_pipeline_vdp.py @@ -0,0 +1,142 @@ +"""Example of how to use the Koopman pipeline.""" + +import numpy as np +from matplotlib import pyplot as plt + +import pykoop +import pykoop.dynamic_models + + +def main() -> None: + """Demonstrate how to use the Koopman pipeline.""" + # Get sample data + X_vdp = pykoop.example_data_vdp() + + # Create pipeline + kp = pykoop.KoopmanPipeline( + lifting_functions=[( + 'sp', + pykoop.SplitPipeline( + lifting_functions_state=[ + ('pl', pykoop.PolynomialLiftingFn(order=3)) + ], + lifting_functions_input=None, + ), + )], + regressor=pykoop.Edmd(), + ) + + # Take last episode for validation + X_train = X_vdp[X_vdp[:, 0] < 4] + X_valid = X_vdp[X_vdp[:, 0] == 4] + + # Fit the pipeline + kp.fit(X_train, n_inputs=1, episode_feature=True) + + # Extract initial conditions and input from validation episode + x0 = X_valid[[0], 1:3] + u = X_valid[:, 3:] + + # Predict with re-lifting between timesteps (default) + X_pred_local = kp.predict_state( + x0, + u, + relift_state=True, + episode_feature=False, + ) + + # Predict without re-lifting between timesteps + X_pred_global = kp.predict_state( + x0, + u, + relift_state=False, + episode_feature=False, + ) + + # Plot trajectories in phase space + fig, ax = plt.subplots(constrained_layout=True) + ax.plot( + X_valid[:, 1], + X_valid[:, 2], + label='True trajectory', + ) + ax.plot( + X_pred_local[:, 0], + X_pred_local[:, 1], + '--', + label='Local prediction', + ) + ax.plot( + X_pred_global[:, 0], + X_pred_global[:, 1], + '--', + label='Global prediction', + ) + ax.set_xlabel('$x_1[k]$') + ax.set_ylabel('$x_2[k]$') + ax.legend() + ax.grid(linestyle='--') + + # Lift validation set + Psi_valid = kp.lift(X_valid[:, 1:], episode_feature=False) + + # Predict lifted state with re-lifting between timesteps (default) + Psi_pred_local = kp.predict_state( + x0, + u, + relift_state=True, + return_lifted=True, + return_input=True, + episode_feature=False, + ) + + # Predict lifted state without re-lifting between timesteps + Psi_pred_global = kp.predict_state( + x0, + u, + relift_state=False, + return_lifted=True, + return_input=True, + episode_feature=False, + ) + + fig, ax = plt.subplots( + kp.n_states_out_, + 1, + constrained_layout=True, + sharex=True, + squeeze=False, + ) + for i in range(ax.shape[0]): + ax[i, 0].plot(Psi_valid[:, i], label='True trajectory') + ax[i, 0].plot(Psi_pred_local[:, i], '--', label='Local prediction') + ax[i, 0].plot(Psi_pred_global[:, i], '--', label='Global prediction') + ax[i, 0].grid(linestyle='--') + ax[i, 0].set_ylabel(rf'$\vartheta_{i + 1}[k]$') + + ax[-1, 0].set_xlabel('$k$') + ax[0, 0].legend() + + fig, ax = plt.subplots( + kp.n_inputs_out_, + 1, + constrained_layout=True, + sharex=True, + squeeze=False, + ) + for i in range(ax.shape[0]): + j = kp.n_states_out_ + i + ax[i, 0].plot(Psi_valid[:, j], label='True trajectory') + ax[i, 0].plot(Psi_pred_local[:, j], '--', label='Local prediction') + ax[i, 0].plot(Psi_pred_global[:, j], '--', label='Global prediction') + ax[i, 0].grid(linestyle='--') + + ax[-1, 0].set_xlabel('$k$') + ax[0, 0].legend() + ax[0, 0].set_ylabel('$u[k]$') + + plt.show() + + +if __name__ == '__main__': + main() diff --git a/pykoop/__init__.py b/pykoop/__init__.py index 215c7a0..543ba65 100644 --- a/pykoop/__init__.py +++ b/pykoop/__init__.py @@ -1,13 +1,32 @@ """Koopman operator identification library in Python.""" -from .koopman_pipeline import (EpisodeDependentLiftingFn, - EpisodeIndependentLiftingFn, KoopmanLiftingFn, - KoopmanPipeline, KoopmanRegressor, - SplitPipeline, combine_episodes, shift_episodes, - split_episodes, strip_initial_conditions) -from .lifting_functions import (BilinearInputLiftingFn, DelayLiftingFn, - PolynomialLiftingFn, SkLearnLiftingFn) +from .koopman_pipeline import ( + EpisodeDependentLiftingFn, + EpisodeIndependentLiftingFn, + KoopmanLiftingFn, + KoopmanPipeline, + KoopmanRegressor, + SplitPipeline, + combine_episodes, + extract_initial_conditions, + extract_input, + score_state, + shift_episodes, + split_episodes, + strip_initial_conditions, +) +from .lifting_functions import ( + BilinearInputLiftingFn, + DelayLiftingFn, + PolynomialLiftingFn, + SkLearnLiftingFn, +) from .regressors import Dmd, Dmdc, Edmd from .tsvd import Tsvd -from .util import (AnglePreprocessor, example_data_msd, random_input, - random_state) +from .util import ( + AnglePreprocessor, + example_data_msd, + example_data_vdp, + random_input, + random_state, +) diff --git a/pykoop/koopman_pipeline.py b/pykoop/koopman_pipeline.py index 0787469..5bdc151 100644 --- a/pykoop/koopman_pipeline.py +++ b/pykoop/koopman_pipeline.py @@ -92,15 +92,21 @@ """ import abc +import logging from typing import Any, Callable, Dict, List, Optional, Tuple import numpy as np import pandas import sklearn.base import sklearn.metrics +from deprecated import deprecated from ._sklearn_metaestimators import metaestimators +# Create logger +log = logging.getLogger(__name__) +log.addHandler(logging.NullHandler()) + class _LiftRetractMixin(metaclass=abc.ABCMeta): """Mixin providing more convenient lift/retract functions. @@ -1395,8 +1401,7 @@ def set_params(self, **kwargs) -> 'SplitPipeline': class KoopmanPipeline(metaestimators._BaseComposition, sklearn.base.BaseEstimator, - sklearn.base.TransformerMixin, - _LiftRetractMixin): + sklearn.base.TransformerMixin, _LiftRetractMixin): """Meta-estimator for chaining lifting functions with an estimator. Attributes @@ -1535,9 +1540,11 @@ def fit(self, ValueError If constructor or fit parameters are incorrect. """ - X = sklearn.utils.validation.check_array(X, - ensure_min_samples=2, - **self._check_array_params) + X = sklearn.utils.validation.check_array( + X, + ensure_min_samples=2, + **self._check_array_params, + ) if self.regressor is None: raise ValueError( '`regressor` must be specified in order to use `fit()`.') @@ -1754,6 +1761,7 @@ def score(self, X: np.ndarray, y: np.ndarray = None) -> float: score = scorer(self, X, None) return score + @deprecated('Use `predict_state` instead') def predict_multistep(self, X: np.ndarray) -> np.ndarray: """Perform a multi-step prediction for the first state of each episode. @@ -1829,12 +1837,187 @@ def predict_multistep(self, X: np.ndarray) -> np.ndarray: episode_feature=self.episode_feature_) return X_p + def predict_state( + self, + X_initial: np.ndarray, + U: np.ndarray, + relift_state: bool = True, + return_lifted: bool = False, + return_input: bool = False, + episode_feature: bool = None, + ) -> np.ndarray: + """Predict state given input for each episode. + + Parameters + ---------- + X_initial : np.ndarray + Initial state. + U : np.ndarray + Input. Length of prediction is governed by length of input. + relift_state : bool + If true, retract and re-lift state between prediction steps + (default). Otherwise, only retract the state after all predictions + are made. Correspond to the local and global error definitions of + [local]_. + return_lifted : bool + If true, return the lifted state. If false, return the original + state (default). + return_input : bool + If true, return the input as well as the state. If false, return + only the original state (default). + episode_feature : bool + True if first feature indicates which episode a timestep is from. + If ``None``, ``self.episode_feature_`` is used. + + Returns + ------- + np.ndarray + Predicted state. If ``return_input``, input is appended to the + array. If ``return_lifted``, the predicted state (and input) are + returned in the lifted space. + + Raises + ------ + ValueError + If an episode is shorter than ``min_samples_``. + """ + # Check fit + sklearn.utils.validation.check_is_fitted(self, 'regressor_fit_') + # Set episode feature if unspecified + if episode_feature is None: + episode_feature = self.episode_feature_ + # Get Koopman ``A`` and ``B`` matrices + koop_mat = self.regressor_.coef_.T + A = koop_mat[:, :koop_mat.shape[0]] + B = koop_mat[:, koop_mat.shape[0]:] + # Split episodes + ep_X0 = split_episodes(X_initial, episode_feature=episode_feature) + ep_U = split_episodes(U, episode_feature=episode_feature) + episodes = [(ex[0], ex[1], eu[1]) for (ex, eu) in zip(ep_X0, ep_U)] + # Predict for each episode + predictions: List[Tuple[float, np.ndarray]] = [] + for (i, X0_i, U_i) in episodes: + # Check length of episode. + if X0_i.shape[0] != self.min_samples_: + raise ValueError(f'Initial condition in episode {i} has ' + f'{X0_i.shape[0]} samples but `min_samples_`=' + f'{self.min_samples_} samples are required.') + if U_i.shape[0] < self.min_samples_: + raise ValueError(f'Input in episode {i} has {U_i.shape[0]} ' + 'samples but at least `min_samples_`=' + f'{self.min_samples_} samples are required.') + crash_index = None + # Iterate over episode and make predictions + if relift_state: + # Number of steps in episode + n_steps_i = U_i.shape[0] + # Initial conditions + X_i = np.zeros((n_steps_i, self.n_states_in_)) + X_i[:self.min_samples_, :] = X0_i + for k in range(self.min_samples_, n_steps_i): + try: + # Lift state and input + window = np.s_[(k - self.min_samples_):k] + Theta_ikm1 = self.lift_state( + X_i[window, :], + episode_feature=False, + ) + Upsilon_ikm1 = self.lift_input( + np.hstack((X_i[window, :], U_i[window, :])), + episode_feature=False, + ) + # Predict + Theta_ik = Theta_ikm1 @ A.T + Upsilon_ikm1 @ B.T + # Retract. If more than one sample is returned by + # ``retract_state``, take only the last one. This will + # happen if there's a delay lifting function. + X_i[[k], :] = self.retract_state( + Theta_ik, + episode_feature=False, + )[[-1], :] + except ValueError as ve: + if (np.all(np.isfinite(Theta_ikm1)) + and np.all(np.isfinite(X_i)) + and np.all(np.isfinite(U_i)) + and np.all(np.isfinite(Upsilon_ikm1)) + and np.all(np.isfinite(Theta_ik))): + raise ve + else: + crash_index = k - 1 + X_i[crash_index:, :] = 0 + break + Theta_i = self.lift_state(X_i, episode_feature=False) + Upsilon_i = self.lift_input( + np.hstack((X_i, U_i)), + episode_feature=False, + ) + else: + # Number of steps in episode + n_steps_i = U_i.shape[0] - self.min_samples_ + # Initial conditions + Theta_i = np.zeros((n_steps_i, self.n_states_out_)) + Upsilon_i = np.zeros((n_steps_i, self.n_inputs_out_)) + Theta_i[[0], :] = self.lift_state(X0_i, episode_feature=False) + for k in range(1, n_steps_i + 1): + try: + X_ikm1 = self.retract_state( + Theta_i[[k - 1], :], + episode_feature=False, + ) + window = np.s_[k:(k + self.min_samples_)] + Upsilon_i[[k - 1], :] = self.lift_input( + np.hstack((X_ikm1, U_i[window, :])), + episode_feature=False, + ) + # Predict + if k < n_steps_i: + Theta_i[[k], :] = (Theta_i[[k - 1], :] @ A.T + + Upsilon_i[[k - 1], :] @ B.T) + except ValueError as ve: + if (np.all(np.isfinite(X_ikm1)) + and np.all(np.isfinite(Theta_i)) + and np.all(np.isfinite(Upsilon_i)) + and np.all(np.isfinite(U_i))): + raise ve + else: + crash_index = k - 1 + Theta_i[crash_index:, :] = 0 + Upsilon_i[crash_index:, :] = 0 + break + X_i = self.retract_state(Theta_i, episode_feature=False) + # If prediction crashed, set remaining entries to NaN + if crash_index is not None: + log.warning(f'Prediction diverged at index {crash_index}. ' + 'Remaining entries set to `NaN`.') + # Don't set ``U_i`` to NaN since it's a known input + X_i[crash_index:, :] = np.nan + Theta_i[crash_index:, :] = np.nan + Upsilon_i[crash_index:, :] = np.nan + # Choose what to return + if return_lifted: + if return_input: + predictions.append((i, np.hstack((Theta_i, Upsilon_i)))) + else: + predictions.append((i, Theta_i)) + else: + if return_input: + predictions.append((i, np.hstack((X_i, U_i)))) + else: + predictions.append((i, X_i)) + # Combine episodes + combined_episodes = combine_episodes( + predictions, + episode_feature=episode_feature, + ) + return combined_episodes + @staticmethod def make_scorer( n_steps: int = None, discount_factor: float = 1, regression_metric: str = 'neg_mean_squared_error', - multistep: bool = True + multistep: bool = True, + relift_state: bool = True, ) -> Callable[['KoopmanPipeline', np.ndarray, Optional[np.ndarray]], float]: """Make a Koopman pipeline scorer. @@ -1864,11 +2047,16 @@ def make_scorer( 'neg_mean_squared_log_error', 'neg_median_absolute_error', 'r2', 'neg_mean_absolute_percentage_error']. See [scorers]_. multistep : bool - If true, predict using :func:`predict_multistep`. Otherwise, + If true, predict using :func:`predict_state`. Otherwise, predict using :func:`predict` (one-step-ahead prediction). Multistep prediction is highly recommended unless debugging. If one-step-ahead prediciton is used, `n_steps` and `discount_factor` are ignored. + relift_state : bool + If true, retract and re-lift state between prediction steps + (default). Otherwise, only retract the state after all predictions + are made. Correspond to the local and global error definitions of + [local]_. Ignored if ``multistep`` is false. Returns ------- @@ -1879,85 +2067,72 @@ def make_scorer( ------ ValueError If ``discount_factor`` is negative or greater than one. - - References - ---------- - .. [scorers] https://scikit-learn.org/stable/modules/model_evaluation.html # noqa: E501 """ - # Check discount factor - if (discount_factor < 0) or (discount_factor > 1): - raise ValueError('`discount_factor` must be positive and less ' - 'than one.') - - # Valid ``regression_metric`` values: - regression_metrics = { - 'explained_variance': - sklearn.metrics.explained_variance_score, - 'r2': - sklearn.metrics.r2_score, - 'neg_mean_absolute_error': - sklearn.metrics.mean_absolute_error, - 'neg_mean_squared_error': - sklearn.metrics.mean_squared_error, - 'neg_mean_squared_log_error': - sklearn.metrics.mean_squared_log_error, - 'neg_median_absolute_error': - sklearn.metrics.median_absolute_error, - 'neg_mean_absolute_percentage_error': - sklearn.metrics.mean_absolute_percentage_error, - } - # Scores that do not need inversion - greater_is_better = ['explained_variance', 'r2'] - def koopman_pipeline_scorer(estimator: KoopmanPipeline, - X: np.ndarray, - y: np.ndarray = None) -> float: + def koopman_pipeline_scorer( + estimator: KoopmanPipeline, + X: np.ndarray, + y: np.ndarray = None, + ) -> float: # Shift episodes X_unshifted, X_shifted = shift_episodes( X, n_inputs=estimator.n_inputs_in_, - episode_feature=estimator.episode_feature_) + episode_feature=estimator.episode_feature_, + ) # Predict if multistep: - X_predicted = estimator.predict_multistep(X_unshifted) + # Get initial conditions for each episode + x0 = extract_initial_conditions( + X_unshifted, + min_samples=estimator.min_samples_, + n_inputs=estimator.n_inputs_in_, + episode_feature=estimator.episode_feature_, + ) + # Get inputs for each episode + u = extract_input( + X_unshifted, + n_inputs=estimator.n_inputs_in_, + episode_feature=estimator.episode_feature_, + ) + # Predict state for each episode + X_predicted = estimator.predict_state( + x0, + u, + relift_state=relift_state, + ) + # Score prediction + score = score_state( + X_predicted, + X_shifted, + n_steps=n_steps, + discount_factor=discount_factor, + regression_metric=regression_metric, + min_samples=estimator.min_samples_, + episode_feature=estimator.episode_feature_, + ) else: + # Warn about ignored non-default arguments + if not relift_state: + log.info('Ignoring `relift_state` since `multistep` is ' + 'false.') + if n_steps is not None: + log.info('Ignoring `n_steps` since `multistep` is false.') + if discount_factor != 1: + log.info('Ignoring `discount_factor` since `multistep` is ' + 'false.') + # Perform single-step prediction X_predicted = estimator.predict(X_unshifted) - # Strip episode feature and initial conditions - X_shifted = strip_initial_conditions(X_shifted, - estimator.min_samples_, - estimator.episode_feature_) - X_predicted = strip_initial_conditions(X_predicted, - estimator.min_samples_, - estimator.episode_feature_) - # Strip episode feature if present - if estimator.episode_feature_: - X_shifted = X_shifted[:, 1:] - X_predicted = X_predicted[:, 1:] - # Compute weights - weights: Optional[np.ndarray] - if multistep: - # Compute number of weights needed - n_samples = X_shifted.shape[0] - if (n_steps is None) or (n_steps > n_samples): - n_weights = n_samples - else: - n_weights = n_steps - # Compute weights. Weights after ``n_steps`` are 0. - weights = np.array( - [discount_factor**k for k in range(n_weights)] - + [0] * (n_samples - n_weights)) - else: - weights = None - # Calculate score - score = regression_metrics[regression_metric]( - X_shifted, - X_predicted, - sample_weight=weights, - multioutput='uniform_average', - ) - # Invert losses - if regression_metric not in greater_is_better: - score *= -1 + # Score prediction + score = score_state( + X_predicted, + X_shifted, + n_steps=None, + discount_factor=1, + regression_metric=regression_metric, + min_samples=estimator.min_samples_, + episode_feature=estimator.episode_feature_, + ) return score return koopman_pipeline_scorer @@ -1972,6 +2147,174 @@ def set_params(self, **kwargs) -> 'KoopmanPipeline': return self +def score_state( + X_predicted: np.ndarray, + X_expected: np.ndarray, + n_steps: int = None, + discount_factor: float = 1, + regression_metric: str = 'neg_mean_squared_error', + min_samples: int = 1, + episode_feature: bool = False, +) -> float: + """Score a predicted data matrix compared to an expected data matrix. + + Parameters + ---------- + X_predicted : np.ndarray + Predicted state data matrix. + X_expected : np.ndarray + Expected state data matrix. + n_steps : int + Number of steps ahead to predict. If ``None`` or longer than the + episode, will score the entire episode. + discount_factor : float + Discount factor used to weight the error timeseries. Should be + positive, with magnitude 1 or slightly less. The error at each + timestep is weighted by ``discount_factor**k``, where ``k`` is the + timestep. + regression_metric : str + Regression metric to use. One of ['explained_variance', + 'neg_mean_absolute_error', 'neg_mean_squared_error', + 'neg_mean_squared_log_error', 'neg_median_absolute_error', 'r2', + 'neg_mean_absolute_percentage_error']. See [scorers]_. + min_samples : int + Number of samples in initial condition. + episode_feature : bool + True if first feature indicates which episode a timestep is from. + + Returns + ------- + float + Score (greater is better). + """ + # Valid ``regression_metric`` values: + regression_metrics = { + 'explained_variance': + sklearn.metrics.explained_variance_score, + 'r2': + sklearn.metrics.r2_score, + 'neg_mean_absolute_error': + sklearn.metrics.mean_absolute_error, + 'neg_mean_squared_error': + sklearn.metrics.mean_squared_error, + 'neg_mean_squared_log_error': + sklearn.metrics.mean_squared_log_error, + 'neg_median_absolute_error': + sklearn.metrics.median_absolute_error, + 'neg_mean_absolute_percentage_error': + sklearn.metrics.mean_absolute_percentage_error, + } + # Scores that do not need inversion + greater_is_better = ['explained_variance', 'r2'] + # Strip episode feature and initial conditions + X_expected = strip_initial_conditions( + X_expected, + min_samples=min_samples, + episode_feature=episode_feature, + ) + X_predicted = strip_initial_conditions( + X_predicted, + min_samples=min_samples, + episode_feature=episode_feature, + ) + # Compute weights + weights = _weights_from_data_matrix( + X_expected, + n_steps=n_steps, + discount_factor=discount_factor, + episode_feature=episode_feature, + ) + # Strip episode feature if present + if episode_feature: + X_expected = X_expected[:, 1:] + X_predicted = X_predicted[:, 1:] + # Calculate score + score = regression_metrics[regression_metric]( + X_expected, + X_predicted, + sample_weight=weights, + multioutput='uniform_average', + ) + # Invert losses + if regression_metric not in greater_is_better: + score *= -1 + return score + + +def extract_initial_conditions( + X: np.ndarray, + min_samples: int = 1, + n_inputs: int = 0, + episode_feature: bool = False, +) -> np.ndarray: + """Extract initial conditions from each episode. + + Parameters + ---------- + X : np.ndarray + Data matrix. + min_samples : int + Number of samples in initial condition. + n_inputs : int + Number of input features at the end of ``X``. + episode_feature : bool + True if first feature indicates which episode a timestep is from. + + Returns + ------- + np.ndarray + Initial conditions from each episode. + """ + episodes = split_episodes(X, episode_feature=episode_feature) + # Strip each episode + initial_conditions = [] + for (i, X_i) in episodes: + if n_inputs == 0: + initial_condition = X_i[:min_samples, :] + else: + initial_condition = X_i[:min_samples, :-n_inputs] + initial_conditions.append((i, initial_condition)) + # Concatenate the initial conditions + X0 = combine_episodes(initial_conditions, episode_feature=episode_feature) + return X0 + + +def extract_input( + X: np.ndarray, + n_inputs: int = 0, + episode_feature: bool = False, +) -> np.ndarray: + """Extract input from a data matrix. + + Parameters + ---------- + X : np.ndarray + Data matrix. + n_inputs : int + Number of input features at the end of ``X``. + episode_feature : bool + True if first feature indicates which episode a timestep is from. + + Returns + ------- + np.ndarray + Input extracted from data matrix. + """ + episodes = split_episodes(X, episode_feature=episode_feature) + # Strip each episode + inputs = [] + for (i, X_i) in episodes: + if n_inputs == 0: + input_ = np.zeros((X_i.shape[0], 0)) + else: + n_states = X_i.shape[1] - n_inputs + input_ = X_i[:, n_states:] + inputs.append((i, input_)) + # Concatenate the inputs + u = combine_episodes(inputs, episode_feature=episode_feature) + return u + + def strip_initial_conditions(X: np.ndarray, min_samples: int = 1, episode_feature: bool = False) -> np.ndarray: @@ -1985,6 +2328,11 @@ def strip_initial_conditions(X: np.ndarray, Number of samples in initial condition. episode_feature : bool True if first feature indicates which episode a timestep is from. + + Returns + ------- + np.ndarray + Data matrix with initial conditions removed. """ episodes = split_episodes(X, episode_feature=episode_feature) # Strip each episode @@ -2117,3 +2465,58 @@ def combine_episodes(episodes: List[Tuple[float, np.ndarray]], # Concatenate the combined episodes Xc = np.vstack(combined_episodes) return Xc + + +def _weights_from_data_matrix( + X: np.ndarray, + n_steps: int = None, + discount_factor: float = 1, + episode_feature: bool = False, +) -> np.ndarray: + """Create an array of scoring weights from a data matrix. + + Parameters + ---------- + X : np.ndarray + Data matrix + n_steps : int + Number of steps ahead to predict. If ``None`` or longer than the + episode, will weight the entire episode. + discount_factor : float + Discount factor used to weight the error timeseries. Should be + positive, with magnitude 1 or slightly less. The error at each + timestep is weighted by ``discount_factor**k``, where ``k`` is the + timestep. + episode_feature : bool + True if first feature indicates which episode a timestep is from. + + Returns + ------- + np.ndarray + Array of weights use for scoring. + + Raises + ------ + ValueError + If ``discount_factor`` is not in [0, 1]. + """ + # Check discount factor + if (discount_factor < 0) or (discount_factor > 1): + raise ValueError('`discount_factor` must be positive and less ' + 'than one.') + weights_list = [] + episodes = split_episodes(X, episode_feature=episode_feature) + for i, X_i in episodes: + # Compute number of nonzero weights needed + n_samples_i = X_i.shape[0] + if n_steps is None: + n_nonzero_weights_i = n_samples_i + else: + n_nonzero_weights_i = min(n_steps, n_samples_i) + # Compute weights. Weights after ``n_steps`` are 0. + weights_i = np.array( + [discount_factor**k for k in range(n_nonzero_weights_i)] + + [0] * (n_samples_i - n_nonzero_weights_i)) + weights_list.append(weights_i) + weights = np.concatenate(weights_list) + return weights diff --git a/pykoop/lmi_regressors.py b/pykoop/lmi_regressors.py index 111ea89..8f608e6 100644 --- a/pykoop/lmi_regressors.py +++ b/pykoop/lmi_regressors.py @@ -25,6 +25,7 @@ # Create logger log = logging.getLogger(__name__) +log.addHandler(logging.NullHandler()) # Create temporary cache directory for memoized computations _cachedir = tempfile.TemporaryDirectory(prefix='pykoop_') @@ -339,16 +340,18 @@ def _create_base_problem( # Choose method to handle inverse of H if inv_method == 'inv': H_inv = picos.Constant('H^-1', _calc_Hinv(H)) - problem.add_constraint(picos.block([ - [Z, U], - [U.T, H_inv], - ]) >> picos_eps) + problem.add_constraint( + picos.block([ + [Z, U], + [U.T, H_inv], + ]) >> picos_eps) elif inv_method == 'pinv': H_inv = picos.Constant('H^+', _calc_Hpinv(H)) - problem.add_constraint(picos.block([ - [Z, U], - [U.T, H_inv], - ]) >> picos_eps) + problem.add_constraint( + picos.block([ + [Z, U], + [U.T, H_inv], + ]) >> picos_eps) elif inv_method == 'eig': VsqrtLmb = picos.Constant('(V Lambda^(1/2))', _calc_VsqrtLmb(H)) problem.add_constraint( diff --git a/pykoop/util.py b/pykoop/util.py index 9eea1ca..42a6556 100644 --- a/pykoop/util.py +++ b/pykoop/util.py @@ -289,3 +289,46 @@ def example_data_msd() -> np.ndarray: # Stack data and return X_msd = np.vstack(X_msd_lst) return X_msd + + +def example_data_vdp() -> np.ndarray: + """Get example Van der Pol oscillator data. + + Returns + ------- + np.ndarray + Sample Van der Pol oscillator data. + """ + # Create Van der Pol object + t_step = 0.01 + vdp = dynamic_models.DiscreteVanDerPol(t_step, mu=2) + # Initial conditions and inputs + t_range = (0, 10) + t = np.arange(*t_range, t_step) + conditions = [ + (0, np.array([1, 0]), 0.1 * np.sin(t)), + (1, np.array([0, 1]), 0.2 * np.cos(t)), + (2, np.array([-1, 0]), -0.2 * np.sin(t)), + (3, np.array([0, -1]), -0.1 * np.cos(t)), + (4, np.array([0.5, 0]), 0.1 * np.cos(t)), + ] + X_vdp_lst = [] + # Loop over initial conditions and inputs + for ep, x0, u in conditions: + # Simulate ODE + t, x = vdp.simulate( + t_range=t_range, + t_step=t_step, + x0=vdp.x0(x0), + u=u, + ) + # Format the data + X_vdp_lst.append( + np.hstack(( + ep * np.ones((t.shape[0], 1)), + x, + np.reshape(u, (-1, 1)), + ))) + # Stack data and return + X_vdp = np.vstack(X_vdp_lst) + return X_vdp diff --git a/requirements.txt b/requirements.txt index 8ca26a3..d374dff 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,6 +4,7 @@ scikit-learn PICOS pandas optht +Deprecated pytest matplotlib diff --git a/setup.py b/setup.py index ca00ec9..4a00c8a 100644 --- a/setup.py +++ b/setup.py @@ -33,6 +33,7 @@ 'picos>=2.2.52', 'pandas>=1.3.1', 'optht>=0.2.0', + 'Deprecated>=1.2.13', ], extra_require={ 'MOSEK solver': ['mosek>=9.2.49'], diff --git a/tests/test_koopman_pipeline.py b/tests/test_koopman_pipeline.py index a67224c..221e1c6 100644 --- a/tests/test_koopman_pipeline.py +++ b/tests/test_koopman_pipeline.py @@ -4,6 +4,7 @@ import pykoop import pykoop.dynamic_models +import pykoop.koopman_pipeline def test_kp_transform_no_lf(): @@ -322,6 +323,411 @@ def test_combine_episodes(X, episodes, episode_feature): np.testing.assert_allclose(X_actual, X) +@pytest.mark.parametrize( + 'params', + [ + { + 'X_predicted': np.array([ + [1, 2, 3, 4], + [2, 3, 3, 2], + ]).T, + 'X_expected': np.array([ + [1, 2, 3, 4], + [2, 3, 3, 2], + ]).T, + 'n_steps': None, + 'discount_factor': 1, + 'regression_metric': 'neg_mean_squared_error', + 'min_samples': 1, + 'episode_feature': False, + 'score_exp': 0, + }, + { + 'X_predicted': np.array([ + [1, 2], + [2, 3], + ]).T, + 'X_expected': np.array([ + [1, 4], + [2, 2], + ]).T, + 'n_steps': None, + 'discount_factor': 1, + 'regression_metric': 'neg_mean_squared_error', + 'min_samples': 1, + 'episode_feature': False, + 'score_exp': -np.mean([2**2, 1]), + }, + { + 'X_predicted': np.array([ + [1, 2, 3, 5], + ]).T, + 'X_expected': np.array([ + [1, 2, 3, 3], + ]).T, + 'n_steps': None, + 'discount_factor': 1, + 'regression_metric': 'neg_mean_squared_error', + 'multistep': True, + 'min_samples': 1, + 'episode_feature': False, + 'score_exp': -np.mean([0, 0, 2**2]), + }, + { + 'X_predicted': np.array([ + [1, 2, 3, 5], + ]).T, + 'X_expected': np.array([ + [1, 2, 3, 3], + ]).T, + 'n_steps': None, + 'discount_factor': 1, + 'regression_metric': 'neg_mean_absolute_error', + 'multistep': True, + 'min_samples': 1, + 'episode_feature': False, + 'score_exp': -np.mean([0, 0, 2]), + }, + { + 'X_predicted': np.array([ + [1, 2, 3, 5], + ]).T, + 'X_expected': np.array([ + [1, 2, 3, 4], + ]).T, + 'n_steps': None, + 'discount_factor': 1, + 'regression_metric': 'neg_mean_squared_error', + 'min_samples': 2, + 'episode_feature': False, + 'score_exp': -np.mean([0, 1]), + }, + { + 'X_predicted': np.array([ + [1, 2, 3, 5], + ]).T, + 'X_expected': np.array([ + [1, 2, 3, 4], + ]).T, + 'n_steps': 2, + 'discount_factor': 1, + 'regression_metric': 'neg_mean_squared_error', + 'min_samples': 1, + 'episode_feature': False, + 'score_exp': 0, + }, + { + 'X_predicted': np.array([ + [1, 2, 3, 5], + ]).T, + 'X_expected': np.array([ + [1, 2, 3, 4], + ]).T, + 'n_steps': None, + 'discount_factor': 0.5, + 'regression_metric': 'neg_mean_squared_error', + 'min_samples': 1, + 'episode_feature': False, + # (0 * 1 + 0 * 0.5 + 1 * 0.25) / (1 + 0.5 + 0.25) + 'score_exp': -0.25 / 1.75, + }, + { + 'X_predicted': np.array([ + [0, 0, 0, 1, 1, 1], + [1, 2, 3, 4, 5, 6], + ]).T, + 'X_expected': np.array([ + [0, 0, 0, 1, 1, 1], + [1, 2, 4, 4, 6, 6], + ]).T, + 'n_steps': None, + 'discount_factor': 1, + 'regression_metric': 'neg_mean_squared_error', + 'min_samples': 1, + 'episode_feature': True, + 'score_exp': -np.mean([0, 1, 1, 0]), + }, + { + 'X_predicted': np.array([ + [0, 0, 0, 1, 1, 1], + [1, 2, 3, 4, 5, 6], + ]).T, + 'X_expected': np.array([ + [0, 0, 0, 1, 1, 1], + [1, 2, 4, 4, 6, 6], + ]).T, + 'n_steps': 1, + 'discount_factor': 1, + 'regression_metric': 'neg_mean_squared_error', + 'min_samples': 1, + 'episode_feature': True, + 'score_exp': -np.mean([0, 1]), + }, + { + 'X_predicted': np.array([ + [0, 0, 0, 1, 1, 1], + [1, 2, 3, 4, 5, 6], + ]).T, + 'X_expected': np.array([ + [0, 0, 0, 1, 1, 1], + [1, 2, 4, 4, 6, 6], + ]).T, + 'n_steps': None, + 'discount_factor': 0.5, + 'regression_metric': 'neg_mean_squared_error', + 'min_samples': 1, + 'episode_feature': True, + 'score_exp': -(0.5 + 1) / (1 + 0.5 + 1 + 0.5), + }, + { + 'X_predicted': np.array([ + [0, 0, 0, 1, 1, 1], + [1, 2, 3, 4, 5, 6], + ]).T, + 'X_expected': np.array([ + [0, 0, 0, 1, 1, 1], + [1, 2, 4, 4, 6, 6], + ]).T, + 'n_steps': 1, + 'discount_factor': 0.5, + 'regression_metric': 'neg_mean_squared_error', + 'min_samples': 1, + 'episode_feature': True, + 'score_exp': -(0 + 1) / (1 + 0 + 1 + 0), + }, + ]) +def test_score_state(params): + score = pykoop.score_state( + params['X_predicted'], + params['X_expected'], + params['n_steps'], + params['discount_factor'], + params['regression_metric'], + params['min_samples'], + params['episode_feature'], + ) + np.testing.assert_allclose(score, params['score_exp']) + + +@pytest.mark.parametrize('X, min_samples, n_inputs, episode_feature, ic_exp', [ + ( + np.array([ + [1, 2, 3, 4], + [4, 5, 6, 7], + ]).T, + 1, + 0, + False, + np.array([ + [1], + [4], + ]).T, + ), + ( + np.array([ + [1, 2, 3, 4], + [4, 5, 6, 7], + ]).T, + 2, + 0, + False, + np.array([ + [1, 2], + [4, 5], + ]).T, + ), + ( + np.array([ + [1, 2, 3, 4], + [4, 5, 6, 7], + [5, 5, 5, 5], + ]).T, + 1, + 1, + False, + np.array([ + [1], + [4], + ]).T, + ), + ( + np.array([ + [0, 0, 1, 1], + [1, 2, 3, 4], + [4, 5, 6, 7], + ]).T, + 1, + 0, + True, + np.array([ + [0, 1], + [1, 3], + [4, 6], + ]).T, + ), + ( + np.array([ + [0, 0, 1, 1], + [1, 2, 3, 4], + [4, 5, 6, 7], + [9, 9, 9, 9], + ]).T, + 1, + 1, + True, + np.array([ + [0, 1], + [1, 3], + [4, 6], + ]).T, + ), + ( + np.array([ + [0, 0, 0, 1, 1, 1], + [1, 2, 2, 3, 4, 5], + [4, 5, 5, 6, 7, 6], + [9, 9, 9, 9, 9, 6], + ]).T, + 2, + 1, + True, + np.array([ + [0, 0, 1, 1], + [1, 2, 3, 4], + [4, 5, 6, 7], + ]).T, + ), +]) +def test_extract_initial_conditions( + X, + min_samples, + n_inputs, + episode_feature, + ic_exp, +): + ic = pykoop.extract_initial_conditions( + X, + min_samples, + n_inputs, + episode_feature, + ) + np.testing.assert_allclose(ic, ic_exp) + + +@pytest.mark.parametrize('X, n_inputs, episode_feature, u_exp', [ + ( + np.array([ + [1, 2, 3, 4], + [6, 7, 8, 9], + ]).T, + 1, + False, + np.array([ + [6, 7, 8, 9], + ]).T, + ), + ( + np.array([ + [1, 2, 3, 4], + [6, 7, 8, 9], + ]).T, + 0, + False, + np.array([]).reshape((0, 4)).T, + ), + ( + np.array([ + [0, 0, 1, 1], + [1, 2, 3, 4], + [6, 7, 8, 9], + ]).T, + 1, + True, + np.array([ + [0, 0, 1, 1], + [6, 7, 8, 9], + ]).T, + ), +]) +def test_extract_input(X, n_inputs, episode_feature, u_exp): + u = pykoop.extract_input(X, n_inputs, episode_feature) + np.testing.assert_allclose(u, u_exp) + + +@pytest.mark.parametrize( + 'X, n_steps, discount_factor, episode_feature, w_exp', + [ + ( + np.array([ + [1, 2, 3, 4], + [5, 6, 7, 8], + ]).T, + 2, + 1, + False, + np.array([1, 1, 0, 0]), + ), + ( + np.array([ + [1, 2, 3, 4], + [5, 6, 7, 8], + ]).T, + 3, + 0.5, + False, + np.array([1, 0.5, 0.25, 0]), + ), + ( + np.array([ + [0, 0, 0, 1, 1, 1, 1], + [1, 2, 3, 4, 5, 6, 7], + [5, 6, 7, 8, 9, 10, 11], + ]).T, + 2, + 0.5, + True, + np.array([1, 0.5, 0, 1, 0.5, 0, 0]), + ), + ( + np.array([ + [0, 0, 0, 1, 1, 1, 1], + [1, 2, 3, 4, 5, 6, 7], + [5, 6, 7, 8, 9, 10, 11], + ]).T, + 10, + 1, + True, + np.array([1, 1, 1, 1, 1, 1, 1]), + ), + ( + np.array([ + [0, 0, 0, 1, 1, 1, 1], + [1, 2, 3, 4, 5, 6, 7], + [5, 6, 7, 8, 9, 10, 11], + ]).T, + 10, + 0.1, + True, + np.array([1, 0.1, 0.01, 1, 0.1, 0.01, 0.001]), + ), + ], +) +def test_weights_from_data_matrix( + X, + n_steps, + discount_factor, + episode_feature, + w_exp, +): + w = pykoop.koopman_pipeline._weights_from_data_matrix( + X, + n_steps, + discount_factor, + episode_feature, + ) + np.testing.assert_allclose(w, w_exp) + + @pytest.mark.parametrize('kp', [ pykoop.KoopmanPipeline( lifting_functions=[ @@ -613,6 +1019,110 @@ def u(t): ] +@pytest.mark.parametrize('kp', [ + pykoop.KoopmanPipeline( + lifting_functions=[ + ('dl', pykoop.DelayLiftingFn(n_delays_state=1, n_delays_input=1)) + ], + regressor=pykoop.Edmd(), + ), + pykoop.KoopmanPipeline( + lifting_functions=[ + ('dla', pykoop.DelayLiftingFn(n_delays_state=2, n_delays_input=2)), + ('dlb', pykoop.DelayLiftingFn(n_delays_state=2, n_delays_input=2)), + ], + regressor=pykoop.Edmd(), + ), + pykoop.KoopmanPipeline( + lifting_functions=[ + ('dla', pykoop.DelayLiftingFn(n_delays_state=2, n_delays_input=1)), + ('ply', pykoop.PolynomialLiftingFn(order=2)), + ('dlb', pykoop.DelayLiftingFn(n_delays_state=1, n_delays_input=2)), + ], + regressor=pykoop.Edmd(), + ), + pykoop.KoopmanPipeline( + lifting_functions=[ + ( + 'sp', + pykoop.SplitPipeline( + lifting_functions_state=[ + ('pl', pykoop.PolynomialLiftingFn(order=2)), + ], + lifting_functions_input=None, + ), + ), + ('dl', pykoop.DelayLiftingFn(n_delays_state=2, n_delays_input=2)), + ], + regressor=pykoop.Edmd(), + ), + pykoop.KoopmanPipeline( + lifting_functions=[ + ('dla', pykoop.DelayLiftingFn(n_delays_state=1, n_delays_input=1)), + ( + 'sp', + pykoop.SplitPipeline( + lifting_functions_state=[ + ('pla', pykoop.PolynomialLiftingFn(order=2)), + ('plb', pykoop.PolynomialLiftingFn(order=2)), + ], + lifting_functions_input=[ + ('plc', pykoop.PolynomialLiftingFn(order=2)), + ], + ), + ), + ('dlb', pykoop.DelayLiftingFn(n_delays_state=1, n_delays_input=1)), + ], + regressor=pykoop.Edmd(), + ), +]) +def test_predict_state(kp): + # Set up problem + t_range = (0, 5) + t_step = 0.1 + msd = pykoop.dynamic_models.MassSpringDamper(0.5, 0.7, 0.6) + + def u(t): + return 0.1 * np.sin(t) + + # Solve ODE for training data + x0 = np.array([0, 0]) + t, x = msd.simulate(t_range, t_step, x0, u, rtol=1e-8, atol=1e-8) + + # Compute input at every training point + u_sim = np.reshape(u(t), (-1, 1)) + # Format data + X = np.hstack(( + np.zeros((t.shape[0] - 1, 1)), + x[:-1, :], + u_sim[:-1, :], + )) + # Fit estimator + kp.fit(X, n_inputs=1, episode_feature=True) + n_samp = kp.min_samples_ + + X_sim = kp.predict_state( + x[:n_samp, :], + u_sim, + episode_feature=False, + relift_state=True, + ) + + # Predict manually + X_sim_exp = np.empty(x.shape) + X_sim_exp[:, :n_samp] = x[:, :n_samp] + for k in range(n_samp, t.shape[0]): + X = np.hstack(( + np.zeros((n_samp, 1)), + X_sim_exp[(k - n_samp):k, :], + u_sim[(k - n_samp):k, :], + )) + Xp = kp.predict(X) + X_sim_exp[[k], :] = Xp[[-1], 1:] + + np.testing.assert_allclose(X_sim, X_sim_exp) + + @pytest.mark.parametrize('lf, X, Xt_exp, n_inputs, episode_feature, attr_exp', test_split_lf_params) def test_split_lifting_fn_attrs(lf, X, Xt_exp, n_inputs, episode_feature, @@ -661,17 +1171,15 @@ def test_strip_initial_conditons(): np.testing.assert_allclose(X1s, X2s) -@pytest.fixture( - params=[ - pykoop.PolynomialLiftingFn(order=2), - pykoop.KoopmanPipeline( - lifting_functions=[ - ('p', pykoop.PolynomialLiftingFn(order=2)), - ], - regressor=pykoop.Edmd(), - ) - ] -) +@pytest.fixture(params=[ + pykoop.PolynomialLiftingFn(order=2), + pykoop.KoopmanPipeline( + lifting_functions=[ + ('p', pykoop.PolynomialLiftingFn(order=2)), + ], + regressor=pykoop.Edmd(), + ) +]) def lift_retract_fixture(request): lf = request.param X = np.array([ diff --git a/tests/test_regressors.py b/tests/test_regressors.py index c7bd56b..df9b011 100644 --- a/tests/test_regressors.py +++ b/tests/test_regressors.py @@ -10,7 +10,6 @@ # TODO This file is a nightmare - @pytest.fixture( params=[ (pykoop.Edmd(), 'msd-no-input', 1e-5, 1e-5, 'exact'),