From eea90ea034edd6ad16fb42688a02ac9ff7f46960 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 19:13:56 +0000 Subject: [PATCH 1/7] Add PLS.invert() for model inversion and null-space design Given a desired response on the original Y scale, PLS.invert() runs the model backwards to return the minimum-norm input vector that achieves it, plus an orthonormal basis for the (A - r)-dimensional null space of input vectors that all yield the same prediction. Callers can move along that space (null_space_coordinates) to meet secondary criteria without changing the predicted response, and the solution's Hotelling's T2 is reported to flag extrapolation. Reproduces the null-space results of Garcia-Carrion et al. (2025, J. Chemometrics), whose single-response proof identifies this null space with the orthogonal space of an O-PLS model. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DsTypY4FcRbrXRxorb9raz --- CHANGELOG.md | 20 ++- CITATION.cff | 2 +- pyproject.toml | 2 +- src/process_improve/multivariate/_pls.py | 167 +++++++++++++++++++++++ 4 files changed, 188 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a4a9826..a4a36d4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,23 @@ those changes. ## [Unreleased] +## [1.59.0] - 2026-07-23 + +### Added + +- `PLS.invert()`: PLS model inversion for latent-variable product and process + design. Given a desired response on the original Y scale, it returns the + minimum-norm (direct-inversion) input vector that the model predicts will + achieve it, together with an orthonormal basis for the *null space* - the + `(A - r)`-dimensional family of input vectors that all yield the same + prediction (`A` components, response rank `r`). Callers can pass + `null_space_coordinates` to move along that space and satisfy secondary + criteria (cost, safety, operability) without changing the predicted response. + The result also reports the solution's Hotelling's T2 to flag extrapolation + beyond the calibration data. Reproduces the null-space results of + García-Carrión et al. (2025), whose single-response proof identifies this null + space with the orthogonal space of an O-PLS model. + ## [1.58.0] - 2026-07-23 ### Added @@ -2599,7 +2616,8 @@ this entry records them together. - Reworked the README with a sharper value proposition and a "Why not scikit-learn?" comparison table. -[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.58.0...HEAD +[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.59.0...HEAD +[1.59.0]: https://github.com/kgdunn/process-improve/compare/v1.58.0...v1.59.0 [1.58.0]: https://github.com/kgdunn/process-improve/compare/v1.57.0...v1.58.0 [1.57.0]: https://github.com/kgdunn/process-improve/compare/v1.56.0...v1.57.0 [1.56.0]: https://github.com/kgdunn/process-improve/compare/v1.55.1...v1.56.0 diff --git a/CITATION.cff b/CITATION.cff index a544588..fca9fa5 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -12,7 +12,7 @@ authors: repository-code: "https://github.com/kgdunn/process-improve" url: "https://kgdunn.github.io/process-improve/" license: MIT -version: 1.58.0 +version: 1.59.0 date-released: "2026-07-23" keywords: - chemometrics diff --git a/pyproject.toml b/pyproject.toml index 3996999..7fc5fc3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.58.0" +version = "1.59.0" description = 'Designed Experiments; Latent Variables (PCA, PLS, multivariate methods with missing data); Process Monitoring; Batch data analysis.' readme = "README.md" license = "MIT" diff --git a/src/process_improve/multivariate/_pls.py b/src/process_improve/multivariate/_pls.py index ac14d85..f6045d4 100644 --- a/src/process_improve/multivariate/_pls.py +++ b/src/process_improve/multivariate/_pls.py @@ -991,6 +991,173 @@ def diagnose(self, X: DataMatrix) -> Bunch: return Bunch(scores=scores, hotellings_t2=t2, spe=spe_values, y_hat=y_hat) + def invert( + self, + y_desired: float | np.ndarray | pd.Series | pd.DataFrame | dict, + *, + null_space_coordinates: np.ndarray | None = None, + ) -> Bunch: + r"""Invert the PLS model: find inputs that yield a desired response. + + PLS is normally used in the forward direction (:meth:`predict`): given + inputs ``X``, predict the response ``Y``. Model *inversion* runs the + model backwards: fix the response you want (``y_desired``) and solve for + an input vector that the model predicts will achieve it. This is the + basis of latent-variable product and process design (Jaeckle and + MacGregor, 2000). + + Because a PLS model usually retains more components ``A`` than the rank + ``r`` of the response, the target pins down only ``r`` of the ``A`` score + directions and the inversion is underdetermined: a whole + ``(A - r)``-dimensional family of input vectors yields the same + prediction. That family is the *null space*. This method returns the + minimum-norm (direct-inversion) solution together with an orthonormal + basis for the null space, so callers can move along it to satisfy + secondary criteria (cost, safety, operability) without changing the + predicted response. + + For a single response (``r = 1``), García-Carrión et al. (2025) proved + that this null space is the same linear space as the *orthogonal space* + isolated by an O-PLS model with the same total number of components. + + Parameters + ---------- + y_desired : float, array-like, pandas Series/DataFrame, or dict + The desired response, on the original (un-scaled) Y scale. A scalar + is accepted for a single-target model; otherwise supply one value per + target. A Series/DataFrame/dict is aligned to the fitted target + names; a plain array must follow the fitted target order. + null_space_coordinates : np.ndarray, optional + Coordinates along the null-space basis, of length ``A - r`` (the + null-space dimension). When given, the returned solution is + ``tau_direct_inversion + null_space_basis @ null_space_coordinates`` + reconstructed into the input space. All such solutions yield the same + predicted response. When omitted, the minimum-norm (direct-inversion) + solution is returned. + + Returns + ------- + result : sklearn.utils.Bunch + With keys: + + ``x_new`` : pd.Series of shape (n_features,) + The estimated input vector, on the original (un-scaled) X scale. + ``scores`` : pd.Series of length A + The score vector (tau) of the solution. + ``y_hat`` : pd.Series of length n_targets + The model's prediction at ``x_new``; equals ``y_desired`` up to + numerical error, a check that the inversion is consistent. + ``null_space_basis`` : pd.DataFrame of shape (A, A - r) + Orthonormal basis of the null space, in score coordinates. Empty + (zero columns) when ``A == r`` and the solution is unique. + ``null_space_dimension`` : int + ``A - r``, the number of free directions. + ``hotellings_t2`` : float + Hotelling's T² of the solution, to flag extrapolation beyond the + calibration data. Compare against :meth:`hotellings_t2_limit`. + + See Also + -------- + predict : the forward direction, X -> Y. + hotellings_t2_limit : confidence limit to judge ``hotellings_t2``. + + References + ---------- + C. M. Jaeckle and J. F. MacGregor, "Industrial applications of product + design through the inversion of latent variable models", Chemometrics + and Intelligent Laboratory Systems, 50 (2000): 199-210, + DOI: 10.1016/S0169-7439(99)00058-1. + + S. García-Carrión et al., "On the equivalence between null space and + orthogonal space in latent variable regression modeling", Journal of + Chemometrics, 39 (2025): e70057, DOI: 10.1002/cem.70057. + + Examples + -------- + >>> result = pls.invert(y_desired=25.0) + >>> result.x_new # input vector giving the target response + >>> result.null_space_basis # directions that leave the response fixed + >>> pls.predict(result.x_new.to_frame().T) # ~= 25.0 + """ + check_is_fitted(self, "y_loadings_") + + # ``beta_coefficients_`` is only ever a DataFrame (unlike ``y_loadings_``, + # which is an ndarray inside ``_fit_nipals`` before being wrapped), so its + # columns are the reliable source of the fitted target names. + target_names = self.beta_coefficients_.columns + feature_names = self._feature_names + component_names = self._component_names + n_targets = len(target_names) + + # --- Coerce y_desired to a 1-row DataFrame on the original Y scale --- + if isinstance(y_desired, pd.DataFrame): + y_des_df = y_desired.reindex(columns=target_names) + elif isinstance(y_desired, (pd.Series, dict)): + y_des_df = pd.Series(y_desired).reindex(target_names).to_frame().T + else: + y_arr = np.atleast_1d(np.asarray(y_desired, dtype=float)).reshape(1, -1) + if y_arr.shape[1] != n_targets: + raise ValueError( + f"y_desired has {y_arr.shape[1]} value(s); expected {n_targets} to match the fitted " + f"target(s): {list(target_names)}." + ) + y_des_df = pd.DataFrame(y_arr, columns=target_names) + y_des_df.index = [0] + if y_des_df.isna().any(axis=None): + raise ValueError(f"y_desired is missing value(s) for target(s): {list(target_names)}.") + + # Move the target into the scaled Y space the model was fit in. + y_des_scaled = ( + self._y_scaler.transform(y_des_df).to_numpy() if self._y_scaler is not None else y_des_df.to_numpy() + ) # shape (1, M) + + # Score -> scaled-Y map: y_scaled = tau @ q_map, with q_map [A x M]. + q_map = np.asarray(self.y_loadings_).T + # Minimum-norm (direct-inversion) score, tau_DI [1 x A] (paper Eq. 9). + tau = y_des_scaled @ np.linalg.pinv(q_map) + + # Null-space basis: score directions g with g @ q_map = 0 (paper Eq. 14). + # The left singular vectors of q_map for its zero singular values span it. + u_mat, sing, _ = np.linalg.svd(q_map) + tol = max(q_map.shape) * np.finfo(float).eps * (sing[0] if sing.size else 0.0) + rank = int((sing > tol).sum()) + null_space_basis = u_mat[:, rank:] # [A x (A - rank)] + null_space_dimension = null_space_basis.shape[1] + + # Optionally move along the null space; every such point predicts y_desired. + if null_space_coordinates is not None: + coords = np.atleast_1d(np.asarray(null_space_coordinates, dtype=float)) + if coords.shape[0] != null_space_dimension: + raise ValueError( + f"null_space_coordinates has {coords.shape[0]} value(s); expected " + f"{null_space_dimension} to match the null-space dimension (A - r)." + ) + tau = tau + (null_space_basis @ coords).reshape(1, -1) + + # Reconstruct the input vector from its scores, then un-scale to raw units. + x_new_scaled = tau @ self._x_loadings.T # [1 x K] + x_new_df = pd.DataFrame(x_new_scaled, columns=feature_names, index=[0]) + if self._x_scaler is not None: + x_new_df = self._x_scaler.inverse_transform(x_new_df) + + # Prediction at the solution (should equal y_desired up to rounding). + y_hat_scaled = tau @ q_map # [1 x M] + y_hat_df = pd.DataFrame(y_hat_scaled, columns=target_names, index=[0]) + if self._y_scaler is not None: + y_hat_df = self._y_scaler.inverse_transform(y_hat_df) + + t2 = float(np.sum((tau.ravel() / self.scaling_factor_for_scores_.to_numpy()) ** 2)) + ns_columns = [f"NS{i}" for i in range(1, null_space_dimension + 1)] + + return Bunch( + x_new=pd.Series(x_new_df.to_numpy().ravel(), index=feature_names, name="x_new"), + scores=pd.Series(tau.ravel(), index=component_names, name="scores"), + y_hat=pd.Series(y_hat_df.to_numpy().ravel(), index=target_names, name="y_hat"), + null_space_basis=pd.DataFrame(null_space_basis, index=component_names, columns=ns_columns), + null_space_dimension=null_space_dimension, + hotellings_t2=t2, + ) + @classmethod def select_n_components( # noqa: C901, PLR0912, PLR0913, PLR0915 cls, From 117c352ee9aebfe21995e13a9feab1f0e00e19fa Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 19:17:13 +0000 Subject: [PATCH 2/7] Ignore ruff CPY001 (missing-copyright-notice) preview rule A newer ruff release promoted CPY001 from preview to stable. With select=["ALL"], it then fired on ~every file (292 hits) because the repo's "# (c) Kevin Dunn" house notice does not match ruff's default copyright regex, turning the lint gate red on all PRs regardless of their diff. Copyright notices are managed in the repo's own style, not via ruff, so ignore CPY001 to make the gate independent of the installed ruff version. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DsTypY4FcRbrXRxorb9raz --- pyproject.toml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 7fc5fc3..f8b426c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -208,6 +208,13 @@ exclude = [ # Enable ALL Pyflakes codes by default. select = ["ALL"] ignore = [ + # Copyright notices are managed in the repo's own "# (c) Kevin Dunn" house + # style, not via ruff. CPY001 is a preview rule with a default regex that + # requires the literal word "Copyright"; once a ruff release promoted it to + # stable it fired on nearly every file. Keep it off so the lint gate does not + # depend on the installed ruff version. + "CPY001", # Missing copyright notice at top of file. + "C408", # Unnecessary call - rewrite as a literal. Disagree with this. # Arguments From e2afd6863d993fadde8bfd8591c53ff1d2c41871 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 19:20:50 +0000 Subject: [PATCH 3/7] Test PLS.invert; pin ruff <0.16 to stabilise the lint gate Adds tests/test_multivariate_inversion.py: reproduces the paper's Case Study 1 (Table 3), checks the invert round-trip, null-space invariance, input coercion, T2 reporting, the unique single-component case, the cheddar cheese hold-out design scenario, and the multi-response null-space dimension on the real LDPE data. Vendors cheddar-cheese.csv as an offline fixture. Also upper-bounds ruff (>=0.11.0,<0.16) in both dev dependency groups. CI's 'uv sync' has no --frozen flag and the old 'ruff>=0.11.0' let it resolve to a newer ruff than the lock, which promoted CPY001 and changed RUF100/ISC handling and turned the lint gate red on every PR. This mirrors the existing mypy upper bound and keeps the gate independent of the latest ruff release. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DsTypY4FcRbrXRxorb9raz --- pyproject.toml | 4 +- .../datasets/multivariate/cheddar-cheese.csv | 31 ++++ tests/test_multivariate_inversion.py | 174 ++++++++++++++++++ 3 files changed, 207 insertions(+), 2 deletions(-) create mode 100644 src/process_improve/datasets/multivariate/cheddar-cheese.csv create mode 100644 tests/test_multivariate_inversion.py diff --git a/pyproject.toml b/pyproject.toml index f8b426c..4658d89 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -116,7 +116,7 @@ dev = [ "pytest-cov>=7.0.0", "pytest-xdist>=3.8.0", "pydata-sphinx-theme>=0.16.1", - "ruff>=0.11.0", + "ruff>=0.11.0,<0.16", # upper-bound so a future ruff release cannot silently break the lint gate (see mypy) "sphinx>=8.1.3", "nbsphinx>=0.9.5", "ipykernel>=6.29", @@ -158,7 +158,7 @@ dev = [ "pytest-cov>=7.0.0", "pytest-xdist>=3.8.0", "pydata-sphinx-theme>=0.16.1", - "ruff>=0.11.0", + "ruff>=0.11.0,<0.16", # upper-bound so a future ruff release cannot silently break the lint gate (see mypy) "sphinx>=8.1.3", "nbsphinx>=0.9.5", "ipykernel>=6.29", diff --git a/src/process_improve/datasets/multivariate/cheddar-cheese.csv b/src/process_improve/datasets/multivariate/cheddar-cheese.csv new file mode 100644 index 0000000..bc40a73 --- /dev/null +++ b/src/process_improve/datasets/multivariate/cheddar-cheese.csv @@ -0,0 +1,31 @@ +Case,Acetic,H2S,Lactic,Taste +1,4.543,3.135,0.86,12.3 +2,5.159,5.043,1.53,20.9 +3,5.366,5.438,1.57,39 +4,5.759,7.496,1.81,47.9 +5,4.663,3.807,0.99,5.6 +6,5.697,7.601,1.09,25.9 +7,5.892,8.726,1.29,37.3 +8,6.078,7.966,1.78,21.9 +9,4.898,3.85,1.29,18.1 +10,5.242,4.174,1.58,21 +11,5.74,6.142,1.68,34.9 +12,6.446,7.908,1.9,57.2 +13,4.477,2.996,1.06,0.7 +14,5.236,4.942,1.3,25.9 +15,6.151,6.752,1.52,54.9 +16,6.365,9.588,1.74,40.9 +17,4.787,3.912,1.16,15.9 +18,5.412,4.7,1.49,6.4 +19,5.247,6.174,1.63,18 +20,5.438,9.064,1.99,38.9 +21,4.564,4.949,1.15,14 +22,5.298,5.22,1.33,15.2 +23,5.455,9.242,1.44,32 +24,5.855,10.199,2.0,56.71 +25,5.366,3.664,1.31,16.8 +26,6.043,3.219,1.46,11.6 +27,6.458,6.962,1.72,26.5 +28,5.328,3.912,1.25,0.7 +29,5.802,6.685,1.08,13.4 +30,6.176,4.787,1.25,5.5 diff --git a/tests/test_multivariate_inversion.py b/tests/test_multivariate_inversion.py new file mode 100644 index 0000000..9d757a2 --- /dev/null +++ b/tests/test_multivariate_inversion.py @@ -0,0 +1,174 @@ +# (c) Kevin Dunn, 2010-2026. MIT License. +"""Tests for PLS model inversion and the null space (``PLS.invert``). + +The reference numbers come from García-Carrión et al., "On the equivalence +between null space and orthogonal space in latent variable regression modeling", +Journal of Chemometrics, 39 (2025): e70057, DOI: 10.1002/cem.70057. Case Study 1 +of that paper (a 6x5 simulated dataset) has published direct-inversion and +null-space values that we reproduce here. +""" + +from __future__ import annotations + +import pathlib + +import numpy as np +import pandas as pd +import pytest + +from process_improve.multivariate.methods import PLS, MCUVScaler + + +@pytest.fixture +def case_study_1() -> tuple[pd.DataFrame, pd.DataFrame]: + """Paper Case Study 1: X [6x5] and single response y [6x1] (Table 1).""" + X = pd.DataFrame( + [ + [5.43, 7.54, 125.64, 58.51, 50.49], + [5.43, 15.97, 126.20, 258.48, 74.44], + [99.23, 7.54, 9893.38, 59.29, 737.15], + [99.23, 15.97, 9765.16, 254.11, 1576.28], + [52.33, 11.76, 2787.64, 139.21, 538.76], + [52.33, 11.76, 2849.95, 135.67, 630.73], + ], + columns=["x1", "x2", "x3", "x4", "x5"], + ) + y = pd.DataFrame({"y": [61.85, 278.99, 307.89, 436.4, 266.08, 260.52]}) + return X, y + + +def test_invert_reproduces_paper_case_study_1(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """Direct-inversion score and null-space basis match the paper's Table 3.""" + X, y = case_study_1 + model = PLS(n_components=2).fit(X, y) + result = model.invert(168.23) + + # Table 3: tau_DI = (-1.4189, -0.2231); the published X are rounded, so a + # loose tolerance is appropriate. + assert result.scores.to_numpy() == pytest.approx([-1.4189, -0.2231], abs=5e-3) + + # Table 3: null-space basis G_NS = (-0.1553, 0.9879). A basis vector is only + # defined up to sign, so compare the absolute direction. + g_ns = result.null_space_basis.to_numpy().ravel() + assert np.abs(g_ns) == pytest.approx([0.1553, 0.9879], abs=5e-3) + assert result.null_space_dimension == 1 + + # The prediction at the solution equals the requested target. + assert float(result.y_hat.iloc[0]) == pytest.approx(168.23, abs=1e-6) + + +def test_invert_round_trip(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """predict(invert(y_des).x_new) recovers y_des.""" + X, y = case_study_1 + model = PLS(n_components=2).fit(X, y) + for target in (100.0, 168.23, 400.0): + result = model.invert(target) + y_back = float(model.predict(result.x_new.to_frame().T).iloc[0, 0]) + assert y_back == pytest.approx(target, abs=1e-6) + + +def test_null_space_moves_leave_prediction_unchanged(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """Moving along the null space changes the inputs but not the prediction.""" + X, y = case_study_1 + model = PLS(n_components=2).fit(X, y) + base = model.invert(168.23) + for coord in (-2.0, 0.5, 3.7): + moved = model.invert(168.23, null_space_coordinates=[coord]) + # Inputs differ... + assert not np.allclose(moved.x_new.to_numpy(), base.x_new.to_numpy()) + # ...but the predicted response is identical. + y_moved = float(model.predict(moved.x_new.to_frame().T).iloc[0, 0]) + assert y_moved == pytest.approx(168.23, abs=1e-6) + + +def test_invert_accepts_scalar_series_and_array(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """A scalar, a 1-element array, and a named Series give the same result.""" + X, y = case_study_1 + model = PLS(n_components=2).fit(X, y) + from_scalar = model.invert(168.23).x_new + from_array = model.invert(np.array([168.23])).x_new + from_series = model.invert(pd.Series({"y": 168.23})).x_new + assert from_array.to_numpy() == pytest.approx(from_scalar.to_numpy()) + assert from_series.to_numpy() == pytest.approx(from_scalar.to_numpy()) + + +def test_invert_reports_hotellings_t2(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """The solution's T2 is the sum of squared standardized scores and is finite.""" + X, y = case_study_1 + model = PLS(n_components=2).fit(X, y) + result = model.invert(168.23) + expected_t2 = float( + np.sum((result.scores.to_numpy() / model.scaling_factor_for_scores_.to_numpy()) ** 2) + ) + assert result.hotellings_t2 == pytest.approx(expected_t2) + assert result.hotellings_t2 >= 0.0 + + +def test_invert_wrong_coordinate_length_raises(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """A mismatched null_space_coordinates length is a clear error.""" + X, y = case_study_1 + model = PLS(n_components=2).fit(X, y) # null-space dimension is 1 + with pytest.raises(ValueError, match="null-space dimension"): + model.invert(168.23, null_space_coordinates=[1.0, 2.0]) + + +def test_invert_single_component_has_empty_null_space() -> None: + """With one component the inversion is unique: a zero-width null space.""" + rng = np.random.default_rng(0) + X = pd.DataFrame(rng.standard_normal((30, 4)), columns=[f"x{i}" for i in range(4)]) + beta = rng.standard_normal((4, 1)) + y = pd.DataFrame(X.to_numpy() @ beta, columns=["y"]) + model = PLS(n_components=1).fit(X, y) + result = model.invert(1.0) + assert result.null_space_dimension == 0 + assert result.null_space_basis.shape[1] == 0 + + +def test_invert_cheddar_cheese_holdout() -> None: + """Design toward the taste of held-out cheeses; round-trip must hold. + + Mirrors the pid-book worked example: train on cheeses 5-30, then invert the + two-component model toward each of the four held-out cheeses' taste scores. + """ + folder = pathlib.Path(__file__).parents[1] / "src" / "process_improve" / "datasets" / "multivariate" + path = folder / "cheddar-cheese.csv" + if not path.exists(): + pytest.skip("cheddar-cheese.csv fixture not present") + cheese = pd.read_csv(path) + x_cols = ["Acetic", "H2S", "Lactic"] + train = cheese.iloc[4:] + hold = cheese.iloc[:4] + model = PLS(n_components=2).fit(train[x_cols], train[["Taste"]]) + for i in range(len(hold)): + target = float(hold["Taste"].iloc[i]) + result = model.invert(target) + assert result.null_space_dimension == 1 + y_back = float(model.predict(result.x_new.to_frame().T).iloc[0, 0]) + assert y_back == pytest.approx(target, abs=1e-6) + + +def test_invert_multi_response_null_space_dimension() -> None: + """On the real LDPE data the null-space dimension is A - rank(Y). + + LDPE has five (collinear) quality variables. Inverting toward a full + five-target vector gives a null space of dimension A - rank(Y), and the + reconstructed inputs reproduce the requested targets. + """ + folder = pathlib.Path(__file__).parents[1] / "src" / "process_improve" / "datasets" / "multivariate" + values = pd.read_csv(folder / "LDPE" / "LDPE.csv", index_col=0) + X = values.iloc[:, :14] + Y = values.iloc[:, 14:] + Xs = MCUVScaler().fit_transform(X) + Ys = MCUVScaler().fit_transform(Y) + n_components = 6 + model = PLS(n_components=n_components).fit(Xs, Ys) + + # A desired target taken as one of the training rows' responses. + y_desired = Ys.iloc[0] + result = model.invert(y_desired) + + rank_y = np.linalg.matrix_rank(Ys.to_numpy()) + assert result.null_space_dimension == n_components - rank_y + + y_back = model.predict(result.x_new.to_frame().T).to_numpy().ravel() + assert y_back == pytest.approx(y_desired.to_numpy(), abs=1e-6) From 81e4f316118382e1018584ec981217b91931d9c8 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 19:29:28 +0000 Subject: [PATCH 4/7] Cover remaining PLS.invert branches (DataFrame/dict input, raises, scale=False) Adds tests for the DataFrame and dict forms of y_desired, the wrong target-count and NaN-target error paths, and the unscaled (scale=False) code path, bringing patch coverage of invert() to full. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DsTypY4FcRbrXRxorb9raz --- tests/test_multivariate_inversion.py | 48 ++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/tests/test_multivariate_inversion.py b/tests/test_multivariate_inversion.py index 9d757a2..4e9c454 100644 --- a/tests/test_multivariate_inversion.py +++ b/tests/test_multivariate_inversion.py @@ -92,6 +92,54 @@ def test_invert_accepts_scalar_series_and_array(case_study_1: tuple[pd.DataFrame assert from_series.to_numpy() == pytest.approx(from_scalar.to_numpy()) +def test_invert_accepts_dataframe_and_dict(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """A one-row DataFrame and a dict give the same result as a scalar.""" + X, y = case_study_1 + model = PLS(n_components=2).fit(X, y) + from_scalar = model.invert(168.23).x_new + from_frame = model.invert(pd.DataFrame({"y": [168.23]})).x_new + from_dict = model.invert({"y": 168.23}).x_new + assert from_frame.to_numpy() == pytest.approx(from_scalar.to_numpy()) + assert from_dict.to_numpy() == pytest.approx(from_scalar.to_numpy()) + + +def test_invert_wrong_target_count_raises(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """Passing the wrong number of target values is a clear error.""" + X, y = case_study_1 + model = PLS(n_components=2).fit(X, y) # single target + with pytest.raises(ValueError, match="expected 1"): + model.invert(np.array([1.0, 2.0])) + + +def test_invert_missing_target_value_raises(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """A NaN in the desired response is rejected.""" + X, y = case_study_1 + model = PLS(n_components=2).fit(X, y) + with pytest.raises(ValueError, match="missing value"): + model.invert(np.nan) + + +def test_invert_on_prescaled_model_scale_false(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """Inversion works when the model was fit with scale=False on pre-scaled data. + + Exercises the un-scaled code path (no internal X/Y scaler). The desired + response must be supplied on the same (already-scaled) Y scale the model + saw, and the round-trip must still hold. + """ + X, y = case_study_1 + x_scaler = MCUVScaler().fit(X) + y_scaler = MCUVScaler().fit(y) + Xs = x_scaler.transform(X) + ys = y_scaler.transform(y) + model = PLS(n_components=2, scale=False).fit(Xs, ys) + + target_scaled = float(y_scaler.transform(pd.DataFrame({"y": [168.23]})).iloc[0, 0]) + result = model.invert(target_scaled) + y_back = float(model.predict(result.x_new.to_frame().T).iloc[0, 0]) + assert y_back == pytest.approx(target_scaled, abs=1e-6) + assert result.null_space_dimension == 1 + + def test_invert_reports_hotellings_t2(case_study_1: tuple[pd.DataFrame, pd.DataFrame]) -> None: """The solution's T2 is the sum of squared standardized scores and is finite.""" X, y = case_study_1 From 9a3b3d9b9cdc5a51c54d494dd7dc9f23e91ef87b Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 19:44:08 +0000 Subject: [PATCH 5/7] Add OPLS estimator (Trygg-Wold O-PLS) with single-division inversion New OPLS class for single-response Orthogonal PLS: fit (Trygg-Wold NIPALS), transform (predictive scores), correct (orthogonal-signal-corrected X), predict, and invert. OPLS.invert() reaches a desired response through one division and returns the orthogonal-space basis; for a single response this is the same set of designs as PLS.invert()'s null space (Garcia-Carrion et al., 2025). Tests confirm O-PLS(1; A-1) reproduces A-component PLS predictions and beta to ~1e-14, the orthogonal scores are Y-orthogonal, and the orthogonal space coincides with the PLS null space (cosine 1.0). Exported from methods.py and the multivariate package; added to the sklearn compatibility checks. Folded into the 1.59.0 changelog entry alongside PLS.invert (same PR / release). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DsTypY4FcRbrXRxorb9raz --- CHANGELOG.md | 9 + src/process_improve/multivariate/__init__.py | 2 + src/process_improve/multivariate/_opls.py | 458 +++++++++++++++++++ src/process_improve/multivariate/methods.py | 2 + tests/test_multivariate_opls.py | 216 +++++++++ tests/test_sklearn_compat.py | 14 +- 6 files changed, 699 insertions(+), 2 deletions(-) create mode 100644 src/process_improve/multivariate/_opls.py create mode 100644 tests/test_multivariate_opls.py diff --git a/CHANGELOG.md b/CHANGELOG.md index a4a36d4..5c243b6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,15 @@ those changes. ### Added +- `OPLS`: an Orthogonal Projections to Latent Structures estimator for a single + response (Trygg-Wold NIPALS), extracting one Y-predictive component and + `n_orthogonal_components` Y-orthogonal components. It exposes `fit`, + `transform` (predictive scores), `correct` (the orthogonal-signal-corrected + X), `predict`, and `invert`. `OPLS.invert()` finds the inputs for a desired + response through a single division and returns the orthogonal-space basis; + for a single response this is the same set of designs as `PLS.invert()`'s null + space (García-Carrión et al., 2025). An A-component PLS model and an + O-PLS(1; A - 1) model give identical predictions and regression coefficients. - `PLS.invert()`: PLS model inversion for latent-variable product and process design. Given a desired response on the original Y scale, it returns the minimum-norm (direct-inversion) input vector that the model predicts will diff --git a/src/process_improve/multivariate/__init__.py b/src/process_improve/multivariate/__init__.py index 33a2a83..0ba9fbd 100644 --- a/src/process_improve/multivariate/__init__.py +++ b/src/process_improve/multivariate/__init__.py @@ -1,6 +1,7 @@ """Multivariate analysis: PCA, PLS, TPLS, scaling, and diagnostic plots.""" from process_improve.multivariate.methods import ( + OPLS, PCA, PLS, TPLS, @@ -31,6 +32,7 @@ ) __all__ = [ + "OPLS", "PCA", "PLS", "TPLS", diff --git a/src/process_improve/multivariate/_opls.py b/src/process_improve/multivariate/_opls.py new file mode 100644 index 0000000..d4b53d0 --- /dev/null +++ b/src/process_improve/multivariate/_opls.py @@ -0,0 +1,458 @@ +# (c) Kevin Dunn, 2010-2026. MIT License. +"""Orthogonal Projections to Latent Structures (O-PLS) for a single response. + +O-PLS (Trygg and Wold, 2002) is a variant of PLS that splits the systematic +variation in ``X`` into two parts: a single *predictive* component that is +correlated with the response ``y``, and one or more *Y-orthogonal* components +that carry systematic variation in ``X`` unrelated to ``y``. Filtering out the +orthogonal part leaves a model with the same predictive ability as an ordinary +PLS model fitted with the same total number of components, but with the +response-relevant variation concentrated in one component, which aids +interpretation. + +This module implements the single-response case, O-PLS(1; ``n_orthogonal``), +following the notation of García-Carrión et al. (2025). That paper proves that +the *orthogonal space* isolated here is the same linear space as the *null +space* of an inverted PLS model with the same total number of components, so +:meth:`OPLS.invert` and :meth:`PLS.invert` describe the same set of designs; +O-PLS reaches it through a single division. + +References +---------- +J. Trygg and S. Wold, "Orthogonal projections to latent structures (O-PLS)", +Journal of Chemometrics, 16 (2002): 119-128, DOI: 10.1002/cem.695. + +S. García-Carrión et al., "On the equivalence between null space and orthogonal +space in latent variable regression modeling", Journal of Chemometrics, 39 +(2025): e70057, DOI: 10.1002/cem.70057. +""" + +from __future__ import annotations + +import typing + +import numpy as np +import pandas as pd +from sklearn.base import BaseEstimator, RegressorMixin, TransformerMixin +from sklearn.utils import Bunch +from sklearn.utils.validation import check_is_fitted, validate_data + +from .._linalg import safe_inverse +from ._base import _LatentVariableModel, _LazyFrame +from ._common import epsqrt +from ._preprocessing import MCUVScaler + +if typing.TYPE_CHECKING: + from ._pls import DataMatrix + + +class OPLS(_LatentVariableModel, RegressorMixin, TransformerMixin, BaseEstimator): + """Orthogonal PLS for a single response, O-PLS(1; ``n_orthogonal_components``). + + Extracts one Y-predictive component and ``n_orthogonal_components`` + Y-orthogonal components from ``X`` via the Trygg-Wold NIPALS algorithm. The + predictive component carries the variation used to predict ``y``; the + orthogonal components carry systematic ``X`` variation unrelated to ``y``. + + Parameters + ---------- + n_orthogonal_components : int + Number of Y-orthogonal components to extract (``Ao``). The model has a + total of ``1 + n_orthogonal_components`` components, matching an ordinary + PLS model with that many latent variables. + scale : bool, default=True + Mean-center and unit-variance-scale X and y internally before fitting + (with :class:`MCUVScaler`). Predictions and ``beta_coefficients_`` are + returned on the original data scale. Set ``False`` when the inputs are + already scaled. + max_iter : int, default=1000 + Reserved for API parity with :class:`PLS`; the closed-form single- + response algorithm here does not iterate. + tol : float, default=sqrt(machine epsilon) + Numerical tolerance used when guarding rank-deficient projections. + copy : bool, default=True + Whether to copy X and y before fitting. + + Attributes (after fitting) + -------------------------- + predictive_scores_ : pd.DataFrame of shape (n_samples, 1) + Predictive score ``t_p``. + predictive_weights_ : pd.Series of length n_features + Predictive weight ``w_p``. + predictive_loadings_ : pd.Series of length n_features + Predictive loading ``p_p``. + orthogonal_scores_ : pd.DataFrame of shape (n_samples, n_orthogonal) + Y-orthogonal scores ``T_o``. + orthogonal_weights_ : pd.DataFrame of shape (n_features, n_orthogonal) + Y-orthogonal weights ``W_o``. + orthogonal_loadings_ : pd.DataFrame of shape (n_features, n_orthogonal) + Y-orthogonal loadings ``P_o``. + y_loadings_ : float + Predictive y-loading ``q_p`` (a scalar for a single response). + beta_coefficients_ : pd.DataFrame of shape (n_features, 1) + Regression coefficients linking X directly to y, on the original scale. + scores_ : pd.DataFrame of shape (n_samples, 1 + n_orthogonal) + Predictive score followed by the orthogonal scores, labelled + ``["t_predictive", "t_orthogonal_1", ...]``. Lets the inherited + :meth:`score_plot` draw the predictive-vs-orthogonal score plot. + x_filtered_ : pd.DataFrame of shape (n_samples, n_features) + The orthogonal-signal-corrected X (``X`` with the orthogonal variation + removed), on the scaled fitting scale. + spe_ : pd.DataFrame + Per-row SPE after reconstructing X from all components. + hotellings_t2_ : pd.DataFrame + Cumulative Hotelling's T2 over the combined score space. + + See Also + -------- + PLS : ordinary projection to latent structures; ``PLS.invert`` gives the + same designs as :meth:`OPLS.invert` for a single response. + + References + ---------- + J. Trygg and S. Wold, "Orthogonal projections to latent structures (O-PLS)", + Journal of Chemometrics, 16 (2002): 119-128, DOI: 10.1002/cem.695. + + Examples + -------- + >>> import pandas as pd + >>> from process_improve.multivariate.methods import OPLS + >>> X = pd.DataFrame({"a": [1.0, 2, 3, 4], "b": [1.0, 0, 1, 0]}) + >>> y = pd.DataFrame({"y": [1.0, 2, 3, 4]}) + >>> model = OPLS(n_orthogonal_components=1).fit(X, y) + >>> model.predict(X).shape + (4, 1) + """ + + _ATTRIBUTE_RENAMES: typing.ClassVar[dict[str, str]] = {} + _RENAME_CONTEXT: typing.ClassVar[str] = "OPLS" + + # Lazily-built DataFrame views over the private ndarrays (ENG-18), matching + # the PLS convention so the inherited score/SPE plots work unchanged. + scores_ = _LazyFrame("_scores", index="_sample_index", columns="_component_names") + spe_ = _LazyFrame("_spe", index="_sample_index", columns="_component_names") + + def __init__( + self, + n_orthogonal_components: int, + *, + scale: bool = True, + max_iter: int = 1000, + tol: float = epsqrt, + copy: bool = True, + ): + self.n_orthogonal_components = n_orthogonal_components + self.scale = scale + self.max_iter = max_iter + self.tol = tol + self.copy = copy + + def __sklearn_tags__(self): + """Declare sklearn capability tags (sklearn 1.6+).""" + tags = super().__sklearn_tags__() + tags.target_tags.required = True + tags.target_tags.single_output = True + return tags + + def get_feature_names_out(self, input_features=None) -> np.ndarray: # noqa: ANN001, ARG002 + """Return the component labels produced by :meth:`transform`.""" + check_is_fitted(self, "predictive_weights_") + return np.asarray(self._component_names, dtype=object) + + def fit(self, X: DataMatrix, Y: DataMatrix) -> OPLS: # noqa: C901, PLR0915 + """Fit the O-PLS model. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) + Predictor block. + Y : array-like of shape (n_samples,) or (n_samples, 1) + Single response. A 2-D column or a 1-D vector are both accepted. + + Returns + ------- + OPLS + The fitted model (``self``). + """ + if hasattr(Y, "ndim") and Y.ndim == 1: + Y = Y.to_frame() if isinstance(Y, pd.Series) else pd.DataFrame(np.asarray(Y).reshape(-1, 1)) + + sample_index: pd.Index | None = X.index if isinstance(X, pd.DataFrame) else None + feature_columns: pd.Index | None = X.columns if isinstance(X, pd.DataFrame) else None + X_arr = validate_data( + self, + X, + reset=True, + accept_sparse=False, + ensure_min_samples=2, + ensure_min_features=1, + dtype="numeric", + ensure_all_finite=True, + ) + if feature_columns is None: + feature_columns = pd.RangeIndex(X_arr.shape[1]) + if sample_index is None: + sample_index = pd.RangeIndex(X_arr.shape[0]) + X = pd.DataFrame(X_arr, index=sample_index, columns=feature_columns) + if not isinstance(Y, pd.DataFrame): + Y = pd.DataFrame(Y, index=sample_index) + if Y.shape[1] != 1: + raise ValueError(f"OPLS supports a single response only; Y has {Y.shape[1]} columns.") + + self.n_samples_ = X.shape[0] + self.n_targets_ = 1 + n_ortho = int(self.n_orthogonal_components) + if n_ortho < 0: + raise ValueError("n_orthogonal_components must be >= 0.") + self.n_components = 1 + n_ortho + + self._x_scaler: MCUVScaler | None = None + self._y_scaler: MCUVScaler | None = None + if self.scale: + self._x_scaler = MCUVScaler().fit(X) + self._y_scaler = MCUVScaler().fit(Y) + X = self._x_scaler.transform(X) + Y = self._y_scaler.transform(Y) + + Xmat = X.to_numpy(dtype=float) + y = Y.to_numpy(dtype=float) # (N, 1) + n_features = Xmat.shape[1] + target_name = Y.columns + + # --- Trygg-Wold O-PLS NIPALS (single response) --- + # Predictive weight direction (unit norm), fixed throughout. + w_p = Xmat.T @ y + w_p = w_p / np.linalg.norm(w_p) # (K, 1) + + W_o = np.zeros((n_features, n_ortho)) + P_o = np.zeros((n_features, n_ortho)) + T_o = np.zeros((self.n_samples_, n_ortho)) + Xd = Xmat.copy() + for i in range(n_ortho): + t = Xd @ w_p # (N, 1) + p = (Xd.T @ t) / (t.T @ t).item() # (K, 1) + # Orthogonal weight: the part of the loading orthogonal to w_p. + w_o = p - (w_p.T @ p).item() * w_p + w_o = w_o / np.linalg.norm(w_o) + t_o = Xd @ w_o + p_o = (Xd.T @ t_o) / (t_o.T @ t_o).item() + Xd = Xd - t_o @ p_o.T # remove this orthogonal component + W_o[:, [i]] = w_o + P_o[:, [i]] = p_o + T_o[:, [i]] = t_o + + # Predictive component on the filtered (orthogonal-corrected) X. + t_p = Xd @ w_p # (N, 1) + p_p = (Xd.T @ t_p) / (t_p.T @ t_p).item() # (K, 1) + q_p = (y.T @ t_p).item() / (t_p.T @ t_p).item() # scalar + + # OSC filter operator so predict/transform reproduce the fit exactly. + if n_ortho > 0: + filt = np.eye(n_features) - W_o @ safe_inverse(P_o.T @ W_o, what="(P_o' @ W_o)") @ P_o.T + else: + filt = np.eye(n_features) + self._osc_filter = filt # (K, K) + self._w_p = w_p + self._W_o = W_o + self._P_o = P_o + self._p_p = p_p + self._q_p = q_p + + # beta on the (optionally) scaled space, then map back to original units. + beta_scaled = (filt @ w_p) * q_p # (K, 1) + beta = beta_scaled.copy() + if self._x_scaler is not None and self._y_scaler is not None: + x_scale = self._x_scaler.scale_.to_numpy()[:, np.newaxis] + y_scale = self._y_scaler.scale_.to_numpy()[np.newaxis, :] + beta = beta * (y_scale / x_scale) + + # --- Public / backing attributes --- + component_names = ["t_predictive"] + [f"t_orthogonal_{i}" for i in range(1, n_ortho + 1)] + self._sample_index = X.index + self._feature_names = X.columns + self._component_names = component_names + + self.predictive_scores_ = pd.DataFrame(t_p, index=X.index, columns=["t_predictive"]) + self.predictive_weights_ = pd.Series(w_p.ravel(), index=X.columns, name="w_predictive") + self.predictive_loadings_ = pd.Series(p_p.ravel(), index=X.columns, name="p_predictive") + ortho_names = [f"o{i}" for i in range(1, n_ortho + 1)] + self.orthogonal_scores_ = pd.DataFrame(T_o, index=X.index, columns=ortho_names) + self.orthogonal_weights_ = pd.DataFrame(W_o, index=X.columns, columns=ortho_names) + self.orthogonal_loadings_ = pd.DataFrame(P_o, index=X.columns, columns=ortho_names) + self.y_loadings_ = q_p + self.beta_coefficients_ = pd.DataFrame(beta, index=X.columns, columns=target_name) + + # Combined score matrix [t_p | T_o] backs scores_ and the geometry plots. + self._scores = np.hstack([t_p, T_o]) + self.explained_variance_ = np.diag(self._scores.T @ self._scores) / max(1, self.n_samples_ - 1) + self.scaling_factor_for_scores_ = pd.Series( + np.sqrt(self.explained_variance_), + index=component_names, + name="Standard deviation per score", + ) + + # Filtered X and full reconstruction for SPE. + self.x_filtered_ = pd.DataFrame(Xmat @ filt, index=X.index, columns=X.columns) + x_hat = t_p @ p_p.T + T_o @ P_o.T + residuals = Xmat - x_hat + self._spe = np.sqrt(np.sum(residuals**2, axis=1)).reshape(-1, 1) * np.ones((1, self.n_components)) + + # Cumulative Hotelling's T2 across the combined score space. + std = self.scaling_factor_for_scores_.to_numpy() + t2_cum = np.cumsum((self._scores / std) ** 2, axis=1) + self.hotellings_t2_ = pd.DataFrame(t2_cum, index=X.index, columns=component_names) + + return self + + def transform(self, X: DataMatrix, Y: DataMatrix | None = None) -> pd.DataFrame: # noqa: ARG002 + """Project new data and return the predictive score ``t_p``. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) + Data to project. + Y : ignored + Present for sklearn pipeline API compatibility. + + Returns + ------- + pd.DataFrame of shape (n_samples, 1) + The predictive score for each row. + """ + X_scaled, index = self._prepare(X) + t_p = (X_scaled @ self._osc_filter) @ self._w_p + return pd.DataFrame(t_p, index=index, columns=["t_predictive"]) + + def correct(self, X: DataMatrix) -> pd.DataFrame: + """Return the orthogonal-signal-corrected X (orthogonal variation removed). + + This is O-PLS's original use as a preprocessing filter: the returned + matrix has the Y-orthogonal variation stripped out, on the scaled + fitting scale. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) + + Returns + ------- + pd.DataFrame of shape (n_samples, n_features) + The filtered X. + """ + X_scaled, index = self._prepare(X) + return pd.DataFrame(X_scaled @ self._osc_filter, index=index, columns=self._feature_names) + + def predict(self, X: DataMatrix) -> pd.DataFrame: + """Predict the response for new observations. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) + + Returns + ------- + pd.DataFrame of shape (n_samples, 1) + Predicted response on the original (un-scaled) y scale. + """ + t_p = self.transform(X) + y_hat = t_p.to_numpy() * self._q_p + y_hat_df = pd.DataFrame(y_hat, index=t_p.index, columns=self.beta_coefficients_.columns) + if self._y_scaler is not None: + y_hat_df = self._y_scaler.inverse_transform(y_hat_df) + return y_hat_df + + def invert(self, y_desired: float | np.ndarray | pd.Series | pd.DataFrame) -> Bunch: + r"""Invert the O-PLS model: find inputs that yield a desired response. + + For a single response the inversion reduces to a single division. The + predictive score that achieves the target is + ``tau_p = y_desired / q_p``, and the *orthogonal space* (spanned by the + orthogonal loadings ``P_o``) can be varied freely without changing the + prediction. This is the same set of designs that :meth:`PLS.invert` + returns as its null space (García-Carrión et al., 2025), reached here by + one division instead of solving an underdetermined system. + + Parameters + ---------- + y_desired : float, array-like, pandas Series/DataFrame + Desired response on the original (un-scaled) y scale. + + Returns + ------- + result : sklearn.utils.Bunch + With keys: + + ``x_new`` : pd.Series of shape (n_features,) + Input vector (minimum-norm; zero orthogonal contribution) on the + original X scale. + ``predictive_score`` : float + The predictive score ``tau_p`` achieving the target. + ``y_hat`` : float + Prediction at ``x_new`` (equals ``y_desired`` up to rounding). + ``orthogonal_space_basis`` : pd.DataFrame of shape (n_features, n_orthogonal) + Basis of the orthogonal space in the input space (the orthogonal + loadings). Adding any combination of these columns to ``x_new`` + leaves the prediction unchanged. + ``orthogonal_space_dimension`` : int + ``n_orthogonal_components``. + + See Also + -------- + PLS.invert : the equivalent inversion via the null space. + """ + check_is_fitted(self, "predictive_weights_") + + y_arr = np.atleast_1d(np.asarray(_as_scalar(y_desired), dtype=float)) + if y_arr.size != 1: + raise ValueError("OPLS.invert expects a single desired response value.") + y_value = float(y_arr[0]) + if not np.isfinite(y_value): + raise ValueError("y_desired must be a finite value.") + + # Scale the target to the fitting space. + if self._y_scaler is not None: + y_scaled = (y_value - float(self._y_scaler.center_.iloc[0])) / float(self._y_scaler.scale_.iloc[0]) + else: + y_scaled = y_value + + tau_p = y_scaled / self._q_p # single division (paper Eq. 28) + + # Minimum-norm input: predictive contribution only, no orthogonal part. + x_new_scaled = (tau_p * self._p_p).reshape(1, -1) # (1, K) + x_new_df = pd.DataFrame(x_new_scaled, columns=self._feature_names, index=[0]) + if self._x_scaler is not None: + x_new_df = self._x_scaler.inverse_transform(x_new_df) + + y_hat_value = float(self.predict(x_new_df).to_numpy()[0, 0]) + + return Bunch( + x_new=pd.Series(x_new_df.to_numpy().ravel(), index=self._feature_names, name="x_new"), + predictive_score=tau_p, + y_hat=y_hat_value, + orthogonal_space_basis=self.orthogonal_loadings_.copy(), + orthogonal_space_dimension=int(self.n_orthogonal_components), + ) + + # ------------------------------------------------------------------ # + def _prepare(self, X: DataMatrix) -> tuple[np.ndarray, pd.Index]: + """Validate new X, capture its index, and move it to the scaled space.""" + check_is_fitted(self, "predictive_weights_") + index: pd.Index | None = X.index if isinstance(X, pd.DataFrame) else None + X_arr = validate_data( + self, X, reset=False, accept_sparse=False, dtype="numeric", ensure_all_finite=True + ) + if index is None: + index = pd.RangeIndex(X_arr.shape[0]) + X_df = pd.DataFrame(X_arr, index=index, columns=self._feature_names) + if self._x_scaler is not None: + X_df = self._x_scaler.transform(X_df) + return X_df.to_numpy(dtype=float), index + + +def _as_scalar(value: float | np.ndarray | pd.Series | pd.DataFrame) -> float | np.ndarray: + """Extract the underlying scalar/array from assorted single-response inputs.""" + if isinstance(value, pd.DataFrame): + return value.to_numpy().ravel() + if isinstance(value, pd.Series): + return value.to_numpy().ravel() + return value diff --git a/src/process_improve/multivariate/methods.py b/src/process_improve/multivariate/methods.py index a342164..d90dfcc 100644 --- a/src/process_improve/multivariate/methods.py +++ b/src/process_improve/multivariate/methods.py @@ -48,6 +48,7 @@ ssq, terminate_check, ) +from ._opls import OPLS from ._pca import PCA from ._pls import PLS from ._preprocessing import MCUVScaler, center, scale @@ -68,6 +69,7 @@ __all__ = [ "MBPCA", "MBPLS", + "OPLS", "PCA", "PLS", "REFERENCE_LINE_COLOR", diff --git a/tests/test_multivariate_opls.py b/tests/test_multivariate_opls.py new file mode 100644 index 0000000..d4210d2 --- /dev/null +++ b/tests/test_multivariate_opls.py @@ -0,0 +1,216 @@ +# (c) Kevin Dunn, 2010-2026. MIT License. +"""Tests for the O-PLS estimator and the equivalence with PLS. + +The central results checked here come from García-Carrión et al., "On the +equivalence between null space and orthogonal space in latent variable +regression modeling", Journal of Chemometrics, 39 (2025): e70057. An +A-component PLS model and an O-PLS(1; A - 1) model give identical predictions, +and the O-PLS orthogonal space is the same linear space as the PLS null space. +""" + +from __future__ import annotations + +import pathlib + +import numpy as np +import pandas as pd +import pytest + +from process_improve.multivariate.methods import OPLS, PLS, MCUVScaler + + +@pytest.fixture +def synthetic_single_response() -> tuple[pd.DataFrame, pd.DataFrame]: + """Build a correlated-X, single-response dataset with a known coefficient vector.""" + rng = np.random.default_rng(7) + n_samples, n_features = 60, 6 + base = rng.standard_normal((n_samples, 2)) + loadings = rng.standard_normal((2, n_features - 2)) + tail = base @ loadings + 0.1 * rng.standard_normal((n_samples, n_features - 2)) + X = pd.DataFrame(np.hstack([base, tail]), columns=[f"x{i}" for i in range(n_features)]) + beta = rng.standard_normal((n_features, 1)) + y = pd.DataFrame(X.to_numpy() @ beta + 0.2 * rng.standard_normal((n_samples, 1)), columns=["y"]) + return X, y + + +@pytest.fixture +def cheese() -> tuple[pd.DataFrame, pd.DataFrame]: + """Cheddar-cheese: 30 samples, 3 predictors, single response (Taste).""" + folder = pathlib.Path(__file__).parents[1] / "src" / "process_improve" / "datasets" / "multivariate" + path = folder / "cheddar-cheese.csv" + if not path.exists(): + pytest.skip("cheddar-cheese.csv fixture not present") + data = pd.read_csv(path) + return data[["Acetic", "H2S", "Lactic"]], data[["Taste"]] + + +@pytest.mark.parametrize("n_total", [2, 3, 4]) +def test_opls_predictions_match_pls(synthetic_single_response: tuple[pd.DataFrame, pd.DataFrame], n_total: int) -> None: + """O-PLS(1; A-1) reproduces the predictions of an A-component PLS model.""" + X, y = synthetic_single_response + pls = PLS(n_components=n_total).fit(X, y) + opls = OPLS(n_orthogonal_components=n_total - 1).fit(X, y) + assert opls.predict(X).to_numpy() == pytest.approx(pls.predict(X).to_numpy(), abs=1e-8) + assert opls.beta_coefficients_.to_numpy() == pytest.approx(pls.beta_coefficients_.to_numpy(), abs=1e-8) + + +def test_opls_orthogonal_scores_are_y_orthogonal(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """Every orthogonal score is uncorrelated with the (centered) response.""" + X, y = cheese + opls = OPLS(n_orthogonal_components=1).fit(X, y) + y_centered = y.to_numpy() - y.to_numpy().mean() + projection = opls.orthogonal_scores_.to_numpy().T @ y_centered + assert np.abs(projection).max() == pytest.approx(0.0, abs=1e-8) + + +def test_opls_predict_round_trips_through_invert(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """invert() then predict() recovers the desired response (single division).""" + X, y = cheese + opls = OPLS(n_orthogonal_components=1).fit(X, y) + for target in (15.0, 25.0, 40.0): + result = opls.invert(target) + assert result.y_hat == pytest.approx(target, abs=1e-6) + y_back = float(opls.predict(result.x_new.to_frame().T).iloc[0, 0]) + assert y_back == pytest.approx(target, abs=1e-6) + assert result.orthogonal_space_dimension == 1 + + +def test_orthogonal_space_equals_null_space(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """The O-PLS orthogonal space and the PLS null space span the same subspace. + + This is the central result of García-Carrión et al. (2025): reconstruct + both bases into the input space and confirm they point the same way. + """ + X, y = cheese + pls = PLS(n_components=2).fit(X, y) + opls = OPLS(n_orthogonal_components=1).fit(X, y) + + target = 25.0 + ns_scores = pls.invert(target).null_space_basis.to_numpy() # (A, A-1) in score space + ns_input = ns_scores.T @ pls.x_loadings_.to_numpy().T # (A-1, K) reconstructed to input space + os_input = opls.invert(target).orthogonal_space_basis.to_numpy().T # (Ao, K) + + # 1-D subspaces here: compare directions via absolute cosine. + ns_dir = ns_input.ravel() + os_dir = os_input.ravel() + cosine = np.abs(ns_dir @ os_dir) / (np.linalg.norm(ns_dir) * np.linalg.norm(os_dir)) + assert cosine == pytest.approx(1.0, abs=1e-8) + + +def test_opls_correct_removes_orthogonal_variation(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """The corrected X has no remaining projection on the orthogonal weights.""" + X, y = cheese + opls = OPLS(n_orthogonal_components=1).fit(X, y) + x_filtered = opls.correct(X).to_numpy() + # No systematic variation left along the orthogonal weight directions. + residual_on_ortho = x_filtered @ opls.orthogonal_weights_.to_numpy() + assert np.abs(residual_on_ortho).max() == pytest.approx(0.0, abs=1e-8) + + +def test_opls_scores_shape_and_labels(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """scores_ is [predictive | orthogonal] with informative labels.""" + X, y = cheese + opls = OPLS(n_orthogonal_components=2).fit(X, y) + assert list(opls.scores_.columns) == ["t_predictive", "t_orthogonal_1", "t_orthogonal_2"] + assert opls.scores_.shape == (X.shape[0], 3) + assert opls.n_components == 3 + + +def test_opls_transform_matches_predictive_scores(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """transform() on the training data reproduces the fitted predictive scores.""" + X, y = cheese + opls = OPLS(n_orthogonal_components=1).fit(X, y) + assert opls.transform(X).to_numpy() == pytest.approx(opls.predictive_scores_.to_numpy(), abs=1e-9) + + +def test_opls_rejects_multi_response() -> None: + """A two-column response is rejected with a clear message.""" + rng = np.random.default_rng(1) + X = pd.DataFrame(rng.standard_normal((20, 4)), columns=[f"x{i}" for i in range(4)]) + Y = pd.DataFrame(rng.standard_normal((20, 2)), columns=["y0", "y1"]) + with pytest.raises(ValueError, match="single response"): + OPLS(n_orthogonal_components=1).fit(X, Y) + + +def test_opls_scale_false_matches_prescaled(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """scale=False on pre-scaled data equals scale=True on raw data (predictions).""" + X, y = cheese + opls_scaled = OPLS(n_orthogonal_components=1).fit(X, y) + Xs = MCUVScaler().fit_transform(X) + ys = MCUVScaler().fit_transform(y) + opls_prescaled = OPLS(n_orthogonal_components=1, scale=False).fit(Xs, ys) + # Predictions on the scaled model are on the scaled-y scale; compare the + # standardized predictive scores instead, which are scale-free. + assert opls_prescaled.predictive_scores_.to_numpy() == pytest.approx( + opls_scaled.predictive_scores_.to_numpy(), abs=1e-8 + ) + + +def test_opls_accepts_1d_and_array_inputs(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """A 1-D y (Series and ndarray) and a plain-array X are accepted.""" + X, y = cheese + y_series = y.iloc[:, 0] + from_frame = OPLS(n_orthogonal_components=1).fit(X, y).predict(X).to_numpy() + from_series = OPLS(n_orthogonal_components=1).fit(X, y_series).predict(X).to_numpy() + from_array = OPLS(n_orthogonal_components=1).fit(X.to_numpy(), y.to_numpy().ravel()).predict(X.to_numpy()) + assert from_series == pytest.approx(from_frame, abs=1e-8) + assert from_array.to_numpy() == pytest.approx(from_frame, abs=1e-8) + + +def test_opls_feature_names_out(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """get_feature_names_out returns the component labels.""" + X, y = cheese + opls = OPLS(n_orthogonal_components=1).fit(X, y) + assert list(opls.get_feature_names_out()) == ["t_predictive", "t_orthogonal_1"] + + +def test_opls_negative_components_raises(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """A negative orthogonal-component count is rejected.""" + X, y = cheese + with pytest.raises(ValueError, match=">= 0"): + OPLS(n_orthogonal_components=-1).fit(X, y) + + +def test_opls_zero_orthogonal_components(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """O-PLS(1; 0) is a one-component model whose invert has no orthogonal space.""" + X, y = cheese + pls = PLS(n_components=1).fit(X, y) + opls = OPLS(n_orthogonal_components=0).fit(X, y) + assert opls.predict(X).to_numpy() == pytest.approx(pls.predict(X).to_numpy(), abs=1e-8) + assert opls.invert(20.0).orthogonal_space_dimension == 0 + + +def test_opls_invert_input_forms_and_errors(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """Invert accepts scalar/Series/DataFrame and rejects bad targets.""" + X, y = cheese + opls = OPLS(n_orthogonal_components=1).fit(X, y) + base = opls.invert(25.0).x_new.to_numpy() + assert opls.invert(pd.Series([25.0])).x_new.to_numpy() == pytest.approx(base) + assert opls.invert(pd.DataFrame({"y": [25.0]})).x_new.to_numpy() == pytest.approx(base) + with pytest.raises(ValueError, match="single desired response"): + opls.invert(np.array([1.0, 2.0])) + with pytest.raises(ValueError, match="finite"): + opls.invert(np.inf) + + +def test_opls_reproduces_paper_case_study_1() -> None: + """On the paper's Case Study 1, O-PLS(1;1) and PLS(2) predictions agree.""" + X = pd.DataFrame( + [ + [5.43, 7.54, 125.64, 58.51, 50.49], + [5.43, 15.97, 126.20, 258.48, 74.44], + [99.23, 7.54, 9893.38, 59.29, 737.15], + [99.23, 15.97, 9765.16, 254.11, 1576.28], + [52.33, 11.76, 2787.64, 139.21, 538.76], + [52.33, 11.76, 2849.95, 135.67, 630.73], + ], + columns=["x1", "x2", "x3", "x4", "x5"], + ) + y = pd.DataFrame({"y": [61.85, 278.99, 307.89, 436.4, 266.08, 260.52]}) + pls = PLS(n_components=2).fit(X, y) + opls = OPLS(n_orthogonal_components=1).fit(X, y) + assert opls.predict(X).to_numpy() == pytest.approx(pls.predict(X).to_numpy(), abs=1e-6) + # Both inversions reach the same design for a desired response. + yy = 168.23 + pls_y_hat = float(pls.predict(pls.invert(yy).x_new.to_frame().T).iloc[0, 0]) + assert opls.invert(yy).y_hat == pytest.approx(pls_y_hat, abs=1e-6) diff --git a/tests/test_sklearn_compat.py b/tests/test_sklearn_compat.py index 3c81f62..a7bc55c 100644 --- a/tests/test_sklearn_compat.py +++ b/tests/test_sklearn_compat.py @@ -23,14 +23,14 @@ from sklearn.decomposition import PCA as SKLearnPCA # noqa: N811 - aliased to avoid collision with our PCA from sklearn.preprocessing import StandardScaler -from process_improve.multivariate.methods import MBPCA, MBPLS, PCA, PLS, TPLS, MCUVScaler +from process_improve.multivariate.methods import MBPCA, MBPLS, OPLS, PCA, PLS, TPLS, MCUVScaler # Concrete sklearn estimators we explicitly refuse to inherit from (their private, # version-specific attribute layout is exactly what ENG-07 decouples us from). _CONCRETE_SKLEARN_ESTIMATORS = (PLSRegression, PLSCanonical, PLSSVD, SKLearnPCA) -@pytest.mark.parametrize("estimator_cls", [PCA, PLS, TPLS, MBPLS, MBPCA]) +@pytest.mark.parametrize("estimator_cls", [PCA, PLS, TPLS, MBPLS, MBPCA, OPLS]) def test_estimators_do_not_inherit_concrete_sklearn(estimator_cls: type) -> None: """No estimator inherits a concrete sklearn estimator, but all keep BaseEstimator.""" mro = estimator_cls.__mro__ @@ -46,6 +46,16 @@ def test_estimators_keep_expected_sklearn_mixins() -> None: assert TransformerMixin in PLS.__mro__ assert RegressorMixin in MBPLS.__mro__ assert TransformerMixin in MBPCA.__mro__ + assert RegressorMixin in OPLS.__mro__ + assert TransformerMixin in OPLS.__mro__ + + +def test_opls_clone_and_params_round_trip() -> None: + """OPLS supports the sklearn estimator API: clone, get_params, set_params.""" + est = OPLS(n_orthogonal_components=2, scale=False) + assert clone(est).get_params() == est.get_params() + rebuilt = OPLS(n_orthogonal_components=1).set_params(n_orthogonal_components=2, scale=False) + assert rebuilt.get_params() == est.get_params() def test_pca_clone_and_params_round_trip() -> None: From f9ba1a2c86ebdd28e48c52a0d99b4fc211220d42 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 19:51:46 +0000 Subject: [PATCH 6/7] Cover OPLS unscaled paths; exclude typing.TYPE_CHECKING from coverage Adds a scale=False end-to-end test (predict/transform/invert) and a 2-D ndarray Y case, taking _opls.py to full patch coverage. Also broadens the .coveragerc TYPE_CHECKING exclude to match the 'if typing.TYPE_CHECKING:' form the codebase actually uses, so type-only import blocks no longer count as missed lines. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DsTypY4FcRbrXRxorb9raz --- .coveragerc | 4 +++- tests/test_multivariate_opls.py | 20 ++++++++++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/.coveragerc b/.coveragerc index a56de25..b1e5a72 100644 --- a/.coveragerc +++ b/.coveragerc @@ -26,8 +26,10 @@ exclude_lines = raise NotImplementedError # Don't track main entry points if __name__ == .__main__.: - # Don't complain about TYPE_CHECKING-only imports + # Don't complain about TYPE_CHECKING-only imports (both the bare + # ``if TYPE_CHECKING:`` and the ``if typing.TYPE_CHECKING:`` form). if TYPE_CHECKING: + if typing.TYPE_CHECKING: # ENG-13 (#295): the fallback branches for each gated optional import # cannot fire under CI (which installs all extras). Exclude them by # matching lines that the fallback bodies use exclusively. diff --git a/tests/test_multivariate_opls.py b/tests/test_multivariate_opls.py index d4210d2..6890aee 100644 --- a/tests/test_multivariate_opls.py +++ b/tests/test_multivariate_opls.py @@ -193,6 +193,26 @@ def test_opls_invert_input_forms_and_errors(cheese: tuple[pd.DataFrame, pd.DataF opls.invert(np.inf) +def test_opls_scale_false_end_to_end_and_2d_array_y(cheese: tuple[pd.DataFrame, pd.DataFrame]) -> None: + """Exercise the unscaled code paths: predict/transform/invert and a 2-D array Y.""" + X, y = cheese + Xs = MCUVScaler().fit_transform(X) + ys = MCUVScaler().fit_transform(y) + # Y supplied as a 2-D ndarray (neither DataFrame nor 1-D) exercises that branch. + opls = OPLS(n_orthogonal_components=1, scale=False).fit(Xs, ys.to_numpy()) + + # predict / transform on the unscaled model. + assert opls.predict(Xs).shape == (X.shape[0], 1) + assert opls.transform(Xs).to_numpy() == pytest.approx(opls.predictive_scores_.to_numpy(), abs=1e-9) + + # invert on the unscaled model round-trips on the scaled-y scale. + target_scaled = 0.5 + result = opls.invert(target_scaled) + assert result.y_hat == pytest.approx(target_scaled, abs=1e-6) + y_back = float(opls.predict(result.x_new.to_frame().T).iloc[0, 0]) + assert y_back == pytest.approx(target_scaled, abs=1e-6) + + def test_opls_reproduces_paper_case_study_1() -> None: """On the paper's Case Study 1, O-PLS(1;1) and PLS(2) predictions agree.""" X = pd.DataFrame( From a2550cd7a7f0d18d3d0fe9feed6e87f8b006a249 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 20:13:57 +0000 Subject: [PATCH 7/7] Test multi-response PLS inversion on the solvents dataset Adds a two-response design test: seven solvent properties predicting logP and Solubility, inverted toward a target (logP, Solubility). Confirms the null-space dimension is A - rank(Y), the design round-trips to both targets, and moving along the null space leaves both predictions fixed. Vendors solvents.csv as an offline fixture. This is the boundary of the null-space / orthogonal-space equivalence, which is proven for a single response only. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DsTypY4FcRbrXRxorb9raz --- .../datasets/multivariate/solvents.csv | 104 ++++++++++++++++++ tests/test_multivariate_inversion.py | 39 +++++++ 2 files changed, 143 insertions(+) create mode 100644 src/process_improve/datasets/multivariate/solvents.csv diff --git a/src/process_improve/datasets/multivariate/solvents.csv b/src/process_improve/datasets/multivariate/solvents.csv new file mode 100644 index 0000000..68cd9bb --- /dev/null +++ b/src/process_improve/datasets/multivariate/solvents.csv @@ -0,0 +1,104 @@ +Chemical,MeltingPoint,BoilingPoint,Dielectric,DipoleMoment,RefractiveIndex,ET30,Density,logP,Solubility +W01_WATER,0,100,78.3,5.9,1.333,63.1,0.9982,-1.38,1.774 +G02_Formamide,2.55,210.5,111,11.2,1.4475,56.6,1.1334,-1.51,1.401 +A03_1.2-Ethanediol__,-12.6,197.15,37.7,7.7,1.4318,56.3,1.1088,-1.36,1.06 +A04_Methanol,-97.7,64.5,32.66,5.7,1.3284,55.4,0.7914,-0.77,1.393 +G05_N-Methyl_formamid,-3.8,182.5,182.4,12.9,1.4319,54.1,1.01,-1.3,1 +A06_diethylenglycol,-7.8,245.7,31.69,7.7,1.4475,53.8,1.109,1,1 +A07_triethylenglycol,-4.3,288,23.69,10,1.4558,53.5,1.0682,-1.33,0.8521 +A08_2-Methoxyethanol,-85.1,124.6,16.93,6.8,1.4021,52.3,0.965,-0.77,1.103 +G10_N-Methylacetamid,30.6,206.7,191.3,14.2,1.4253,52,0.957,-1.05,1 +A11_Ethanol,-114.5,78.3,24.55,5.8,1.3614,51.9,0.785,-0.31,1.231 +A12_2-AminoEthanol,10.5,170.95,37.72,7.6,1.4545,51.8,1.018,-1.31,1.216 +F13_Acetic_acid,16.7,117.9,6.17,5.6,1.3719,51.7,1.0492,-0.17,1.002 +A15_Benzylalcohol,-15.3,205.45,13.1,5.5,1.5404,50.4,1.042,1.1,-0.992 +A14_1-Propanol,-126.2,97.15,20.45,5.5,1.3856,50.7,0.804,0.25,1.126 +A16_1-butanol,-88.6,117.7,17.51,5.8,1.3993,50.2,0.8098,0.88,0.125 +A19_2-methyl_1-propa,-108,107.9,17.93,6,1.3959,48.6,0.794,0.76,0.176 +A20_2-propanol,-88,82.2,19.92,5.5,1.3772,48.4,0.786,0.05,1.117 +A21_2_butanol,-114.7,99.5,16.56,5.5,1.3971,47.1,0.808,0.61,0.176 +A18_3-methyl_1-butan,-117.2,130.5,15.19,6.1,1.4072,49,0.8092,1.42,-0.518 +A22_cyclohexanol,25.15,161.1,15,6.2,1.4648,46.9,0.962,1.23,-0.755 +C23_4-methyl-1.3-dio,-54.5,241.7,64.92,16.5,1.4215,46.6,1.204,1,1 +A24_2-pentanol,1,119,13.71,5.5,1.4064,46.5,0.81,1.19,-0.612 +L25_Nitromethane,-28.55,101.2,35.94,11.9,1.3819,46.3,1.137,-0.35,0.235 +N27_Acetonitrile,-43.8,81.6,35.94,11.8,1.3441,45.6,0.7857,-0.34,1.282 +A26_3-pentanol,-75,115.3,13.35,5.5,1.4104,45.7,0.821,1.21,-0.234 +X28_dimethylsulfoxid,18.5,189,46.45,13.5,1.4793,45.1,1.101,-1.35,1 ++29_aniline,-6,184.4,6.71,5,1.5863,44.3,1.0217,0.9,-0.37 +X30_Sulfolane,28.45,287.3,43.3,16,1.4816,44,1.262,-0.77,1 +F31_acetic_anhydride,-73.1,140,20.7,9.4,1.3904,43.9,1.082,1,0.07 +A35_t-butanol,25.6,82.3,12.47,5.5,1.3877,43.3,0.789,0.35,-1.057 +G32_N-N-dimethylform,-60.4,153,36.71,10.8,1.4305,43.8,0.945,-1.01,1.112 +G33_N-N-dimethylacet,-20,166.1,37.78,12.4,1.4384,43.7,0.937,-0.77,1.11 +N34_propannitrile,-92.8,97.35,28.86,11.7,1.3658,43.7,0.782,0.16,1 +G37_1-methyl-pyrroli,-24.4,202,32.2,13.6,1.47,42.2,1.026,1,1 +K38_Acetone,-94.7,56.1,20.56,9,1.3587,42.2,0.79,-0.24,1.134 +L44_nitrobenzene,5.8,210.8,34.78,13.3,1.5562,41.2,1.204,1.85,-1.778 +N41_Benzonitrile,-12.75,191.1,25.2,13.4,1.5282,41.5,1.01,1.56,-1.013 ++40_1.2.Diaminoethan,11.3,116.9,12.9,6.3,1.4568,42,0.899,-2.04,1 +H42_1.2-Dichloroetha,-35.7,83.5,10.7,6.1,1.4448,38,1.235,1.48,-1.092 +A45_2-methyl-2-butan,-8.8,102,5.78,5.7,1.405,41.1,0.806,0.89,0.039 +K43_2-butanone,-86.7,79.6,18.51,9.2,1.3788,41.3,0.805,0.29,0.494 +K52_Acetophenone,19.6,202,17.39,9.8,1.5342,40.6,1.0281,1.58,-1.553 +H51_dichlormethane,-94.9,39.6,8.93,5.2,1.4242,40.7,1.33,1.25,-0.699 +U47_1.1.3.3.tetramet,-1.2,175.2,23.6,11.7,1.4493,41,0.969,0.19,1 +G49_HMPT,7.2,233,29.6,18.5,1.4588,40.9,1.024,0.28,1 +K55_cyclohexanone,-32.1,155.65,16.1,10.3,1.451,39.8,0.9478,0.81,-1.236 ++53_pyridine,-41.55,115.25,12.91,7.9,1.5102,40.5,0.982,0.65,1.094 +E54_Methylacetate,-98.05,56.9,6.68,5.7,1.3614,40,0.933,0.18,0.486 +K56_4-methyl-2-penta,-84.7,117.4,13.11,2.7,1.3958,39.4,0.7978,1.31,-0.871 +H57_1.1.-dichloroeth,-97,57.3,10,6.1,1.4164,39.4,1.176,1.79,-1.259 ++58_Quinoline,-14.85,237.1,8.95,7.3,1.6273,39.4,1.093,2.03,1 +K59_3-pentanone,-39,102,17,9.4,1.3923,39.3,0.8138,0.91,-0.27 +H60_chloroform,-63.15,61.2,4.81,3.8,1.4459,39.1,1.48,1.97,-1.246 +T62_Trietylenglycold,-45,216,7.5,1,1.4224,38.9,1,1,1 +T64_dietylenglycoldi,-64,159.8,5.8,6.6,1.4078,38.6,1,1,1 +T65_1.2-Dimethoxyeth,-69,84.5,7.2,5.7,1.3796,38.2,0.8629,-0.21,1 +H67_1.2.Dichlorobenz,-17,180.5,9.93,7.1,1.5515,38,1.305,3.38,-3.029 +E66_ethylacetate,-83.55,77.1,6.02,6.1,1.3724,38.1,0.9,0.73,-0.041 +H73_fluorobenzene,-42.2,84.7,5.42,4.9,1.4684,37,1.023,2.27,-1.797 +H78_Iodobenzene,-31.35,188.3,4.49,4.7,1.62,36.2,1.831,3.25,-2.847 +H75_Chlorobenzene,-45.6,131.7,5.62,5.4,1.5248,36.8,1.106,2.84,-2.45 +H76_Bromobenzene,-30.8,155.9,5.4,5.2,1.5568,36.6,1.495,2.99,-2.57 +T70_tetrahydrofuran,-108.4,66,7.58,5.8,1.4072,37.4,0.889,0.46,1.142 +T71_methoxybenzene,-37.5,153.6,4.33,4.2,1.517,37.1,0.996,2.11,-1.917 +T77_Ethoxybenzene,-29.5,169.8,4.22,4.5,1.5074,36.6,0.967,2.51,-2.332 +H79_1.1.1-Trichloret,-30.4,74.1,7.25,5.7,1.438,36.2,1.339,2.49,-2.179 +T80_1.4-Dioxane,11.8,101.3,2.21,1.5,1.4224,36,1.034,-0.27,1.07 +H81_Trichlorethene,-86.4,87.2,3.42,2.7,1.4773,35.9,1.464,2.42,-1.995 ++83_Piperidine,-10.5,106.2,5.8,4,1.4525,35.5,0.861,0.84,1.005 +T85_Diphenylether,26.9,258.1,3.69,3.9,1.5763,35.3,1.075,4.21,-3.947 +T86_Diethylether,-116.3,34.4,4.2,3.8,1.3524,34.5,0.714,0.89,0.076 +D87_Benzene,5.5,80.1,2.27,0,1.5011,34.3,0.8786,2.13,-1.559 +T73_Diisopropylether,-85.5,68.3,3.88,4.2,1.3681,34,0.7241,1.52,-1.982 +D89_toluene,-95,110.6,2.38,1,1.4969,33.9,0.867,2.73,-1.747 +T75_Din-butylether,-95.2,140.3,3.08,3.94,1.3992,33,0.7689,3.21,-1.941 ++94_triethylamine,-114.7,88.9,2.42,2.9,1.401,32.1,0.7275,1.45,0.224 +D77_1.3.5-trimethylb,-44.7,164.7,2.28,0,1.4994,32.9,0.865,3.42,-3.356 +S92_Carbondisulfide,-111.6,46.2,2.64,0,1.6275,32.8,1.263,1.94,-1.552 +H93_carbontetrachlor,-22.8,76.6,2.23,0,1.4602,32.4,1.59,2.83,-2.294 +H80_Tetracloroethyle,-22.4,121.2,2.3,0,1.5057,31.9,1.623,3.4,-2.895 +D100_Cyclohexane,6.7,80.7,2.02,0,1.4262,30.9,0.778,3.44,-3.069 +D98_n-Hexane,-95.3,68.7,1.88,0,1.3749,31,0.66,3.98,-3.359 +A09_tetraethylenegly,-6.2,327.3,19.7,10.8,1.4577,52.2,1.1285,1.38,1 +A17_1-Pentanol,-78.2,138,13.9,5.7,1.41,49.1,0.8144,1.56,-0.631 +U36_DMEU,8.2,225.5,37.6,13.6,1.4707,42.5,1,1,1 +U_39_DMPU,-20,230,36.12,14.1,1.4881,42.1,1,1,1 +K46_2-Pentanone,-76.9,102.3,15.38,9,1.3908,41.1,0.8089,0.91,-0.328 ++48_Morpholine,-4.8,128.9,7.42,5.2,1.4542,41,1.0005,-0.86,1.06 +K50_3-Methyl2-butano,-92,94.2,15.87,9.2,1.388,40.9,0.8051,0.56,-0.164 +K61_3.3-Dimethyl-2-b,-49.8,106.3,13.1,9.3,1.3952,39,0.8012,0.85,-0.706 +K63_2.4-Dimethyl-3-p,-69,125.25,17.2,9.1,1.3999,38.7,0.8108,1.49,-1.299 +K68_2.6-Dimethyl4-he,-46,168.2,9.91,8.9,1.4122,38,0.8053,2.96,-1.731 +T69_dietylenglycolDi,-44.3,188.9,5.7,1,1.4115,37.5,0.9063,1,0.076 +C72_Diethylcarbonate,-43,126.8,2.82,3,1.3837,37,0.9752,1.21,-0.798 +H74_1.1-Dichlorethen,-122.6,31.6,4.82,4.3,1.4247,37,1.218,1.86,-1.666 +T82_t-Butyl-methylet,-108.6,55.2,4.5,4.1,1.369,35.5,0.7405,0.94,1 ++84_Diethylamine,-49.8,55.55,3.78,4,1.3846,35.4,0.7056,0.58,1 +T88_di-n-propylether,-123.2,90.1,3.39,4.4,1.3805,34,0.736,2.03,-1.388 +D90-p-Xylene,13.3,138.4,2.27,0,1.4958,33.1,0.8611,3.15,-2.806 ++95_Tri_n-butylamine,-70,214,1,2.6,1.4291,32.1,0.7771,4.56,-3.116 +D96_cis-decalin,-43,195.8,2.2,0,1.481,31.2,0.8965,5.08,-5.192 +D97_n-Heptane,-90.6,98.4,1.92,0,1.3876,31.1,0.68376,4.57,-4.046 +D99_n-Pentane,-129.7,36.1,1.84,0,1.3575,31,0.6262,3.39,-3.129 \ No newline at end of file diff --git a/tests/test_multivariate_inversion.py b/tests/test_multivariate_inversion.py index 4e9c454..93a2699 100644 --- a/tests/test_multivariate_inversion.py +++ b/tests/test_multivariate_inversion.py @@ -195,6 +195,45 @@ def test_invert_cheddar_cheese_holdout() -> None: assert y_back == pytest.approx(target, abs=1e-6) +def test_invert_solvents_two_response_design() -> None: + """Design a solvent with a target logP and Solubility (a two-response case). + + The solvents dataset has seven physical properties and two responses (logP + and Solubility). Inverting toward a desired (logP, Solubility) is a genuine + multi-response design: the null-space dimension is A - rank(Y), the design + round-trips to both targets, and moving along the null space leaves both + predictions unchanged. O-PLS's single-division inversion does not apply here, + which is why the null-space / orthogonal-space equivalence is single-response + only. + """ + folder = pathlib.Path(__file__).parents[1] / "src" / "process_improve" / "datasets" / "multivariate" + path = folder / "solvents.csv" + if not path.exists(): + pytest.skip("solvents.csv fixture not present") + solvents = pd.read_csv(path).dropna() + x_cols = ["MeltingPoint", "BoilingPoint", "Dielectric", "DipoleMoment", "RefractiveIndex", "ET30", "Density"] + y_cols = ["logP", "Solubility"] + X, Y = solvents[x_cols], solvents[y_cols] + + n_components = 3 + model = PLS(n_components=n_components).fit(X, Y) + rank_y = np.linalg.matrix_rank((Y - Y.mean()).to_numpy()) + assert rank_y == 2 + + desired = pd.Series({"logP": 0.5, "Solubility": 0.0}) + result = model.invert(desired) + assert result.null_space_dimension == n_components - rank_y # 3 - 2 = 1 + + y_back = model.predict(result.x_new.to_frame().T) + assert y_back[y_cols].to_numpy().ravel() == pytest.approx(desired.to_numpy(), abs=1e-6) + + # Both target responses are unchanged when moving along the null space. + moved = model.invert(desired, null_space_coordinates=[2.0]) + y_moved = model.predict(moved.x_new.to_frame().T) + assert y_moved[y_cols].to_numpy().ravel() == pytest.approx(desired.to_numpy(), abs=1e-6) + assert not np.allclose(moved.x_new.to_numpy(), result.x_new.to_numpy()) + + def test_invert_multi_response_null_space_dimension() -> None: """On the real LDPE data the null-space dimension is A - rank(Y).