Simplify

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.

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.

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:

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