Merge pull request #11 from HauserGroup/actual-actual-o2pls

Actual actual o2pls

Author: jsture

CHANGES.md

@@ -56,6 +56,12 @@ and default-value changes will be documented here.
 
 ### Added
 
+- `O2PLS`, a dense two-block O2PLS estimator with X/Y preprocessing, sequential
+  X- and Y-orthogonal filtering, final joint-subspace re-estimation, and
+  bidirectional `predict` / `predict_x` methods. The v1 implementation is dense
+  only and exposes `coef_filtered_` for scaled, X-filtered inputs rather than a
+  raw-space `coef_` alias.
+
 - `OPLS.coef_raw_` / `OPLS.intercept_raw_`: linear coefficients on the original raw
   input feature space, collapsing scaling, the orthogonal filter and the predictive
   PLS into one map, so `X @ coef_raw_.T + intercept_raw_` reproduces `predict(X)`.

docs/api/o2pls.md

@@ -0,0 +1,3 @@
+# O2PLS
+
+::: scikit_opls.O2PLS

docs/index.md

@@ -11,6 +11,8 @@ cleaned `X`. With `n_orthogonal=0` it reduces exactly to `PLSRegression`.
 ## Highlights
 
 - [`OPLS`](api/opls.md) — regressor and supervised transformer.
+- [`O2PLS`](api/o2pls.md) — two-block joint/orthogonal decomposition with
+  bidirectional prediction.
 - Cross-validated `n_orthogonal` selection via scikit-learn's `GridSearchCV`
   (see [Quickstart](quickstart.md)).
 - [`OPLSDA`](api/opls_da.md) — binary classifier composing `OPLS`.

docs/quickstart.md

@@ -17,6 +17,32 @@ model.transform_orthogonal(X) # orthogonal scores
 model.r2y_, model.rmse_       # training-fit summaries
 ```
 
+## Two-block O2PLS
+
+```python
+import numpy as np
+from scikit_opls import O2PLS
+
+rng = np.random.default_rng(1)
+Z = rng.normal(size=(80, 2))
+X = Z @ rng.normal(size=(2, 20)) + 0.1 * rng.normal(size=(80, 20))
+Y = Z @ rng.normal(size=(2, 10)) + 0.1 * rng.normal(size=(80, 10))
+
+model = O2PLS(n_components=2, n_x_orthogonal=1, n_y_orthogonal=1).fit(X, Y)
+model.predict(X)       # predicted joint Y structure, raw Y units
+model.predict_x(Y)     # predicted joint X structure, raw X units
+model.transform(X)     # X-side joint scores
+model.transform_y(Y)   # Y-side joint scores
+```
+
+`O2PLS.coef_filtered_` maps scaled, X-orthogonally-filtered `X` to scaled
+predicted `Y`; it is not a raw-space sklearn `coef_` alias. The joint loadings use
+the O2PLS orthonormal convention, so joint loadings equal joint weights.
+Requested orthogonal components may be truncated with a `ConvergenceWarning` when
+the enlarged preliminary subspace leaves no numerically resolvable block-specific
+residual variation, especially when the requested component total approaches a
+block's rank or feature dimension.
+
 ## Choosing `n_orthogonal` by cross-validation
 
 Use scikit-learn's `GridSearchCV` directly — `OPLS` has no path structure, so a

examples/o2pls_synthetic.py

@@ -0,0 +1,53 @@
+"""Synthetic two-block O2PLS example."""
+
+from __future__ import annotations
+
+import numpy as np
+from sklearn.metrics import r2_score
+
+from scikit_opls import O2PLS
+
+
+def _make_blocks(seed: int = 0):
+    """Generate two blocks with joint and block-specific latent structure."""
+    rng = np.random.default_rng(seed)
+    n_samples = 100
+    n_joint = 2
+    latent, _ = np.linalg.qr(rng.normal(size=(n_samples, 4)))
+    joint = latent[:, :n_joint] * np.array([5.0, 2.0])
+    x_specific = latent[:, 2:3] * 4.0
+    y_specific = latent[:, 3:4] * 3.0
+
+    x_basis, _ = np.linalg.qr(rng.normal(size=(12, 3)))
+    y_basis, _ = np.linalg.qr(rng.normal(size=(8, 3)))
+    x_joint = joint @ x_basis[:, :n_joint].T
+    y_joint = joint @ y_basis[:, :n_joint].T
+    X = x_joint + x_specific @ x_basis[:, 2:3].T
+    Y = y_joint + y_specific @ y_basis[:, 2:3].T
+    X += 0.01 * rng.normal(size=X.shape)
+    Y += 0.01 * rng.normal(size=Y.shape)
+    return X, Y, x_joint, y_joint
+
+
+def main() -> None:
+    """Fit O2PLS and report joint-structure recovery on synthetic data."""
+    X, Y, x_joint, y_joint = _make_blocks()
+    model = O2PLS(
+        n_components=2,
+        n_x_orthogonal=1,
+        n_y_orthogonal=1,
+        scale="center",
+    ).fit(X, Y)
+
+    print(f"X shape: {X.shape}")
+    print(f"Y shape: {Y.shape}")
+    print(f"X-orthogonal components: {model.n_x_orthogonal_}")
+    print(f"Y-orthogonal components: {model.n_y_orthogonal_}")
+    print(f"R2X joint / orthogonal: {model.r2x_:.3f} / {model.r2x_ortho_:.3f}")
+    print(f"R2Y joint / orthogonal: {model.r2y_:.3f} / {model.r2y_ortho_:.3f}")
+    print(f"Predicted joint Y R2: {r2_score(y_joint, model.predict(X)):.3f}")
+    print(f"Predicted joint X R2: {r2_score(x_joint, model.predict_x(Y)):.3f}")
+
+
+if __name__ == "__main__":
+    main()

src/scikit_opls/__init__.py

@@ -1,11 +1,13 @@
 """Scikit-OPLS: Orthogonal Projections to Latent Structures for scikit-learn."""
 
+from scikit_opls._o2pls import O2PLS
 from scikit_opls._opls import OPLS
 from scikit_opls._opls_da import OPLSDA
 
 __version__ = "0.1.0"
 
 __all__ = [
+    "O2PLS",
     "OPLS",
     "OPLSDA",
     "__version__",