diff --git a/README.md b/README.md index 574a947..0c5becd 100644 --- a/README.md +++ b/README.md @@ -130,8 +130,8 @@ permutation_test(OPLS(n_orthogonal=2), X, y) > [!NOTE] > **Pipeline support in plotting:** Diagnostic plotting displays support `OPLS`, -> `OPLSDA`, pipelines ending in one, and fitted search meta-estimators exposing -> `best_estimator_` around either shape. When passing a pipeline, pass raw `X` as +> `OPLSDA`, and pipelines ending in one. For tuned models, pass +> `search.best_estimator_` explicitly. When passing a pipeline, pass raw `X` as > expected by the pipeline. When passing the final OPLS step directly, pass the > already transformed matrix. For pipeline S-plots, points are in the transformed > feature space received by the final OPLS step. diff --git a/docs/quickstart.md b/docs/quickstart.md index b71b5d5..9036082 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -125,9 +125,8 @@ permutation_test(OPLS(n_orthogonal=2), X, y) ``` !!! note "Pipeline support in plotting" -Diagnostic plotting displays support `OPLS`, `OPLSDA`, pipelines ending in one, -and fitted search meta-estimators exposing `best_estimator_` around either shape. -When passing a pipeline, pass raw `X` as expected by the pipeline. When passing the -final OPLS step directly, pass the already transformed matrix. For pipeline -S-plots, points are in the transformed feature space received by the final OPLS -step. +Diagnostic plotting displays support `OPLS`, `OPLSDA`, and pipelines ending in +one. For tuned models, pass `search.best_estimator_` explicitly. When passing a +pipeline, pass raw `X` as expected by the pipeline. When passing the final OPLS +step directly, pass the already transformed matrix. For pipeline S-plots, points +are in the transformed feature space received by the final OPLS step. diff --git a/src/scikit_opls/_inspection.py b/src/scikit_opls/_inspection.py index a711f97..3db1676 100644 --- a/src/scikit_opls/_inspection.py +++ b/src/scikit_opls/_inspection.py @@ -1,19 +1,9 @@ -"""Internal stateless math for OPLS VIP scores and explained-variance metrics. +"""Stateless math helpers for OPLS explained-variance and VIP diagnostics. -Private module — not part of the public API. The VIP scores are exposed as lazy -``vip_`` / ``ortho_vip_`` properties on :class:`~scikit_opls.OPLS` and -:class:`~scikit_opls.OPLSDA`; these functions compute them from fitted weights. - -VIP (Variable Importance in Projection) is defined in the style of Galindo-Prieto -et al. (2014); these are not intended to reproduce ropls VIP values exactly: - -- predictive VIP is the standard PLS VIP of the predictive model fitted on the - orthogonally filtered X, weighting each component by the Y variance it explains; -- orthogonal VIP is an X-variance-weighted score for the removed orthogonal - components, weighting each component by the X variance it explains. - -For non-empty blocks with positive explained variance, VIP is normalized so that -sum(vip**2) == n_features. Empty or degenerate blocks return zeros. +Private module — not part of the public API. Used by the fitted attributes of +:class:`~scikit_opls.OPLS` and :class:`~scikit_opls.OPLSDA`. VIP scores are +normalized so ``sum(vip**2) == n_features`` when component importance is +positive; degenerate inputs return zeros. """ from __future__ import annotations @@ -27,53 +17,47 @@ def _safe_total_ss(X: NDArray[np.float64]) -> float: """Total sum of squares with a nonzero guard.""" total = float(np.sum(np.asarray(X, dtype=np.float64) ** 2)) - return max(total, np.finfo(np.float64).eps) + return max(total, _EPS) -def component_explained_x_variance( +def _validate_x_scores_loadings( X: NDArray[np.float64], scores: NDArray[np.float64], loadings: NDArray[np.float64], -) -> NDArray[np.float64]: - """Per-component fraction of X sum-of-squares explained. - - X must already be in the model space used to fit the scores/loadings. - """ - X_arr = np.asarray(X, dtype=np.float64) - T = np.asarray(scores, dtype=np.float64) - P = np.asarray(loadings, dtype=np.float64) - if T.ndim != 2 or P.ndim != 2: - raise ValueError("scores and loadings must be 2D arrays.") - if T.shape[1] != P.shape[1]: +) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: + X = np.asarray(X, dtype=np.float64) + scores = np.asarray(scores, dtype=np.float64) + loadings = np.asarray(loadings, dtype=np.float64) + if X.ndim != 2 or scores.ndim != 2 or loadings.ndim != 2: + raise ValueError("X, scores and loadings must all be 2D arrays.") + if scores.shape[0] != X.shape[0]: + raise ValueError("scores must have one row per sample of X.") + if loadings.shape[0] != X.shape[1]: + raise ValueError("loadings must have one row per feature of X.") + if scores.shape[1] != loadings.shape[1]: raise ValueError("scores and loadings must have the same number of components.") - total = _safe_total_ss(X_arr) - out = np.empty(T.shape[1], dtype=np.float64) - for i in range(T.shape[1]): - Xi = T[:, [i]] @ P[:, [i]].T - out[i] = np.sum(Xi**2) / total - return out + if not np.all(np.isfinite(X)): + raise ValueError("X must contain only finite values.") + if not np.all(np.isfinite(scores)): + raise ValueError("scores must contain only finite values.") + if not np.all(np.isfinite(loadings)): + raise ValueError("loadings must contain only finite values.") + return X, scores, loadings -def cumulative_r2_from_residuals( - original: NDArray[np.float64], - residuals_by_component: list[NDArray[np.float64]], -) -> NDArray[np.float64]: - """Cumulative R² from a sequence of residual matrices.""" - total = _safe_total_ss(original) - return np.asarray( - [1.0 - float(np.sum(resid**2)) / total for resid in residuals_by_component], - dtype=np.float64, - ) - - -def component_r2_from_cumulative( - cumulative: NDArray[np.float64], +def component_explained_x_variance( + X: NDArray[np.float64], + scores: NDArray[np.float64], + loadings: NDArray[np.float64], ) -> NDArray[np.float64]: - """Convert cumulative R² to per-component increments.""" - cumulative = np.asarray(cumulative, dtype=np.float64) - if cumulative.size == 0: - return cumulative - return np.diff(np.r_[0.0, cumulative]) + """Per-component ``SS(t_i @ p_i.T) / SS(X)`` for fitted arrays.""" + X, scores, loadings = _validate_x_scores_loadings(X, scores, loadings) + total = _safe_total_ss(X) + out = np.empty(scores.shape[1], dtype=np.float64) + for i in range(scores.shape[1]): + Xi = scores[:, [i]] @ loadings[:, [i]].T + out[i] = np.sum(Xi**2) / total + return out def component_r2y_from_scores( @@ -89,10 +73,28 @@ def component_r2y_from_scores( y_arr = np.asarray(y, dtype=np.float64) if y_arr.ndim == 1: y_arr = y_arr.reshape(-1, 1) + if y_arr.ndim != 2: + raise ValueError(f"y must be 1D or 2D, got shape {y_arr.shape}.") T = np.asarray(scores, dtype=np.float64) + if T.ndim != 2: + raise ValueError(f"scores must be 2D, got shape {T.shape}.") Q = np.asarray(y_loadings, dtype=np.float64) if Q.ndim == 1: - Q = Q.reshape(-1, 1) + # A 1D y_loadings is one value per component (single target), matching the + # (n_targets, n_components) convention used elsewhere (predictive_vip). + Q = Q.reshape(1, -1) + elif Q.ndim != 2: + raise ValueError(f"y_loadings must be 1D or 2D, got shape {Q.shape}.") + if T.shape[0] != y_arr.shape[0]: + raise ValueError("scores must have one row per sample of y.") + if Q.shape[1] != T.shape[1]: + raise ValueError("y_loadings must have one column per component.") + if not np.all(np.isfinite(y_arr)): + raise ValueError("y must contain only finite values.") + if not np.all(np.isfinite(T)): + raise ValueError("scores must contain only finite values.") + if not np.all(np.isfinite(Q)): + raise ValueError("y_loadings must contain only finite values.") total = _safe_total_ss(y_arr - y_arr.mean(axis=0, keepdims=True)) out = np.empty(T.shape[1], dtype=np.float64) for i in range(T.shape[1]): @@ -106,61 +108,22 @@ def explained_x_variance( scores: NDArray[np.float64], loadings: NDArray[np.float64], ) -> float: - """Fraction of the (preprocessed) ``X`` sum-of-squares captured by a block. - - Used for both the predictive block (``R2X``) and the orthogonal block - (``R2X_ortho``): ``SS(scores @ loadingsᵀ) / SS(X)``. - - Parameters - ---------- - X : ndarray of shape (n_samples, n_features) - Preprocessed predictor matrix. - scores : ndarray of shape (n_samples, n_components) - Block scores. - loadings : ndarray of shape (n_features, n_components) - Block loadings. - - Returns - ------- - fraction : float - Captured nominal fraction; ``0.0`` if the block is empty or ``X`` - has zero sum-of-squares. Can exceed 1.0 if the supplied scores/loadings do - not form an orthogonal sum-of-squares decomposition; callers should treat - it as a diagnostic summary. - """ - if X.ndim != 2 or scores.ndim != 2 or loadings.ndim != 2: - raise ValueError("X, scores and loadings must all be 2D arrays.") - if scores.shape[0] != X.shape[0]: - raise ValueError("scores must have one row per sample of X.") - if loadings.shape[0] != X.shape[1]: - raise ValueError("loadings must have one row per feature of X.") - if scores.shape[1] != loadings.shape[1]: - raise ValueError("scores and loadings must have the same number of components.") + """Nominal ``SS(T @ P.T) / SS(X)``; not clipped to ``[0, 1]``.""" + X, scores, loadings = _validate_x_scores_loadings(X, scores, loadings) + if scores.shape[1] == 0: + return 0.0 total = float(np.sum(X**2)) - if total <= 0.0 or scores.shape[1] == 0: + if total <= 0.0: return 0.0 - # Rebuild the part of X represented by this score/loading block and compare - # its sum of squares with the full preprocessed X block. - approx = scores @ loadings.T - return float(np.sum(approx**2) / total) + return float(np.sum((scores @ loadings.T) ** 2) / total) def _weighted_vip( weights: NDArray[np.float64], ss_per_component: NDArray[np.float64] ) -> NDArray[np.float64]: - """VIP from per-component weight vectors and their importance weights. - - Parameters - ---------- - weights : ndarray of shape (n_features, n_components) - Per-component weight vectors. - ss_per_component : ndarray of shape (n_components,) - Non-negative variance explained by each component. + """Return VIP scores from component weights and importance values. - Returns - ------- - vip : ndarray of shape (n_features,) - VIP scores; all-zero when there are no components or zero total variance. + Zeros for empty components or zero total importance. """ if weights.ndim != 2: raise ValueError(f"weights must be 2D, got shape {weights.shape}.") @@ -195,26 +158,15 @@ def predictive_vip( x_scores: NDArray[np.float64], y_loadings: NDArray[np.float64], ) -> NDArray[np.float64]: - """Predictive VIP from the engine's weights/scores/Y-loadings. - - Parameters - ---------- - x_weights : ndarray of shape (n_features, n_components) - Predictive weight vectors. - x_scores : ndarray of shape (n_samples, n_components) - Predictive scores. - y_loadings : ndarray of shape (n_components,) or (1, n_components) - Y-loadings of the predictive components. - - Returns - ------- - vip : ndarray of shape (n_features,) - Predictive VIP scores. - """ + """Return predictive PLS VIP from weights, scores and Y-loadings.""" if x_weights.ndim != 2: raise ValueError(f"x_weights must be 2D, got shape {x_weights.shape}.") if x_scores.ndim != 2: raise ValueError(f"x_scores must be 2D, got shape {x_scores.shape}.") + if not np.all(np.isfinite(x_weights)): + raise ValueError("x_weights must contain only finite values.") + if not np.all(np.isfinite(x_scores)): + raise ValueError("x_scores must contain only finite values.") _, n_components = x_weights.shape if x_scores.shape[1] != n_components: @@ -237,6 +189,8 @@ def predictive_vip( y_loadings_2d = y_loadings else: raise ValueError(f"y_loadings must be 1D or 2D, got shape {y_loadings.shape}.") + if not np.all(np.isfinite(y_loadings_2d)): + raise ValueError("y_loadings must contain only finite values.") # Standard PLS VIP weights each component by the Y sum of squares explained by # that component: loading strength times score energy. @@ -249,22 +203,7 @@ def orthogonal_vip( x_ortho_scores: NDArray[np.float64], x_ortho_loadings: NDArray[np.float64], ) -> NDArray[np.float64]: - """Orthogonal VIP, each component weighted by the X variance it captures. - - Parameters - ---------- - x_ortho_weights : ndarray of shape (n_features, n_orthogonal) - Orthogonal weight vectors. - x_ortho_scores : ndarray of shape (n_samples, n_orthogonal) - Orthogonal scores. - x_ortho_loadings : ndarray of shape (n_features, n_orthogonal) - Orthogonal loadings. - - Returns - ------- - vip : ndarray of shape (n_features,) - Orthogonal VIP scores. - """ + """Return orthogonal VIP weighted by removed X variance.""" if x_ortho_weights.ndim != 2: raise ValueError( f"x_ortho_weights must be 2D, got shape {x_ortho_weights.shape}." @@ -277,6 +216,12 @@ def orthogonal_vip( raise ValueError( f"x_ortho_loadings must be 2D, got shape {x_ortho_loadings.shape}." ) + if not np.all(np.isfinite(x_ortho_weights)): + raise ValueError("x_ortho_weights must contain only finite values.") + if not np.all(np.isfinite(x_ortho_scores)): + raise ValueError("x_ortho_scores must contain only finite values.") + if not np.all(np.isfinite(x_ortho_loadings)): + raise ValueError("x_ortho_loadings must contain only finite values.") n_features, n_components = x_ortho_weights.shape if x_ortho_scores.shape[1] != n_components: diff --git a/src/scikit_opls/_o2pls.py b/src/scikit_opls/_o2pls.py index 849ccc2..e151638 100644 --- a/src/scikit_opls/_o2pls.py +++ b/src/scikit_opls/_o2pls.py @@ -29,7 +29,7 @@ o2pls_fit, ) from scikit_opls._preprocessing import VALID_SCALING, apply_scaling, compute_scaling -from scikit_opls._utils import _has_nonzero_variation +from scikit_opls._utils import _has_nonzero_variation, _reject_bool_param class O2PLS(RegressorMixin, TransformerMixin, BaseEstimator): @@ -166,12 +166,9 @@ def fit(self, X: ArrayLike, Y: ArrayLike) -> O2PLS: self : object Fitted estimator. """ - if isinstance(self.n_components, bool): - raise ValueError("n_components must be an integer, not bool.") - if isinstance(self.n_x_orthogonal, bool): - raise ValueError("n_x_orthogonal must be an integer, not bool.") - if isinstance(self.n_y_orthogonal, bool): - raise ValueError("n_y_orthogonal must be an integer, not bool.") + _reject_bool_param("n_components", self.n_components) + _reject_bool_param("n_x_orthogonal", self.n_x_orthogonal) + _reject_bool_param("n_y_orthogonal", self.n_y_orthogonal) self._validate_params() # Preserve whether the user supplied 1D or 2D Y so predict() can mirror the diff --git a/src/scikit_opls/_o2pls_core.py b/src/scikit_opls/_o2pls_core.py index 3e5f5da..5fdab73 100644 --- a/src/scikit_opls/_o2pls_core.py +++ b/src/scikit_opls/_o2pls_core.py @@ -4,14 +4,13 @@ import warnings from dataclasses import dataclass -from numbers import Integral import numpy as np from numpy.typing import NDArray from sklearn.exceptions import ConvergenceWarning from sklearn.utils.extmath import svd_flip -from scikit_opls._utils import _has_nonzero_variation +from scikit_opls._utils import _has_nonzero_variation, _validate_int _TOL = 1e-12 @@ -73,22 +72,11 @@ def _effective_rank(s: NDArray[np.float64], tol: float) -> int: def _validate_positive_int(name: str, value: int) -> int: - # ``bool`` is a subclass of ``int``; reject it explicitly so True/False are not - # silently accepted as component counts. - if isinstance(value, bool) or not isinstance(value, Integral): - raise TypeError(f"{name} must be an integer, got {type(value).__name__}.") - if value < 1: - raise ValueError(f"{name} must be >= 1, got {value}.") - return int(value) + return _validate_int(name, value, minimum=1) def _validate_nonnegative_int(name: str, value: int) -> int: - # Keep the same bool handling as the positive-integer validator above. - if isinstance(value, bool) or not isinstance(value, Integral): - raise TypeError(f"{name} must be an integer, got {type(value).__name__}.") - if value < 0: - raise ValueError(f"{name} must be >= 0, got {value}.") - return int(value) + return _validate_int(name, value, minimum=0) def _validate_tol(tol: float) -> float: @@ -106,25 +94,10 @@ def _validate_tol(tol: float) -> float: def _cross_cov_svd_x_to_y( Xs: NDArray[np.float64], Ys: NDArray[np.float64], k: int ) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: - """SVD of ``Xs.T @ Ys`` with paired deterministic signs. + """Deterministic-sign SVD of ``Xs.T @ Ys``. - Parameters - ---------- - Xs : ndarray of shape (n_samples, n_x_features) - Scaled X block. - Ys : ndarray of shape (n_samples, n_y_features) - Scaled Y block. - k : int - Number of components to compute. - - Returns - ------- - W : ndarray of shape (n_x_features, k) - X-side weights. - C : ndarray of shape (n_y_features, k) - Y-side weights. - s : ndarray - Singular values from the thin cross-covariance SVD. + Returns X-side weights ``W`` of shape (n_x_features, k), Y-side weights + ``C`` of shape (n_y_features, k), and all singular values from the thin SVD. """ k = _validate_nonnegative_int("k", k) X = np.asarray(Xs, dtype=np.float64) @@ -159,22 +132,10 @@ def _cross_cov_svd_x_to_y( def _lstsq_map( scores: NDArray[np.float64], block: NDArray[np.float64] ) -> NDArray[np.float64]: - """Return the least-squares map B solving scores @ B ≈ block. - - The returned array has shape (n_score_columns, n_block_columns). - This is not transposed into a feature-by-component loading matrix. + """Least-squares map ``B`` solving ``scores @ B ≈ block``. - Parameters - ---------- - scores : ndarray of shape (n_samples, n_score_columns) - Predictor scores. - block : ndarray of shape (n_samples, n_block_columns) - Target block. - - Returns - ------- - coef : ndarray of shape (n_score_columns, n_block_columns) - Least-squares coefficient matrix. + Shape is (n_score_columns, n_block_columns); not transposed into a + feature-by-component loading matrix. """ scores = np.asarray(scores, dtype=np.float64) block = np.asarray(block, dtype=np.float64) @@ -199,7 +160,7 @@ def _extract_one_orthogonal_component( *, tol: float = _TOL, ) -> OrthogonalBlockComponent | None: - """Extract one sequential O2PLS orthogonal component from ``block``. + """Extract one replayable sequential orthogonal component from ``block``. The residual first removes the enlarged preliminary joint subspace from the current block. The leading left singular vector of @@ -207,22 +168,7 @@ def _extract_one_orthogonal_component( block-specific variation most associated with the preliminary joint score space. The resulting score/loading pair is deflated from the current block and stored so the same sequential filter can be replayed on new data. - - Parameters - ---------- - block : ndarray of shape (n_samples, n_features) - Current data block to deflate. - joint_scores : ndarray of shape (n_samples, n_components) - Preliminary joint scores. - joint_weights : ndarray of shape (n_features, n_components) - Preliminary joint weights. - tol : float, default=_TOL - Tolerance for numerical rank and variance deflation. - - Returns - ------- - component : OrthogonalBlockComponent or None - The extracted component, or None if no resolvable variation remains. + Returns ``None`` if no resolvable variation remains. """ tol = _validate_tol(tol) X = np.asarray(block, dtype=np.float64) @@ -267,7 +213,8 @@ def _extract_one_orthogonal_component( # sequential filter that will later be replayed on new samples. score = X @ weight score_ssq = float(score @ score) - block_ssq = max(_ssq(X), 1.0) + x_ssq = _ssq(X) + block_ssq = max(x_ssq, 1.0) if score_ssq <= tol * block_ssq: return None @@ -279,7 +226,7 @@ def _extract_one_orthogonal_component( if not np.all(np.isfinite(filtered)): return None # Refuse components that do not measurably reduce the block sum of squares. - if _ssq(filtered) >= _ssq(X) - tol * max(_ssq(X), 1.0): + if _ssq(filtered) >= x_ssq - tol * block_ssq: return None return OrthogonalBlockComponent( @@ -348,12 +295,12 @@ def _replay_orthogonal_filter( def _stack_columns( - values: list[NDArray[np.float64]], n_rows: int, n_columns: int + values: list[NDArray[np.float64]], n_rows: int ) -> NDArray[np.float64]: """Column-stack values or return a correctly shaped empty matrix.""" if values: return np.column_stack(values) - return np.zeros((n_rows, n_columns), dtype=np.float64) + return np.zeros((n_rows, 0), dtype=np.float64) def _reconstruction_r2( @@ -375,27 +322,10 @@ def o2pls_fit( *, tol: float = _TOL, ) -> O2PLSComponents: - """Fit dense O2PLS components on already preprocessed blocks. + """Fit dense O2PLS components on already preprocessed X/Y blocks. - Parameters - ---------- - Xs : ndarray of shape (n_samples, n_x_features) - Preprocessed X block. - Ys : ndarray of shape (n_samples, n_y_features) - Preprocessed Y block. - n_components : int - Number of joint components. - n_x_orthogonal : int - Number of X-specific orthogonal components. - n_y_orthogonal : int - Number of Y-specific orthogonal components. - tol : float, default=_TOL - Numerical tolerance. - - Returns - ------- - fit : O2PLSComponents - Dataclass containing all fitted matrices and diagnostics. + Uses an enlarged preliminary joint subspace, sequential X/Y orthogonal + filtering, and final joint SVD re-estimation on the filtered blocks. """ X0 = np.asarray(Xs, dtype=np.float64) Y0 = np.asarray(Ys, dtype=np.float64) @@ -505,12 +435,12 @@ def o2pls_fit( B_T = _lstsq_map(T, U) B_U = _lstsq_map(U, T) - X_orth_weights = _stack_columns(x_weights, n_x_features, 0) - X_orth_scores = _stack_columns(x_scores, n_samples, 0) - X_orth_loadings = _stack_columns(x_loadings, n_x_features, 0) - Y_orth_weights = _stack_columns(y_weights, n_y_features, 0) - Y_orth_scores = _stack_columns(y_scores, n_samples, 0) - Y_orth_loadings = _stack_columns(y_loadings, n_y_features, 0) + X_orth_weights = _stack_columns(x_weights, n_x_features) + X_orth_scores = _stack_columns(x_scores, n_samples) + X_orth_loadings = _stack_columns(x_loadings, n_x_features) + Y_orth_weights = _stack_columns(y_weights, n_y_features) + Y_orth_scores = _stack_columns(y_scores, n_samples) + Y_orth_loadings = _stack_columns(y_loadings, n_y_features) # Reconstruct nominal joint, orthogonal and residual parts in the original # preprocessed coordinates for diagnostics and estimator attributes. diff --git a/src/scikit_opls/_opls.py b/src/scikit_opls/_opls.py index 517eceb..f865c79 100644 --- a/src/scikit_opls/_opls.py +++ b/src/scikit_opls/_opls.py @@ -1,12 +1,10 @@ """Orthogonal PLS (OPLS) regressor with a scikit-learn interface. -OPLS (Trygg & Wold, 2002) splits the variation in ``X`` into a *predictive* part -correlated with ``y`` and *orthogonal* parts uncorrelated with ``y``. This -estimator removes the orthogonal variation with an OSC-style orthogonal filter -(:mod:`scikit_opls._orthogonal`) and then fits -:class:`sklearn.cross_decomposition.PLSRegression` on the cleaned ``X`` as the -predictive engine. With ``n_orthogonal=0``, this becomes ordinary -``PLSRegression`` after the package's selected X preprocessing. +OPLS removes X-orthogonal variation with an OSC-style filter +(:mod:`scikit_opls._orthogonal`), then fits +:class:`sklearn.cross_decomposition.PLSRegression` on the filtered X block. +With ``n_orthogonal=0``, this reduces to ordinary PLS after this package's +selected X preprocessing. """ # scikit-learn's validate_data uses sentinel-string parameter defaults that lead @@ -23,6 +21,7 @@ from __future__ import annotations +from dataclasses import dataclass from numbers import Integral import numpy as np @@ -46,19 +45,24 @@ ) from scikit_opls._orthogonal import apply_orthogonal_filter, opls_filter from scikit_opls._preprocessing import VALID_SCALING, apply_scaling, compute_scaling -from scikit_opls._utils import _has_nonzero_variation +from scikit_opls._utils import _has_nonzero_variation, _reject_bool_param + + +@dataclass(frozen=True) +class _OPLSProjection: + """All fitted OPLS model-space arrays for one validated raw X block.""" + + Xs: NDArray[np.float64] + X_filtered: NDArray[np.float64] + t_pred: NDArray[np.float64] + t_ortho: NDArray[np.float64] def _orthogonal_filter_matrix( x_ortho_weights: NDArray[np.float64], x_ortho_loadings: NDArray[np.float64], ) -> NDArray[np.float64]: - """Right-side linear operator ``F`` such that ``X_filtered == X_scaled @ F``. - - The replayed orthogonal filter applies ``X <- X - (X w_i) p_iᵀ`` for each - component, i.e. right multiplication by ``I - outer(w_i, p_i)``. Composing them - in order yields the single matrix equivalent of :func:`apply_orthogonal_filter`. - """ + """Right-side matrix ``F`` such that ``X_filtered == X_scaled @ F``.""" W = np.asarray(x_ortho_weights, dtype=np.float64) P = np.asarray(x_ortho_loadings, dtype=np.float64) n_features = W.shape[0] @@ -80,11 +84,8 @@ def _compose_raw_coefficients( ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: """Collapse filtered/scaled-space PLS coefficients into raw-X coefficients. - The fitted prediction is linear: ``X -> (X - mean) / std -> @ F -> @ Bᶠ + b``, - where ``b`` is the predictive engine's prediction offset (``pls_.predict(0)``, - not ``pls_.intercept_``). With ``B_scaled = F @ Bᶠ`` and ``inv_scale = 1 / std`` - this reduces to ``y = X @ B_raw + (b - (mean * inv_scale) @ B_scaled)`` where - ``B_raw = inv_scale[:, None] * B_scaled``. + ``intercept_filtered`` must be the predictive engine's prediction at zero + filtered input (``pls_.predict(0)``), not necessarily its public ``intercept_``. """ coef_arr = np.asarray(coef_filtered, dtype=np.float64) if coef_arr.ndim == 1: @@ -253,21 +254,33 @@ def fit(self, X: ArrayLike, y: ArrayLike) -> OPLS: self : OPLS The fitted estimator. """ - # Refits must invalidate lazily cached VIP arrays because fitted weights may - # change even when the same Python estimator instance is reused. - for _attr in ("_vip_", "_ortho_vip_"): - self.__dict__.pop(_attr, None) - if isinstance(self.n_components, bool): - raise ValueError("n_components must be an integer, not bool.") - if isinstance(self.n_orthogonal, bool): - raise ValueError("n_orthogonal must be an integer, not bool.") + self._clear_fit_caches() + X, y = self._validate_fit_data(X, y) + Xs, X_filtered = self._fit_orthogonal_filter(X, y) + self._fit_predictive_engine(X_filtered, y) + self._set_raw_coefficients(X_filtered) + self._set_fit_diagnostics(X, Xs, X_filtered, y) + return self + + def _clear_fit_caches(self) -> None: + for attr in ("_vip_", "_ortho_vip_"): + self.__dict__.pop(attr, None) + + def _validate_fit_data( + self, + X: ArrayLike, + y: ArrayLike, + ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: + _reject_bool_param("n_components", self.n_components) + _reject_bool_param("n_orthogonal", self.n_orthogonal) self._validate_params() - # multi_output=False ravels a column-vector y (with a DataConversionWarning) - # and rejects multi-column y: OPLS is univariate. - # Note: The internal OSC primitive supports multivariate blocks, but the public - # OPLS estimator currently exposes only univariate OPLS regression. X, y = validate_data( - self, X, y, dtype=np.float64, ensure_min_samples=2, copy=self.copy + self, + X, + y, + dtype=np.float64, + ensure_min_samples=2, + copy=self.copy, ) if not _has_nonzero_variation(y): raise ValueError("OPLS requires a non-constant target y.") @@ -277,26 +290,34 @@ def fit(self, X: ArrayLike, y: ArrayLike) -> OPLS: f"n_components={self.n_components} exceeds the maximum of " f"min(n_samples, n_features)={min(X.shape)}." ) + return X, y + def _fit_orthogonal_filter( + self, + X: NDArray[np.float64], + y: NDArray[np.float64], + ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: self.x_mean_, self.x_std_ = compute_scaling(X, self.scale) Xs = apply_scaling(X, self.x_mean_, self.x_std_) - # The predictive direction and downstream PLS model both require resolvable - # variation after preprocessing. if not _has_nonzero_variation(Xs, axis=0): raise ValueError("X has no non-zero variation after preprocessing.") - - # opls_filter handles n_orthogonal=0 (pass-through) and truncation itself. - # y is centered only for direction extraction; the predictive PLS engine - # below still fits against the original y scale. ofit = opls_filter(Xs, y - y.mean(), self.n_orthogonal) - X_filtered = ofit.x_filtered + self.x_ortho_weights_ = ofit.x_ortho_weights + self.x_ortho_loadings_ = ofit.x_ortho_loadings + self.x_ortho_scores_ = ofit.x_ortho_scores + self.n_orthogonal_ = ofit.n_components + return Xs, ofit.x_filtered + + def _fit_predictive_engine( + self, + X_filtered: NDArray[np.float64], + y: NDArray[np.float64], + ) -> None: if not _has_nonzero_variation(X_filtered, axis=0): raise ValueError( "X has no remaining variation after orthogonal filtering; " "reduce n_orthogonal." ) - - # Validate the numerical rank of the actual matrix passed to PLSRegression. rank_filtered = np.linalg.matrix_rank(X_filtered) if self.n_components > rank_filtered: raise ValueError( @@ -304,29 +325,30 @@ def fit(self, X: ArrayLike, y: ArrayLike) -> OPLS: f"X after orthogonal filtering ({rank_filtered}). " "Reduce n_components or n_orthogonal." ) - - self.x_ortho_weights_ = ofit.x_ortho_weights - self.x_ortho_loadings_ = ofit.x_ortho_loadings - self.x_ortho_scores_ = ofit.x_ortho_scores - self.n_orthogonal_ = ofit.n_components - - # The sklearn PLS engine sees only the preprocessed, orthogonally filtered X. self.pls_ = PLSRegression(n_components=self.n_components, scale=False) self.pls_.fit(X_filtered, y) - # transform() returns the predictive scores: one output column per component. self._n_features_out = self.n_components - - # Surface the predictive model parameters from the engine. self.x_weights_ = self.pls_.x_weights_ self.x_loadings_ = self.pls_.x_loadings_ self.x_scores_ = self.pls_.x_scores_ self.y_loadings_ = self.pls_.y_loadings_ self.coef_filtered_ = self.pls_.coef_ self.intercept_ = self.pls_.intercept_ - # The engine's prediction offset is predict(0), not intercept_: sklearn's - # PLSRegression centers the filtered X internally, so predict(Z) == - # Z @ coef_.T + predict(0) but intercept_ omits that centering term (it only - # coincides with predict(0) when the filtered X is already centered). + # PLSRegression centers the filtered X block internally; reconstructing + # scaled X from predictive scores (``_predictive_x_hat``) needs that fitted + # mean back. sklearn does not expose it publicly, so we intentionally read + # the private ``_x_mean`` attribute and pin the behavior with regression + # tests (test_opls_coefficients.py) rather than relying on it silently. + try: + self._pls_x_mean = np.asarray(self.pls_._x_mean, dtype=np.float64) + except AttributeError as exc: + raise AttributeError( + "Could not access PLSRegression._x_mean. OPLS's reconstruction of " + "scaled X from predictive scores depends on sklearn's fitted PLS " + "centering state." + ) from exc + + def _set_raw_coefficients(self, X_filtered: NDArray[np.float64]) -> None: engine_offset = self.pls_.predict( np.zeros((1, X_filtered.shape[1]), dtype=np.float64) ).ravel() @@ -339,8 +361,13 @@ def fit(self, X: ArrayLike, y: ArrayLike) -> OPLS: self.x_ortho_loadings_, ) - # Diagnostics are computed on the training data in the same coordinate - # systems exposed by the fitted attributes. + def _set_fit_diagnostics( + self, + X: NDArray[np.float64], + Xs: NDArray[np.float64], + X_filtered: NDArray[np.float64], + y: NDArray[np.float64], + ) -> None: y_fit = self.pls_.predict(X_filtered) self.r2x_ = explained_x_variance(Xs, self.x_scores_, self.x_loadings_) self.r2x_ortho_ = explained_x_variance( @@ -365,13 +392,12 @@ def fit(self, X: ArrayLike, y: ArrayLike) -> OPLS: self.y_loadings_, ) - X_model_full, X_hat_full = self._full_reconstruction_validated(X) - self.q_residuals_train_ = np.sum((X_model_full - X_hat_full) ** 2, axis=1) + proj = self._project_validated(X) + self.q_residuals_train_ = self._q_residuals_from_projection(proj, space="full") self.x_residual_ss_ = float(np.sum(self.q_residuals_train_)) - X_model_pred, X_hat_pred = self._predictive_reconstruction_validated(X) - self.q_residuals_predictive_train_ = np.sum( - (X_model_pred - X_hat_pred) ** 2, axis=1 + self.q_residuals_predictive_train_ = self._q_residuals_from_projection( + proj, space="predictive" ) y_fit_arr = np.asarray(y_fit, dtype=np.float64) @@ -380,7 +406,25 @@ def fit(self, X: ArrayLike, y: ArrayLike) -> OPLS: y_arr = y_arr.reshape(-1, 1) self.y_residual_ss_ = float(np.sum((y_arr - y_fit_arr) ** 2)) - return self + def _project_validated(self, X_valid: NDArray[np.float64]) -> _OPLSProjection: + """Project already validated raw X into fitted OPLS model spaces.""" + Xs = apply_scaling(X_valid, self.x_mean_, self.x_std_) + X_filtered, t_ortho = apply_orthogonal_filter( + Xs, + self.x_ortho_weights_, + self.x_ortho_loadings_, + ) + t_pred = self.pls_.transform(X_filtered) + return _OPLSProjection( + Xs=Xs, + X_filtered=X_filtered, + t_pred=t_pred, + t_ortho=t_ortho, + ) + + def _predictive_x_hat(self, t_pred: NDArray[np.float64]) -> NDArray[np.float64]: + """Reconstruct scaled X from predictive scores.""" + return self._pls_x_mean + t_pred @ self.x_loadings_.T def predict(self, X: ArrayLike) -> NDArray[np.float64]: """Predict ``y`` for new samples. @@ -396,8 +440,8 @@ def predict(self, X: ArrayLike) -> NDArray[np.float64]: Predicted target values. """ X_valid = self._validate_X_predict(X) - X_filtered, _ = self._filter_validated(X_valid) - return self.pls_.predict(X_filtered).ravel() + proj = self._project_validated(X_valid) + return self.pls_.predict(proj.X_filtered).ravel() def transform(self, X: ArrayLike) -> NDArray[np.float64]: """Project samples onto the predictive components. @@ -413,8 +457,7 @@ def transform(self, X: ArrayLike) -> NDArray[np.float64]: Predictive scores. """ X_valid = self._validate_X_predict(X) - X_filtered, _ = self._filter_validated(X_valid) - return self.pls_.transform(X_filtered) + return self._project_validated(X_valid).t_pred def transform_orthogonal(self, X: ArrayLike) -> NDArray[np.float64]: """Project samples onto the orthogonal components. @@ -432,7 +475,7 @@ def transform_orthogonal(self, X: ArrayLike) -> NDArray[np.float64]: Orthogonal scores. """ X_valid = self._validate_X_predict(X) - return self._filter_validated(X_valid)[1] + return self._project_validated(X_valid).t_ortho def filter_transform(self, X: ArrayLike) -> NDArray[np.float64]: """Return ``X`` after preprocessing and orthogonal filtering. @@ -454,7 +497,7 @@ def filter_transform(self, X: ArrayLike) -> NDArray[np.float64]: Preprocessed ``X`` with the fitted orthogonal variation removed. """ X_valid = self._validate_X_predict(X) - return self._filter_validated(X_valid)[0] + return self._project_validated(X_valid).X_filtered @property def vip_(self) -> NDArray[np.float64]: @@ -514,29 +557,6 @@ def get_feature_names_out(self, input_features=None) -> NDArray[np.object_]: [f"opls_pred{i}" for i in range(self._n_features_out)], dtype=object ) - def score(self, X: ArrayLike, y: ArrayLike, sample_weight=None) -> float: - """Coefficient of determination R² of the prediction. - - Inherited from :class:`~sklearn.base.RegressorMixin` (``OPLS`` is also a - :class:`~sklearn.base.TransformerMixin`; the regression ``score`` applies, - not a transformer score). - - Parameters - ---------- - X : array-like of shape (n_samples, n_features) - Test samples. - y : array-like of shape (n_samples,) - True target values for ``X``. - sample_weight : array-like of shape (n_samples,), default=None - Sample weights. - - Returns - ------- - score : float - R² of ``self.predict(X)`` against ``y``. - """ - return super().score(X, y, sample_weight) - def _validate_X_predict(self, X: ArrayLike) -> NDArray[np.float64]: # noqa: N802 """Validate prediction/projection input against fitted OPLS metadata.""" check_is_fitted(self) @@ -548,17 +568,6 @@ def _validate_X_predict(self, X: ArrayLike) -> NDArray[np.float64]: # noqa: N80 reset=False, ) - def _filter_validated( - self, X: NDArray[np.float64] - ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: - """Filter an already validated dense array without checking names again.""" - Xs = apply_scaling(X, self.x_mean_, self.x_std_) - # apply_orthogonal_filter returns both the filtered matrix for prediction - # and the replayed orthogonal scores for transform_orthogonal(). - return apply_orthogonal_filter( - Xs, self.x_ortho_weights_, self.x_ortho_loadings_ - ) - def _score_distance_from_scores( self, scores: NDArray[np.float64], @@ -613,68 +622,26 @@ def _score_distance_validated( kind: str = "predictive", ) -> NDArray[np.float64]: """Compute score distance internally, assuming X is already validated.""" + proj = self._project_validated(X_valid) if kind == "predictive": - X_filtered, _ = self._filter_validated(X_valid) - scores = self.pls_.transform(X_filtered) + scores = proj.t_pred reference = self.x_scores_ elif kind == "orthogonal": if self.n_orthogonal_ == 0: return np.zeros(X_valid.shape[0], dtype=np.float64) - scores = self._filter_validated(X_valid)[1] + scores = proj.t_ortho reference = self.x_ortho_scores_ elif kind == "all": - X_filtered, ortho = self._filter_validated(X_valid) - pred = self.pls_.transform(X_filtered) if self.n_orthogonal_ == 0: - scores = pred + scores = proj.t_pred reference = self.x_scores_ else: - scores = np.hstack([pred, ortho]) + scores = np.hstack([proj.t_pred, proj.t_ortho]) reference = np.hstack([self.x_scores_, self.x_ortho_scores_]) else: raise ValueError("kind must be one of {'predictive', 'orthogonal', 'all'}.") return self._score_distance_from_scores(scores, reference) - def _pls_x_mean(self) -> NDArray[np.float64]: - """Return the fitted PLS engine's internal X mean.""" - x_mean = getattr(self.pls_, "_x_mean", None) - if x_mean is None: - raise AttributeError( - "The fitted PLS engine does not expose _x_mean; " - "cannot reconstruct predictive X-space residuals." - ) - return np.asarray(x_mean, dtype=np.float64) - - def _scaled_x_validated(self, X_valid: NDArray[np.float64]) -> NDArray[np.float64]: - """Return X in the fitted scaled/centered model space.""" - return apply_scaling(X_valid, self.x_mean_, self.x_std_) - - def _predictive_reconstruction_validated( - self, - X_valid: NDArray[np.float64], - ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: - """Return scaled X and its predictive-only reconstruction in scaled X-space.""" - Xs = self._scaled_x_validated(X_valid) - X_filtered, _ = self._filter_validated(X_valid) - T = self.pls_.transform(X_filtered) - X_hat = self._pls_x_mean() + T @ self.x_loadings_.T - return Xs, X_hat - - def _full_reconstruction_validated( - self, - X_valid: NDArray[np.float64], - ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: - """Return scaled X and full OPLS reconstruction in scaled X-space.""" - Xs = self._scaled_x_validated(X_valid) - if self.n_orthogonal_ > 0: - T_ortho = self._filter_validated(X_valid)[1] - X_ortho_hat = T_ortho @ self.x_ortho_loadings_.T - else: - X_ortho_hat = np.zeros_like(Xs) - _, X_pred_hat = self._predictive_reconstruction_validated(X_valid) - X_hat = X_ortho_hat + X_pred_hat - return Xs, X_hat - def q_residuals( self, X: ArrayLike, @@ -710,24 +677,31 @@ def _q_residuals_validated( space: str = "full", ) -> NDArray[np.float64]: """Compute Q residuals internally, assuming X is already validated.""" - if space == "full": - X_model, X_hat = self._full_reconstruction_validated(X_valid) - elif space == "predictive": - X_model = self._scaled_x_validated(X_valid) - _, X_hat = self._predictive_reconstruction_validated(X_valid) + return self._q_residuals_from_projection( + self._project_validated(X_valid), space=space + ) + + def _q_residuals_from_projection( + self, + proj: _OPLSProjection, + *, + space: str, + ) -> NDArray[np.float64]: + """Compute Q residuals from an existing OPLS projection.""" + X_pred_hat = self._predictive_x_hat(proj.t_pred) + if space == "predictive": + X_hat = X_pred_hat + elif space == "full": + if self.n_orthogonal_ > 0: + X_ortho_hat = proj.t_ortho @ self.x_ortho_loadings_.T + X_hat = X_ortho_hat + X_pred_hat + else: + X_hat = X_pred_hat else: raise ValueError("space must be one of {'full', 'predictive'}.") - resid = X_model - X_hat + resid = proj.Xs - X_hat return np.sum(resid**2, axis=1) - def _filter(self, X: ArrayLike) -> tuple[NDArray[np.float64], NDArray[np.float64]]: - """Preprocess and orthogonal-filter new ``X`` exactly as at fit time. - - Returns the filtered ``X`` and the orthogonal scores. - """ - X_valid = self._validate_X_predict(X) - return self._filter_validated(X_valid) - def __sklearn_tags__(self): tags = super().__sklearn_tags__() tags.regressor_tags.poor_score = True diff --git a/src/scikit_opls/_opls_da.py b/src/scikit_opls/_opls_da.py index 845c1ce..f859b42 100644 --- a/src/scikit_opls/_opls_da.py +++ b/src/scikit_opls/_opls_da.py @@ -30,6 +30,7 @@ from scikit_opls._opls import OPLS from scikit_opls._preprocessing import VALID_SCALING +from scikit_opls._utils import _reject_bool_param class OPLSDA(ClassifierMixin, BaseEstimator): @@ -95,10 +96,8 @@ def fit(self, X: ArrayLike, y: ArrayLike) -> OPLSDA: self : OPLSDA The fitted estimator. """ - if isinstance(self.n_components, bool): - raise ValueError("n_components must be an integer, not bool.") - if isinstance(self.n_orthogonal, bool): - raise ValueError("n_orthogonal must be an integer, not bool.") + _reject_bool_param("n_components", self.n_components) + _reject_bool_param("n_orthogonal", self.n_orthogonal) self._validate_params() @@ -172,10 +171,6 @@ def _validate_X_predict(self, X: ArrayLike) -> NDArray[np.float64]: # noqa: N80 ) return np.asarray(X_valid, dtype=np.float64) - def _validate_x_predict(self, X: ArrayLike) -> NDArray[np.float64]: - """Backward-compatible alias for the canonical validation helper.""" - return self._validate_X_predict(X) - def decision_function(self, X: ArrayLike) -> NDArray[np.float64]: """Raw signed OPLS regression output; positive favours ``classes_[1]``. @@ -191,8 +186,10 @@ def decision_function(self, X: ArrayLike) -> NDArray[np.float64]: zero are assigned to ``classes_[0]`` by :meth:`predict`. """ X_valid = self._validate_X_predict(X) - X_filtered, _ = self.opls_._filter_validated(X_valid) - return np.asarray(self.opls_.pls_.predict(X_filtered), dtype=np.float64).ravel() + proj = self.opls_._project_validated(X_valid) + return np.asarray( + self.opls_.pls_.predict(proj.X_filtered), dtype=np.float64 + ).ravel() def predict(self, X: ArrayLike) -> NDArray: """Predict class labels. diff --git a/src/scikit_opls/_orthogonal.py b/src/scikit_opls/_orthogonal.py index 62735c5..663ebc1 100644 --- a/src/scikit_opls/_orthogonal.py +++ b/src/scikit_opls/_orthogonal.py @@ -1,58 +1,31 @@ -"""OSC-style orthogonal filtering for OPLS. - -The public OPLS estimator passes a preprocessed ``X`` and a single centered -response ``y``. The low-level predictive-direction helper also accepts a -multivariate response block for tests and internal reuse, but OPLS itself is -univariate. -""" +"""OSC-style orthogonal filtering primitives for OPLS.""" from __future__ import annotations import warnings from dataclasses import dataclass -from numbers import Integral import numpy as np from numpy.typing import ArrayLike, NDArray from sklearn.exceptions import ConvergenceWarning +from scikit_opls._utils import _validate_int + _TOL = 1e-12 def _validate_n_components(n_components: int) -> int: - # ``bool`` behaves like an integer in Python, but True/False are ambiguous as - # component counts and should be treated as invalid user input. - if isinstance(n_components, bool) or not isinstance(n_components, Integral): - raise TypeError( - "n_components must be a non-negative integer, " - f"got {type(n_components).__name__}." - ) - if n_components < 0: - raise ValueError(f"n_components must be >= 0, got {n_components}.") - return int(n_components) + return _validate_int( + "n_components", n_components, minimum=0, type_phrase="a non-negative integer" + ) @dataclass class OrthogonalComponents: """Result of :func:`opls_filter`. - Attributes - ---------- - x_ortho_weights : ndarray of shape (n_features, n_components) - Orthogonal weight vectors ``w_o``. - x_ortho_scores : ndarray of shape (n_samples, n_components) - Orthogonal scores ``t_o``. - x_ortho_loadings : ndarray of shape (n_features, n_components) - Orthogonal loadings ``p_o``. - x_filtered : ndarray of shape (n_samples, n_features) - ``X`` with the orthogonal variation removed. - x_predictive_weight : ndarray of shape (n_features,) - Normalised predictive weight direction ``w_p`` when computed. For - ``opls_filter(..., n_components=0)``, this is a zero vector because the - predictive direction is unused. - n_components : int - Number of orthogonal components actually extracted (may be smaller than - requested if ``X`` ran out of orthogonal variation). + ``n_components`` may be smaller than requested if ``X`` ran out of orthogonal + variation. ``x_predictive_weight`` is a zero vector when ``n_components=0``. """ x_ortho_weights: NDArray[np.float64] @@ -64,36 +37,17 @@ class OrthogonalComponents: def predictive_weight(X: ArrayLike, Y: ArrayLike) -> NDArray[np.float64]: - """Leading joint X–Y direction (unit norm). - - Generalises ``w_p ∝ Xᵀy`` to multivariate ``Y`` via the dominant left singular - vector of ``S = Xᵀ Y``. For a single-column ``Y`` this reduces exactly to the - normalised ``Xᵀy`` (up to sign), so single-``y`` OPLS is unchanged. The public - OPLS estimator passes a single centered response; multivariate ``Y`` is - accepted here only because the predictive-direction computation is - block-agnostic. - - Parameters - ---------- - X : ndarray of shape (n_samples, n_features) - Preprocessed predictor matrix. - Y : ndarray of shape (n_samples,) or (n_samples, n_targets) - Centered response(s). - - Returns - ------- - w_p : ndarray of shape (n_features,) - Unit-normalised predictive direction. - - Raises - ------ - ValueError - If ``X`` is orthogonal to ``Y`` (the predictive direction is undefined). + """Return the unit X-side direction of maximal X/Y covariance. + + For univariate ``Y`` this is normalized ``X.T @ y``; for multivariate ``Y``, + it is the leading left singular vector of ``X.T @ Y``. """ X = np.asarray(X, dtype=np.float64) Y = np.asarray(Y, dtype=np.float64) if X.ndim != 2: raise ValueError(f"X must be 2D, got shape {X.shape}") + if X.shape[1] == 0: + raise ValueError("X must contain at least one feature.") if not np.all(np.isfinite(X)): raise ValueError("X must contain only finite values.") if not np.all(np.isfinite(Y)): @@ -111,7 +65,7 @@ def predictive_weight(X: ArrayLike, Y: ArrayLike) -> NDArray[np.float64]: w = X.T @ y norm_w = float(np.linalg.norm(w)) ref = float(np.linalg.norm(X) * np.linalg.norm(y)) - if norm_w <= _TOL * ref or norm_w == 0.0: + if norm_w <= _TOL * ref: raise ValueError( "X is numerically orthogonal to Y; predictive direction is undefined." ) @@ -130,7 +84,7 @@ def predictive_weight(X: ArrayLike, Y: ArrayLike) -> NDArray[np.float64]: S = X.T @ Y # (n_features, n_targets) s_norm = float(np.linalg.norm(S)) ref = float(np.linalg.norm(X) * np.linalg.norm(Y)) - if s_norm <= _TOL * ref or s_norm == 0.0: + if s_norm <= _TOL * ref: raise ValueError( "X is numerically orthogonal to Y; predictive direction is undefined." ) @@ -148,37 +102,17 @@ def orthogonal_filter( predictive_direction: NDArray[np.float64], n_components: int, ) -> OrthogonalComponents: - """Remove up to ``n_components`` directions in ``block`` orthogonal to a given one. - - OSC-style deflation of one block against a supplied predictive direction - (passed in rather than computed from ``y``), as used by OPLS. - - Parameters - ---------- - block : ndarray of shape (n_samples, n_features) - Preprocessed block to deflate. - predictive_direction : ndarray of shape (n_features,) - Unit-norm direction defining the predictive subspace to preserve. - n_components : int - Number of orthogonal components to extract. - - Returns - ------- - components : OrthogonalComponents - Fitted orthogonal weights/scores/loadings, the filtered block, the - predictive direction, and the number of components actually extracted. - - Notes - ----- - For each component: form the predictive score ``t = X w_p``, its loading - ``p = Xᵀt / (tᵀt)``, then the part of ``p`` orthogonal to ``w_p`` becomes the - orthogonal weight ``w_o``. The orthogonal score ``t_o = X w_o`` and loading - ``p_o = Xᵀt_o / (t_oᵀt_o)`` are deflated out: ``X <- X - t_o p_oᵀ``. + """Sequentially deflate block variation orthogonal to a predictive direction. + + The supplied predictive direction is normalized defensively. Fewer components + may be returned if no numerically resolvable orthogonal variation remains. """ n_components = _validate_n_components(n_components) X = np.asarray(block, dtype=np.float64) if X.ndim != 2: raise ValueError(f"block must be 2D, got shape {X.shape}") + if X.shape[1] == 0: + raise ValueError("block must contain at least one feature.") w_pred = np.asarray(predictive_direction, dtype=np.float64).ravel() n_samples, n_features = X.shape if w_pred.shape != (n_features,): @@ -190,17 +124,15 @@ def orthogonal_filter( if not np.all(np.isfinite(w_pred)): raise ValueError("predictive_direction must contain only finite values.") - # The deflation math below assumes a unit-norm predictive direction. Normalise it - # defensively in case a caller passes an un-normalised direction. The direction is - # unused when no components are requested, so only guard/normalise when it matters. - if n_components > 0: - w_norm = float(np.linalg.norm(w_pred)) - if w_norm <= np.finfo(np.float64).tiny: - raise ValueError( - "predictive_direction must be numerically non-zero when " - "n_components > 0." - ) + # Normalize defensively; the deflation formulas assume unit-norm w_pred. Only + # raise on a zero direction when it will actually be used (n_components > 0). + w_norm = float(np.linalg.norm(w_pred)) + if w_norm > np.finfo(np.float64).tiny: w_pred = w_pred / w_norm + elif n_components > 0: + raise ValueError( + "predictive_direction must be numerically non-zero when n_components > 0." + ) W = np.zeros((n_features, n_components)) T = np.zeros((n_samples, n_components)) @@ -217,14 +149,12 @@ def orthogonal_filter( tt = float(t @ t) if tt <= _TOL * res_norm_sq: break - # ``p`` describes how the current residual X loads onto the predictive - # score. Removing its projection onto w_pred leaves only X variation - # orthogonal to the predictive direction. p = X_res.T @ t / tt - w_o = p - float(w_pred @ p) * w_pred # part of the loading orthogonal to w_pred + # Remove the predictive-direction part of p to obtain an orthogonal weight. + w_o = p - float(w_pred @ p) * w_pred w_norm = float(np.linalg.norm(w_o)) p_norm = float(np.linalg.norm(p)) - if w_norm <= _TOL * p_norm or w_norm == 0.0: + if w_norm <= _TOL * p_norm: break # no orthogonal variation left w_o /= w_norm t_o = X_res @ w_o @@ -232,8 +162,7 @@ def orthogonal_filter( if too <= _TOL * res_norm_sq: break p_o = X_res.T @ t_o / too - # Sequentially deflate the fitted orthogonal score/loading pair so the next - # component is extracted from the remaining X variation. + # Deflate before extracting the next orthogonal component. X_res -= np.outer(t_o, p_o) W[:, i] = w_o T[:, i] = t_o @@ -260,34 +189,13 @@ def orthogonal_filter( def opls_filter(X: ArrayLike, Y: ArrayLike, n_components: int) -> OrthogonalComponents: - """OPLS X-orthogonal filter: predictive direction from ``(X, Y)``, deflate ``X``. - - Parameters - ---------- - X : ndarray of shape (n_samples, n_features) - Preprocessed (centered/scaled) predictor matrix. - Y : ndarray of shape (n_samples,) or (n_samples, n_targets) - Centered response(s). - n_components : int - Number of orthogonal components to extract. - - Returns - ------- - components : OrthogonalComponents - See :func:`orthogonal_filter`. - - Notes - ----- - The predictive direction is computed once from the original ``(X, Y)`` and - reused for every orthogonal component. For univariate ``Y`` this is exact, not a - shortcut: each orthogonal score is constructed exactly orthogonal to ``Y``, so - removing it leaves ``Xᵀy`` (hence the predictive direction ``w_p ∝ Xᵀy``) - unchanged. Recomputing ``w_p`` from each deflated residual would yield the same - direction, so the canonical Trygg-Wold OPLS algorithm coincides with this - fixed-direction filter for single-response OPLS. - - When ``n_components=0``, ``Y`` is not inspected because no predictive direction - is needed; the returned predictive weight is a zero vector. + """Compute the predictive direction from ``(X, Y)`` once, then deflate ``X``. + + Reusing one direction for every component is exact, not a shortcut: each + orthogonal score is built orthogonal to ``Y``, so removing it leaves ``Xᵀy`` + (hence the predictive direction) unchanged — recomputing it from each + deflated residual would give the same answer. When ``n_components=0``, ``Y`` + is not inspected and the returned predictive weight is a zero vector. """ n_components = _validate_n_components(n_components) X = np.asarray(X, dtype=np.float64) @@ -310,23 +218,9 @@ def apply_orthogonal_filter( x_ortho_weights: NDArray[np.float64], x_ortho_loadings: NDArray[np.float64], ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: - """Replay a fitted orthogonal filter on new data. - - Parameters - ---------- - X : ndarray of shape (n_samples, n_features) - Preprocessed predictor matrix. - x_ortho_weights : ndarray of shape (n_features, n_components) - Orthogonal weights from :func:`opls_filter`. - x_ortho_loadings : ndarray of shape (n_features, n_components) - Orthogonal loadings from :func:`opls_filter`. - - Returns - ------- - X_filtered : ndarray of shape (n_samples, n_features) - ``X`` with the fitted orthogonal variation removed. - x_ortho_scores : ndarray of shape (n_samples, n_components) - Orthogonal scores of the new samples. + """Replay fitted sequential orthogonal deflations on new preprocessed X. + + Returns ``(X_filtered, x_ortho_scores)``. """ X = np.asarray(X, dtype=np.float64) W = np.asarray(x_ortho_weights, dtype=np.float64) diff --git a/src/scikit_opls/_preprocessing.py b/src/scikit_opls/_preprocessing.py index 936db4f..7efc3e7 100644 --- a/src/scikit_opls/_preprocessing.py +++ b/src/scikit_opls/_preprocessing.py @@ -1,4 +1,4 @@ -"""Column scaling for OPLS.""" +"""Column scaling helpers for OPLS/O2PLS.""" from __future__ import annotations @@ -10,49 +10,36 @@ _EPS = np.finfo(np.float64).eps +def _as_finite_2d_array(name: str, X: ArrayLike) -> NDArray[np.float64]: + """Return ``X`` as a finite non-empty 2D float64 array.""" + arr = np.asarray(X, dtype=np.float64) + if arr.ndim != 2: + raise ValueError(f"{name} must be 2D, got shape {arr.shape}.") + if arr.shape[0] == 0 or arr.shape[1] == 0: + raise ValueError(f"{name} must have at least one sample and one feature.") + if not np.all(np.isfinite(arr)): + raise ValueError(f"{name} must contain only finite values.") + return arr + + def compute_scaling( X: ArrayLike, mode: str ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: - """Return ``(mean_, scale_)`` column vectors for the requested scaling mode. - - - ``none``: no centering, no scaling. - - ``center``: mean-centering only. - - ``pareto``: mean-centering, divide by ``sqrt(std)``. - - ``standard``: mean-centering, divide by ``std`` (unit variance). - - Standard deviation uses ``ddof=1`` (sample). Columns with - (near) zero variance get a scale of ``1.0`` to avoid division by zero. - - Parameters - ---------- - X : ndarray of shape (n_samples, n_features) - Predictor matrix. - mode : {"none", "center", "pareto", "standard"} - Scaling mode. - - Returns - ------- - mean_ : ndarray of shape (n_features,) - Per-column centering vector. - scale_ : ndarray of shape (n_features,) - Per-column scaling vector. + """Return column means and scales for a supported preprocessing mode. + + Modes are ``"none"``, ``"center"``, ``"pareto"``, and ``"standard"``. + Sample standard deviation uses ``ddof=1``; near-constant columns receive + scale ``1.0`` to preserve feature alignment. """ if mode not in VALID_SCALING: raise ValueError( f"Unknown scaling mode {mode!r}; expected one of {VALID_SCALING}." ) - X = np.asarray(X, dtype=np.float64) - if X.ndim != 2: - raise ValueError(f"X must be 2D, got shape {X.shape}.") - if X.shape[0] == 0 or X.shape[1] == 0: - raise ValueError("X must have at least one sample and one feature.") - if not np.all(np.isfinite(X)): - raise ValueError("X must contain only finite values.") + X = _as_finite_2d_array("X", X) n_features = X.shape[1] if mode == "none": - # Identity preprocessing is represented as zero centering and unit scaling - # so apply_scaling can use one code path for every mode. + # Represent identity preprocessing with the same mean/scale contract. return ( np.zeros(n_features, dtype=np.float64), np.ones(n_features, dtype=np.float64), @@ -65,56 +52,34 @@ def compute_scaling( if X.shape[0] > 1: std = X.std(axis=0, ddof=1) else: - std = np.ones(n_features) - # Keep constant columns aligned with the original feature matrix. They remain - # constant after scaling instead of causing division by zero. + std = np.ones(n_features, dtype=np.float64) + + # Preserve constant columns instead of dropping or dividing by zero. std = np.where(std <= _EPS, 1.0, std) if mode == "pareto": - # Pareto scaling divides by sqrt(sample std), reducing large-variance - # features without forcing every feature to unit variance. + # Pareto scaling divides by sqrt(sample std). return mean_, np.sqrt(std) - return mean_, std # standard + return mean_, std def apply_scaling( X: ArrayLike, mean_: ArrayLike, scale_: ArrayLike ) -> NDArray[np.float64]: - """Apply a previously computed centering/scaling to ``X``. - - Parameters - ---------- - X : ndarray of shape (n_samples, n_features) - Predictor matrix. - mean_ : ndarray of shape (n_features,) - Centering vector from :func:`compute_scaling`. - scale_ : ndarray of shape (n_features,) - Scaling vector from :func:`compute_scaling`. - - Returns - ------- - X_scaled : ndarray of shape (n_samples, n_features) - ``(X - mean_) / scale_``. - """ - X = np.asarray(X, dtype=np.float64) + """Apply fitted column centering/scaling as ``(X - mean_) / scale_``.""" + X = _as_finite_2d_array("X", X) mean_ = np.asarray(mean_, dtype=np.float64) scale_ = np.asarray(scale_, dtype=np.float64) - if X.ndim != 2: - raise ValueError(f"X must be 2D, got shape {X.shape}.") - if X.shape[0] == 0 or X.shape[1] == 0: - raise ValueError("X must have at least one sample and one feature.") if mean_.shape != (X.shape[1],): raise ValueError(f"mean_ must have shape ({X.shape[1]},), got {mean_.shape}.") if scale_.shape != (X.shape[1],): raise ValueError(f"scale_ must have shape ({X.shape[1]},), got {scale_.shape}.") - if not np.all(np.isfinite(X)): - raise ValueError("X must contain only finite values.") if not np.all(np.isfinite(mean_)): raise ValueError("mean_ must contain only finite values.") if not np.all(np.isfinite(scale_)): raise ValueError("scale_ must contain only finite values.") - if np.any(scale_ == 0.0): - raise ValueError("scale_ must not contain zeros.") + if np.any(scale_ <= 0.0): + raise ValueError("scale_ must contain only positive values.") return (X - mean_) / scale_ diff --git a/src/scikit_opls/_utils.py b/src/scikit_opls/_utils.py index 12cab21..d6a234c 100644 --- a/src/scikit_opls/_utils.py +++ b/src/scikit_opls/_utils.py @@ -2,12 +2,44 @@ from __future__ import annotations +from numbers import Integral + import numpy as np from numpy.typing import ArrayLike _EPS = np.finfo(np.float64).eps +def _reject_bool_param(name: str, value: object) -> None: + """Raise if ``value`` is ``bool`` for a sklearn Integral-constrained parameter. + + ``bool`` is a subclass of ``int``, so sklearn's ``Interval(Integral, ...)`` + constraint accepts ``True``/``False`` silently. Call this before + ``self._validate_params()`` for each Integral-constrained hyperparameter. + """ + if isinstance(value, bool): + raise ValueError(f"{name} must be an integer, not bool.") + + +def _validate_int( + name: str, + value: object, + *, + minimum: int, + type_phrase: str = "an integer", +) -> int: + """Reject non-``Integral``, ``bool``, and below-``minimum`` values. + + ``bool`` is a subclass of ``int``; reject it explicitly so ``True``/``False`` + are not silently accepted as counts or indices. + """ + if isinstance(value, bool) or not isinstance(value, Integral): + raise TypeError(f"{name} must be {type_phrase}, got {type(value).__name__}.") + if value < minimum: + raise ValueError(f"{name} must be >= {minimum}, got {value}.") + return int(value) + + def _has_nonzero_variation( values: ArrayLike, *, diff --git a/src/scikit_opls/plotting.py b/src/scikit_opls/plotting.py index 41be736..497d435 100644 --- a/src/scikit_opls/plotting.py +++ b/src/scikit_opls/plotting.py @@ -15,7 +15,6 @@ class with a :meth:`from_estimator` constructor that computes the plotted arrays from __future__ import annotations import warnings -from numbers import Integral from typing import TYPE_CHECKING import numpy as np @@ -28,6 +27,7 @@ class with a :meth:`from_estimator` constructor that computes the plotted arrays from scikit_opls._opls import OPLS from scikit_opls._opls_da import OPLSDA from scikit_opls._preprocessing import apply_scaling +from scikit_opls._utils import _validate_int if TYPE_CHECKING: import matplotlib.axes @@ -52,11 +52,7 @@ def _validate_component_index(value: object, name: str) -> int: Rejects booleans (a subclass of ``int``) and non-integers so a stray float or ``True`` cannot silently index the wrong column or fail later inside NumPy. """ - if isinstance(value, bool) or not isinstance(value, Integral): - raise TypeError(f"{name} must be an integer index, got {type(value).__name__}.") - if value < 0: - raise ValueError(f"{name} must be >= 0, got {value}.") - return int(value) + return _validate_int(name, value, minimum=0, type_phrase="an integer index") def _unwrap_estimator_and_data( @@ -64,16 +60,9 @@ def _unwrap_estimator_and_data( ) -> tuple[OPLS, NDArray[np.float64]]: """Unwrap wrappers, returning fitted base model and the X seen by that model. - Accepts ``OPLS``/``OPLSDA``, a pipeline ending in one, or a fitted search - meta-estimator exposing ``best_estimator_`` around either shape. + Accepts ``OPLS``/``OPLSDA``, or a pipeline ending in one. """ - if hasattr(estimator, "cv_results_") and not hasattr(estimator, "best_estimator_"): - raise TypeError( - "Search meta-estimators must be fitted with refit=True so plotting " - "can use best_estimator_." - ) - # Search estimators keep the selected pipeline/model on best_estimator_. - inner = getattr(estimator, "best_estimator_", estimator) + inner = estimator if isinstance(inner, Pipeline): check_is_fitted(inner) @@ -107,10 +96,8 @@ def _unwrap_estimator_and_data( base = inner else: raise TypeError( - "estimator must be a fitted OPLS, OPLSDA, Pipeline ending in one, " - "or a fitted search meta-estimator wrapping one. Ensemble and " - "meta-classifier wrappers are unsupported because they do not expose " - "one latent OPLS space." + "estimator must be a fitted OPLS, OPLSDA, or Pipeline ending in one. " + "Ensemble, search estimators, and meta-classifier wrappers are unsupported." ) if not isinstance(base, OPLS): @@ -133,9 +120,8 @@ def _unwrap_estimator_and_data( class OPLSScoresDisplay: """Predictive vs orthogonal score scatter for an OPLS-family model. - Works for :class:`~scikit_opls.OPLS`, :class:`~scikit_opls.OPLSDA`, a pipeline - ending in one, or a fitted search meta-estimator exposing ``best_estimator_`` - around either shape. Construct with :meth:`from_estimator`. + Works for :class:`~scikit_opls.OPLS`, :class:`~scikit_opls.OPLSDA`, or a + pipeline ending in one. Construct with :meth:`from_estimator`. Parameters ---------- @@ -198,7 +184,7 @@ def from_estimator( Parameters ---------- - estimator : OPLS, OPLSDA, Pipeline or search meta-estimator + estimator : OPLS, OPLSDA or Pipeline A fitted estimator. When passing a pipeline, pass raw ``X`` as expected by that pipeline. When passing the final OPLS-family step directly, pass the already transformed matrix seen by that step. @@ -253,14 +239,10 @@ def from_estimator( f"estimator with {n_ortho} orthogonal component(s)." ) - # Project supplied data through the fitted filter before asking the PLS - # engine for predictive scores. - X_filtered, t_ortho = base._filter_validated(X_trans) - scores = base.pls_.transform(X_filtered) - if isinstance(scores, tuple): - t_pred_arr = scores[0] - else: - t_pred_arr = scores + # Project supplied data through the fitted preprocessing/filtering path. + proj = base._project_validated(X_trans) + t_pred_arr = proj.t_pred + t_ortho = proj.t_ortho t_pred = t_pred_arr[:, predictive_component] t_o = t_ortho[:, orthogonal_component] if n_ortho > 0 else np.zeros_like(t_pred) display = cls( @@ -329,9 +311,8 @@ class SPlotDisplay: Points may therefore represent transformed features rather than original input columns. - Accepts :class:`~scikit_opls.OPLS`, :class:`~scikit_opls.OPLSDA`, a pipeline - ending in one, or a fitted search meta-estimator exposing ``best_estimator_`` - around either shape. Construct with :meth:`from_estimator`. + Accepts :class:`~scikit_opls.OPLS`, :class:`~scikit_opls.OPLSDA`, or a pipeline + ending in one. Construct with :meth:`from_estimator`. Parameters ---------- @@ -379,7 +360,7 @@ def from_estimator( Parameters ---------- - estimator : OPLS, OPLSDA, Pipeline or search meta-estimator + estimator : OPLS, OPLSDA or Pipeline A fitted estimator. When passing a pipeline, pass raw ``X`` as expected by that pipeline. When passing the final OPLS-family step directly, pass the already transformed matrix seen by that step. @@ -433,13 +414,8 @@ def from_estimator( # Use the fitted predictive score for the selected component as the common # reference vector for both covariance and correlation. - X_filtered, _ = base._filter_validated(X_trans) - scores = base.pls_.transform(X_filtered) - if isinstance(scores, tuple): - t_arr = scores[0] - else: - t_arr = scores - t = np.asarray(t_arr)[:, component] + proj = base._project_validated(X_trans) + t = np.asarray(proj.t_pred)[:, component] t = t - t.mean() n = t.shape[0] diff --git a/src/scikit_opls/validation.py b/src/scikit_opls/validation.py index 3807359..1da9be0 100644 --- a/src/scikit_opls/validation.py +++ b/src/scikit_opls/validation.py @@ -1,8 +1,7 @@ """Permutation testing for OPLS model significance.""" -# sklearn.base.clone, check_array and joblib.Parallel are under-typed (Parallel is -# annotated as returning Optional); suppress the resulting static-checker false -# positives (the test suite is the real correctness gate). +# sklearn/joblib typing is incomplete for clone, check_array and Parallel. +# Runtime validation and tests are the correctness gate. # pyright: reportAttributeAccessIssue=false, reportArgumentType=false # pyright: reportGeneralTypeIssues=false @@ -24,32 +23,43 @@ cross_val_predict, ) from sklearn.utils import check_random_state -from sklearn.utils.validation import check_array, check_consistent_length +from sklearn.utils.validation import check_array, check_consistent_length, column_or_1d -from scikit_opls._utils import _has_nonzero_variation +from scikit_opls._utils import _has_nonzero_variation, _validate_int _CVType = int | BaseCrossValidator | BaseShuffleSplit | Iterable | None +def _as_univariate_array(name: str, values: ArrayLike) -> NDArray[np.float64]: + """Return values as a finite 1D float64 array.""" + try: + arr = column_or_1d(np.asarray(values, dtype=np.float64), warn=False) + except ValueError as exc: + raise ValueError( + f"{name} must be univariate; multi-output targets are not supported." + ) from exc + if not np.all(np.isfinite(arr)): + raise ValueError(f"{name} must contain only finite values.") + return arr + + def _safe_r2_score(y_true: ArrayLike, y_pred: ArrayLike) -> float: - y_true_arr = np.asarray(y_true, dtype=np.float64).ravel() - y_pred_arr = np.asarray(y_pred, dtype=np.float64).ravel() + y_true_arr = _as_univariate_array("y_true", y_true) + y_pred_arr = _as_univariate_array("y_pred", y_pred) if y_true_arr.shape != y_pred_arr.shape: raise ValueError( "y_true and y_pred must have the same flattened shape, " f"got {y_true_arr.shape} and {y_pred_arr.shape}." ) if not _has_nonzero_variation(y_true_arr): - # sklearn's r2_score defines constant-target cases awkwardly for model - # significance; NaN makes the undefined metric explicit downstream. - return np.nan + return float("nan") return float(r2_score(y_true_arr, y_pred_arr)) def _cross_val_q2( estimator: BaseEstimator, X: ArrayLike, y: ArrayLike, cv: _CVType ) -> float: - """Out-of-fold Q2 of ``estimator`` on ``(X, y)`` using the provided ``cv``.""" + """Out-of-fold Q2 of ``estimator`` on ``(X, y)``.""" y_pred = cross_val_predict(clone(estimator), X, y, cv=cv) return _safe_r2_score(y, y_pred) @@ -77,27 +87,29 @@ class PermutationResult: def _fitted_r2y(fitted: BaseEstimator) -> float: - # GridSearchCV and similar search estimators expose the selected model through - # best_estimator_; recurse until we reach the OPLS-like estimator itself. if hasattr(fitted, "r2y_"): return float(getattr(fitted, "r2y_")) - if hasattr(fitted, "cv_results_") and not hasattr(fitted, "best_estimator_"): + + best = getattr(fitted, "best_estimator_", None) + if best is not None: + return _fitted_r2y(best) + + if hasattr(fitted, "cv_results_"): raise TypeError( "Search meta-estimators must use refit=True so permutation_test can " "access best_estimator_." ) - if hasattr(fitted, "best_estimator_"): - return _fitted_r2y(getattr(fitted, "best_estimator_")) + raise TypeError( - "permutation_test requires an OPLS-like regression estimator exposing r2y_, " - "or a GridSearchCV wrapping one." + "permutation_test requires an OPLS-like regression estimator exposing " + "r2y_, or a refit-enabled search estimator wrapping one." ) def _permuted_scores( estimator: BaseEstimator, X: ArrayLike, y_perm: ArrayLike, cv: _CVType ) -> tuple[float, float]: - """R2Y and out-of-fold Q2 for one permuted target (one parallel task).""" + """Return R2Y/Q2 for one permuted target.""" fitted = clone(estimator).fit(X, y_perm) r2y = _fitted_r2y(fitted) q2 = _cross_val_q2(estimator, X, y_perm, cv=cv) @@ -105,12 +117,27 @@ def _permuted_scores( def _contains_classifier(estimator: BaseEstimator) -> bool: - # Walk simple meta-estimators such as CalibratedClassifierCV(estimator=...). - if is_classifier(estimator): - return True - if hasattr(estimator, "estimator"): - return _contains_classifier(getattr(estimator, "estimator")) - return False + """Return whether estimator is tagged as a classifier.""" + return is_classifier(estimator) + + +def _resolve_cv(estimator: BaseEstimator, cv: _CVType, y: NDArray[np.float64]): + if cv is None: + estimator_cv = getattr(estimator, "cv", None) + cv = estimator_cv if estimator_cv is not None else min(5, len(y)) + + # Materialize one-shot split iterables so observed and permuted passes reuse + # the same splits instead of consuming the iterator once. + if cv is not None and not isinstance(cv, Integral) and not hasattr(cv, "split"): + cv = list(cv) + + return check_cv(cv, y=y, classifier=False) + + +def _empirical_p_value(observed: float, permuted: NDArray[np.float64]) -> float: + if np.isnan(observed): + return float("nan") + return float((1 + int(np.sum(permuted >= observed))) / (permuted.size + 1)) def permutation_test( @@ -136,8 +163,8 @@ def permutation_test( An unfitted OPLS-like estimator (cloned internally for each fit). X : array-like of shape (n_samples, n_features) Predictors. - y : array-like of shape (n_samples,) - Response. + y : array-like of shape (n_samples,) or (n_samples, 1) + Univariate response. Multi-output targets are rejected. n_permutations : int, default=20 Number of label permutations. cv : int, cross-validation generator or None, default=None @@ -172,47 +199,24 @@ def permutation_test( "permutation_test is for regression models; " "classifiers like OPLSDA are not supported." ) - if isinstance(n_permutations, bool) or not isinstance(n_permutations, Integral): - raise TypeError( - f"n_permutations must be an integer, got {type(n_permutations).__name__}" - ) - if n_permutations < 1: - raise ValueError(f"n_permutations must be >= 1, got {n_permutations}") + n_permutations = _validate_int("n_permutations", n_permutations, minimum=1) X = check_array(X, dtype=np.float64) - y = np.asarray(y, dtype=np.float64).ravel() + y = _as_univariate_array("y", y) check_consistent_length(X, y) - if not np.all(np.isfinite(y)): - raise ValueError("y must contain only finite values.") if len(y) < 3: raise ValueError( "permutation_test requires at least 3 samples so each CV training " "fold can contain at least 2 samples." ) - if cv is None: - estimator_cv = getattr(estimator, "cv", None) - # Prefer an estimator-owned cv setting when present; otherwise keep folds - # valid for small data by capping the default at n_samples. - cv = estimator_cv if estimator_cv is not None else min(5, len(y)) - # A one-shot iterable of splits would be consumed by the observed-Q2 pass and - # leave nothing for the permutations; materialise it so every pass sees the - # same splits. - if cv is not None and not isinstance(cv, Integral) and not hasattr(cv, "split"): - cv = list(cv) - cv_checked = check_cv(cv, y=y, classifier=False) + cv_checked = _resolve_cv(estimator, cv, y) - # Fit once on the true labels to establish the observed in-sample R2Y. fitted = clone(estimator).fit(X, y) observed_r2y = _fitted_r2y(fitted) - - rng = check_random_state(random_state) - # Q2 is always out-of-fold, so compute it through the same CV object used for - # every permutation. observed_q2 = _cross_val_q2(estimator, X, y, cv=cv_checked) - # Draw all permutations serially from the RNG so the result is independent of - # the execution order the parallel backend chooses. + rng = check_random_state(random_state) perms = [rng.permutation(y) for _ in range(n_permutations)] scored = Parallel(n_jobs=n_jobs)( delayed(_permuted_scores)(estimator, X, y_perm, cv_checked) for y_perm in perms @@ -220,22 +224,11 @@ def permutation_test( permuted_r2y = np.asarray([r2y for r2y, _ in scored], dtype=np.float64) permuted_q2 = np.asarray([q2 for _, q2 in scored], dtype=np.float64) - # An undefined observed metric (NaN) must not masquerade as significant. - r2y_p = ( - np.nan - if np.isnan(observed_r2y) - else (1 + int(np.sum(permuted_r2y >= observed_r2y))) / (n_permutations + 1) - ) - q2_p = ( - np.nan - if np.isnan(observed_q2) - else (1 + int(np.sum(permuted_q2 >= observed_q2))) / (n_permutations + 1) - ) return PermutationResult( r2y=observed_r2y, q2=observed_q2, permuted_r2y=permuted_r2y, permuted_q2=permuted_q2, - r2y_p_value=float(r2y_p), - q2_p_value=float(q2_p), + r2y_p_value=_empirical_p_value(observed_r2y, permuted_r2y), + q2_p_value=_empirical_p_value(observed_q2, permuted_q2), ) diff --git a/tests/test_diagnostics.py b/tests/test_diagnostics.py index c686b4f..0b65efb 100644 --- a/tests/test_diagnostics.py +++ b/tests/test_diagnostics.py @@ -257,7 +257,7 @@ def test_pls_engine_exposes_internal_x_mean_for_reconstruction(): X, y = _regression_data() model = OPLS(n_components=1, n_orthogonal=0, scale="none").fit(X, y) assert hasattr(model.pls_, "_x_mean") - assert model.pls_._x_mean.shape == (X.shape[1],) + assert getattr(model.pls_, "_x_mean").shape == (X.shape[1],) # ============================================================================== @@ -335,22 +335,15 @@ def test_oplsda_diagnostics_expect_raw_x_not_prescaled_x(): ) -def test_component_r2_from_cumulative_empty_and_differences(): - """Verify component_r2_from_cumulative with empty and non-empty inputs.""" - from scikit_opls._inspection import component_r2_from_cumulative +def test_component_r2y_from_scores_1d_y_loadings_matches_2d(): + """A 1D y_loadings is one value per component, matching the 2D (1, n) form.""" + from scikit_opls._inspection import component_r2y_from_scores - assert component_r2_from_cumulative(np.array([])).shape == (0,) - np.testing.assert_allclose( - component_r2_from_cumulative(np.array([0.2, 0.5, 0.8])), - [0.2, 0.3, 0.3], - ) - - -def test_cumulative_r2_from_residuals_decreases_with_smaller_residuals(): - """Verify cumulative_r2_from_residuals decreases as residuals decrease.""" - from scikit_opls._inspection import cumulative_r2_from_residuals + rng = np.random.default_rng(0) + y = rng.normal(size=10) + T = rng.normal(size=(10, 3)) + q_1d = rng.normal(size=3) - X = np.eye(4) - out = cumulative_r2_from_residuals(X, [0.5 * X, 0.1 * X]) - assert out[1] > out[0] - assert np.all((0.0 <= out) & (out <= 1.0)) + out_1d = component_r2y_from_scores(y, T, q_1d) + out_2d = component_r2y_from_scores(y, T, q_1d.reshape(1, -1)) + np.testing.assert_allclose(out_1d, out_2d) diff --git a/tests/test_o2pls.py b/tests/test_o2pls.py index fa54ffd..bc34abf 100644 --- a/tests/test_o2pls.py +++ b/tests/test_o2pls.py @@ -125,6 +125,13 @@ def test_score_matches_r2_score_for_1d_y(): assert model.score(X, y) == pytest.approx(r2_score(y, model.predict(X))) +def test_score_matches_r2_score_uniform_average_for_multi_output_y(): + X, Y = _make_o2pls_data(seed=7) + model = O2PLS(n_components=2).fit(X, Y) + + assert model.score(X, Y) == pytest.approx(r2_score(Y, model.predict(X))) + + @pytest.mark.parametrize( ("method", "arg_name"), [ @@ -246,6 +253,20 @@ def test_o2pls_predict_unscales_y(): assert_allclose(model.predict(X_raw), expected) +@pytest.mark.parametrize("scale", ["none", "center", "pareto", "standard"]) +def test_o2pls_predict_matches_scaled_filtered_coefficient_path_all_scales(scale): + X, y = _regression_data() + rng = np.random.default_rng(0) + Y = np.column_stack([y, y + 0.1 * rng.normal(size=y.shape)]) + X_raw = 10.0 + 3.0 * X + Y_raw = -5.0 + 2.0 * Y + model = O2PLS(scale=scale).fit(X_raw, Y_raw) + expected = ( + model.filter_transform_x(X_raw) @ model.coef_filtered_ + ) * model.y_std_ + model.y_mean_ + assert_allclose(model.predict(X_raw), expected) + + def test_o2pls_predict_x_unscales_x(): X, y = _regression_data() rng = np.random.default_rng(0) diff --git a/tests/test_opls_coefficients.py b/tests/test_opls_coefficients.py index f9f6c41..bd985da 100644 --- a/tests/test_opls_coefficients.py +++ b/tests/test_opls_coefficients.py @@ -65,6 +65,21 @@ def test_raw_coefficients_shape_and_no_coef_alias(): assert not hasattr(model, "coef_") +def test_raw_coefficients_reproduce_predict_scale_none(): + """Regression guard: scale="none" skips the standard-scaling centering path, + so raw-coefficient composition must still match predict() exactly.""" + rng = np.random.default_rng(6) + X = rng.normal(size=(80, 8)) + beta = rng.normal(size=8) + y = X @ beta + 0.1 * rng.normal(size=80) + + model = OPLS(scale="none", n_components=1, n_orthogonal=1).fit(X, y) + + y_predict = model.predict(X) + y_linear = (X @ model.coef_raw_.T + model.intercept_raw_).ravel() + assert_allclose(y_predict, y_linear, rtol=1e-10, atol=1e-10) + + def test_oplsda_inner_opls_has_raw_coefficients(): rng = np.random.default_rng(3) X = rng.normal(size=(40, 6)) diff --git a/tests/test_plotting.py b/tests/test_plotting.py index a1fe121..02505d7 100644 --- a/tests/test_plotting.py +++ b/tests/test_plotting.py @@ -83,14 +83,18 @@ def test_scores_display_replots_on_given_axes(): plt.close("all") -def test_splot_display_from_estimator_with_cv(): +def test_splot_display_requires_best_estimator_for_cv(): from matplotlib.collections import PathCollection X, y = _regression_data() model = GridSearchCV( OPLS(n_components=1), {"n_orthogonal": [0, 1, 2, 3]}, cv=4 ).fit(X, y) - disp = SPlotDisplay.from_estimator(model, X) + + with pytest.raises(TypeError, match="search estimators"): + SPlotDisplay.from_estimator(model, X) + + disp = SPlotDisplay.from_estimator(model.best_estimator_, X) assert isinstance(disp, SPlotDisplay) assert disp.covariance.shape == (X.shape[1],) assert disp.correlation.shape == (X.shape[1],) @@ -234,7 +238,7 @@ def test_splot_pipeline_feature_space_matches_transformed_features(): plt.close("all") -def test_scores_display_search_wrapper_around_pipeline(): +def test_scores_display_requires_best_estimator_for_search_pipeline(): X, y = _regression_data() search = GridSearchCV( Pipeline( @@ -247,7 +251,10 @@ def test_scores_display_search_wrapper_around_pipeline(): cv=3, ).fit(X, y) - disp = OPLSScoresDisplay.from_estimator(search, X) + with pytest.raises(TypeError, match="search estimators"): + OPLSScoresDisplay.from_estimator(search, X) + + disp = OPLSScoresDisplay.from_estimator(search.best_estimator_, X) assert isinstance(disp, OPLSScoresDisplay) assert disp.t_predictive.shape == (X.shape[0],) @@ -302,7 +309,7 @@ def test_plotting_rejects_unsupported_pipeline_shapes(): SPlotDisplay.from_estimator(non_opls_final, X) -def test_plotting_search_refit_false_raises_clear_error(): +def test_plotting_search_wrapper_raises_clear_error(): X, y = _regression_data() search = GridSearchCV( OPLS(n_components=1), @@ -311,7 +318,7 @@ def test_plotting_search_refit_false_raises_clear_error(): refit=False, ).fit(X, y) - with pytest.raises(TypeError, match="refit=True"): + with pytest.raises(TypeError, match="search estimators"): SPlotDisplay.from_estimator(search, X) diff --git a/tests/test_preprocessing.py b/tests/test_preprocessing.py index 5162f73..2ae5aa4 100644 --- a/tests/test_preprocessing.py +++ b/tests/test_preprocessing.py @@ -79,10 +79,7 @@ def test_apply_scaling_validates_representative_bad_inputs(): # 2. wrong mean_ shape with pytest.raises(ValueError, match="mean_ must have shape"): apply_scaling(X, np.zeros(2), np.ones(3)) - # 3. zero scale - with pytest.raises(ValueError, match="scale_ must not contain zeros"): - apply_scaling(X, np.zeros(3), np.array([1.0, 0.0, 1.0])) - # 4. nonfinite input + # 3. nonfinite input with pytest.raises(ValueError, match="finite"): apply_scaling(np.array([[1.0, np.inf, 1.0]]), np.zeros(3), np.ones(3)) @@ -153,3 +150,25 @@ def test_apply_scaling_extra_validation(): apply_scaling(X, np.array([0.0, np.nan, 0.0]), np.ones(3)) with pytest.raises(ValueError, match="finite"): apply_scaling(X, np.zeros(3), np.array([1.0, np.inf, 1.0])) + + +def test_apply_scaling_rejects_negative_scale(): + X = np.ones((3, 2)) + mean = np.zeros(2) + scale = np.array([1.0, -1.0]) + + with pytest.raises(ValueError, match="positive"): + apply_scaling(X, mean, scale) + + +@pytest.mark.parametrize( + "scale", + [ + np.array([1.0, 0.0, 1.0]), + np.array([1.0, -1.0, 1.0]), + ], +) +def test_apply_scaling_rejects_non_positive_scale(scale): + X = np.ones((4, 3)) + with pytest.raises(ValueError, match="scale_ must contain only positive values"): + apply_scaling(X, np.zeros(3), scale) diff --git a/tests/test_validation.py b/tests/test_validation.py index c0fa7a8..d920fc9 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -6,8 +6,12 @@ import pytest from sklearn.utils._testing import assert_allclose -from scikit_opls import OPLS -from scikit_opls.validation import _safe_r2_score, permutation_test +from scikit_opls import OPLS, OPLSDA +from scikit_opls.validation import ( + _contains_classifier, + _safe_r2_score, + permutation_test, +) from ._data import make_regression_data as _regression_data @@ -61,6 +65,25 @@ def test_permutation_test_mismatched_lengths_raise(): permutation_test(OPLS(n_orthogonal=1), X, y[:-1]) +def test_permutation_test_accepts_single_column_y(): + X, y = _regression_data(seed=8) + result = permutation_test( + OPLS(n_orthogonal=1), + X, + y.reshape(-1, 1), + n_permutations=5, + random_state=0, + ) + assert result.r2y is not None + + +def test_permutation_test_rejects_multioutput_y(): + X, y = _regression_data(seed=8) + Y = np.column_stack([y, y]) + with pytest.raises(ValueError, match="univariate response|multi-output"): + permutation_test(OPLS(n_orthogonal=1), X, Y, n_permutations=5, random_state=0) + + def test_permutation_test_n_jobs_is_reproducible(): """Parallel execution must match serial: permutations are drawn up front.""" X, y = _regression_data(seed=9) @@ -72,8 +95,6 @@ def test_permutation_test_n_jobs_is_reproducible(): def test_permutation_test_non_regression_estimator_raises(): - from scikit_opls import OPLSDA - X, y = _regression_data(seed=10) labels = np.where(y > 0.0, "hi", "lo") # OPLSDA is a classifier, raising a clean TypeError immediately @@ -205,8 +226,6 @@ def test_permutation_test_classifier_wrapped_raises(): from sklearn.model_selection import GridSearchCV from sklearn.pipeline import Pipeline - from scikit_opls import OPLSDA - X, y = _regression_data(seed=42) # 1. Pipeline containing a classifier pipe = Pipeline([("clf", OPLSDA())]) @@ -217,3 +236,24 @@ def test_permutation_test_classifier_wrapped_raises(): gs = GridSearchCV(pipe, {"clf__n_components": [1]}) with pytest.raises(TypeError, match="classifiers like OPLSDA are not supported"): permutation_test(gs, X, y) + + +def test_contains_classifier_detects_oplsda(): + assert _contains_classifier(OPLSDA()) + assert not _contains_classifier(OPLS()) + + +def test_contains_classifier_detects_pipeline_ending_in_oplsda(): + from sklearn.pipeline import Pipeline + + assert _contains_classifier(Pipeline([("model", OPLSDA())])) + assert not _contains_classifier(Pipeline([("model", OPLS())])) + + +def test_contains_classifier_detects_grid_search_over_oplsda(): + from sklearn.model_selection import GridSearchCV + + search = GridSearchCV(OPLSDA(), {"n_orthogonal": [0, 1]}, cv=2) + assert _contains_classifier(search) + other = GridSearchCV(OPLS(), {"n_orthogonal": [0, 1]}, cv=2) + assert not _contains_classifier(other)