diff --git a/.coveragerc b/.coveragerc index a56de250..b1e5a72a 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/CHANGELOG.md b/CHANGELOG.md index a4a9826e..5c243b65 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,32 @@ those changes. ## [Unreleased] +## [1.59.0] - 2026-07-23 + +### 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 + 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 +2625,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 a5445884..fca9fa58 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 39969999..4658d89b 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" @@ -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", @@ -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 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 00000000..bc40a737 --- /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/src/process_improve/datasets/multivariate/solvents.csv b/src/process_improve/datasets/multivariate/solvents.csv new file mode 100644 index 00000000..68cd9bb8 --- /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/src/process_improve/multivariate/__init__.py b/src/process_improve/multivariate/__init__.py index 33a2a838..0ba9fbd6 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 00000000..d4b53d01 --- /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/_pls.py b/src/process_improve/multivariate/_pls.py index ac14d85a..f6045d4c 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, diff --git a/src/process_improve/multivariate/methods.py b/src/process_improve/multivariate/methods.py index a342164c..d90dfcc8 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_inversion.py b/tests/test_multivariate_inversion.py new file mode 100644 index 00000000..93a2699b --- /dev/null +++ b/tests/test_multivariate_inversion.py @@ -0,0 +1,261 @@ +# (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_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 + 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_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). + + 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) diff --git a/tests/test_multivariate_opls.py b/tests/test_multivariate_opls.py new file mode 100644 index 00000000..6890aee3 --- /dev/null +++ b/tests/test_multivariate_opls.py @@ -0,0 +1,236 @@ +# (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_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( + [ + [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 3c81f62b..a7bc55cf 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: