From 2415baf051b5a7a9960435a5ce3bd1c8b61edabd Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 20:47:14 +0100
Subject: [PATCH 01/15] ref: remove redundant sign-correction in bootstrap

orthogonal_procrustes returns an unconstrained orthogonal matrix
(reflections allowed), so the per-component sign step never fired
(0 flips across 2000+ bootstraps up to 200x150 columns). Output is
bit-for-bit identical; the Procrustes alignment alone resolves sign.
---
 plsdo/core.py | 13 +++++--------
 1 file changed, 5 insertions(+), 8 deletions(-)

diff --git a/plsdo/core.py b/plsdo/core.py
index 0b29a72..56b7e37 100644
--- a/plsdo/core.py
+++ b/plsdo/core.py
@@ -137,14 +137,11 @@ def bootstrap(self, n_bootstraps: int = 10000) -> None:
             aligned_u_load = boot_u_load @ Q
             aligned_vt_load = Q.T @ boot_vt_load
 
-            # Sign correction
-            signs = np.sign(
-                np.sum(aligned_vt_load * self.vt_loadings, axis=1, keepdims=True)
-            )
-            signs[signs == 0] = 1.0
-
-            u_distribution.append(aligned_u_load * signs.T)
-            vt_distribution.append(aligned_vt_load * signs)
+            # No separate sign correction: orthogonal_procrustes returns an
+            # unconstrained orthogonal matrix (reflections allowed), so the
+            # alignment above already resolves each component's arbitrary sign.
+            u_distribution.append(aligned_u_load)
+            vt_distribution.append(aligned_vt_load)
 
         self.u_se = np.std(np.stack(u_distribution, axis=2), axis=2, ddof=1)
         self.vt_se = np.std(np.stack(vt_distribution, axis=2), axis=2, ddof=1)

From 0ae472e6d85b68ae596bebd324df56517dfc891a Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 20:52:17 +0100
Subject: [PATCH 02/15] test: add known-answer correctness tests for the PLS
 engine

Plant a known rank-1 structure and a near-degenerate two-component
structure and assert the engine recovers them: top-LV direction,
spectrum dominance, permutation significance in both directions, and
bootstrap-ratio reliability incl. the Procrustes alignment. Replaces
the tautological test_sign_consistency and fixes the vacuous 0.001
threshold in test_random_data_not_significant.

Mutation-checked: n-1->n divisor, no-op permutation, and Q=eye each
break an assertion.
---
 tests/test_core.py | 169 +++++++++++++++++++++++++++++++++++++++------
 1 file changed, 147 insertions(+), 22 deletions(-)

diff --git a/tests/test_core.py b/tests/test_core.py
index b3a0994..ecfdc08 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -133,7 +133,10 @@ def test_random_data_not_significant(self):
         model.fit()
         model.permutation_test(n_perms=500)
 
-        assert np.all(model.p_values > 0.001)
+        # > 0.05 is the meaningful non-significance threshold. The previous
+        # 0.001 floor was vacuous: the Phipson & Smyth p-value cannot fall
+        # below 1 / (n_perms + 1) ≈ 0.002, so it held even for a no-op test.
+        assert np.all(model.p_values > 0.05)
 
 
 class TestBootstrap:
@@ -192,27 +195,6 @@ def test_seed_reproducibility(self, x_array, y_array):
 
         np.testing.assert_array_equal(m1.u_bootstrap_ratios, m2.u_bootstrap_ratios)
 
-    def test_sign_consistency(self):
-        """Bootstrap loadings for the dominant feature should not flip sign."""
-        rng = np.random.default_rng(0)
-        n = 30
-        signal = rng.standard_normal(n)
-        X = np.column_stack([signal + 0.1 * rng.standard_normal(n) for _ in range(3)])
-        Y = np.column_stack([signal + 0.1 * rng.standard_normal(n) for _ in range(3)])
-        from plsdo.io import zscore_columns
-
-        X = zscore_columns(X)
-        Y = zscore_columns(Y)
-
-        model = PLS(X, Y, seed=42)
-        model.fit()
-        model.bootstrap(n_bootstraps=100)
-
-        dominant_loading_sign = np.sign(model.u_loadings[0, 0])
-        dominant_bsr_sign = np.sign(model.u_bootstrap_ratios[0, 0])
-        assert dominant_loading_sign == dominant_bsr_sign
-
-
 class TestBootstrapZscoreX:
     def test_zscore_x_false_does_not_alter_dummy_x(self):
         """Bootstrap with zscore_x=False must leave integer dummy X unchanged."""
@@ -304,3 +286,146 @@ def test_filters_on_significance_and_reliability(self):
         # Only LV1 should survive (significant + reliable on both sides)
         expected = np.array([True, False, False])
         np.testing.assert_array_equal(model.final_lvs, expected)
+
+
+class TestKnownAnswer:
+    """Feed the engine data with a planted, known structure and assert it is
+    recovered. The maths is the oracle: these are the automated *correctness*
+    signal for the notebook→class port, and they guard the SVD construction,
+    permutation sensitivity (both directions), and the bootstrap Procrustes
+    alignment against silent regression.
+    """
+
+    @staticmethod
+    def _planted_rank1():
+        """Rank-1 planted signal.
+
+        A single score vector ``t`` drives two X features (loadings ``a``) and
+        two Y features (loadings ``b``), each with equal magnitude and mixed
+        sign so that column z-scoring preserves the structure; the remaining
+        features are noise. After z-scoring the cross-covariance is ≈
+        ``outer(sign(a), sign(b)) · n/(n-1)``, a rank-1 matrix whose top
+        singular vectors recover ``a`` and ``b``.
+        """
+        from plsdo.io import zscore_columns
+
+        n = 60
+        a = np.array([1.0, -1.0, 0.0, 0.0])
+        b = np.array([1.0, -1.0, 0.0])
+        rng = np.random.default_rng(0)
+        t = rng.standard_normal(n)
+        X = np.outer(t, a) + 0.05 * rng.standard_normal((n, len(a)))
+        Y = np.outer(t, b) + 0.05 * rng.standard_normal((n, len(b)))
+        return zscore_columns(X), zscore_columns(Y), a, b, n
+
+    @staticmethod
+    def _near_degenerate():
+        """Two planted components with nearly equal singular values.
+
+        Bootstrap resamples rotate and swap the two near-degenerate components,
+        so the Procrustes alignment in ``bootstrap()`` is essential to keep the
+        loadings stable. Without it the standard errors inflate and the
+        bootstrap ratios collapse.
+        """
+        from plsdo.io import zscore_columns
+
+        n = 80
+        rng = np.random.default_rng(0)
+        t1 = rng.standard_normal(n)
+        t2 = rng.standard_normal(n)
+        t2 = t2 - (t2 @ t1) / (t1 @ t1) * t1  # orthogonalise t2 against t1
+        a1 = np.array([1.0, 1.0, 0, 0, 0, 0])
+        b1 = np.array([1.0, 1.0, 0, 0])
+        a2 = np.array([0, 0, 1.0, 1.0, 0, 0])
+        b2 = np.array([0, 0, 1.0, 1.0])
+        c2 = 0.95  # second component slightly weaker → near-degenerate
+        X = np.outer(t1, a1) + c2 * np.outer(t2, a2) + 0.3 * rng.standard_normal((n, 6))
+        Y = np.outer(t1, b1) + c2 * np.outer(t2, b2) + 0.3 * rng.standard_normal((n, 4))
+        return zscore_columns(X), zscore_columns(Y)
+
+    def test_top_lv_recovers_planted_directions(self):
+        """Top LV aligns with the planted loadings and dominates the spectrum."""
+        X, Y, a, b, n = self._planted_rank1()
+        model = PLS(X, Y, seed=42)
+        model.fit()
+
+        # Top singular value matches the analytic value 2·n/(n-1): two planted
+        # features each side, |correlation| ≈ 1, scaled by the n/(n-1) divisor.
+        # rel=0.01 distinguishes it from the 1/n divisor (which gives ≈ 2.0).
+        assert model.s[0] == pytest.approx(2 * n / (n - 1), rel=0.01)
+
+        # Second singular value is far smaller — the signal is rank-1.
+        assert model.s[1] / model.s[0] < 0.1
+
+        # The top singular vectors align with the planted directions (up to the
+        # arbitrary global sign of an SVD component).
+        ahat = a / np.linalg.norm(a)
+        bhat = b / np.linalg.norm(b)
+        cos_u = model.u[:, 0] @ ahat
+        cos_v = model.vt[0, :] @ bhat
+        assert abs(cos_u) > 0.95
+        assert abs(cos_v) > 0.95
+        # The X and Y sides share the same global sign (joint structure of the
+        # cross-covariance), so the two cosines have the same sign.
+        assert cos_u * cos_v > 0
+
+        # The two largest loadings fall on the planted features.
+        assert set(np.argsort(np.abs(model.u[:, 0]))[-2:]) == {0, 1}
+        assert set(np.argsort(np.abs(model.vt[0, :]))[-2:]) == {0, 1}
+
+    def test_permutation_significant_on_planted_not_on_random(self):
+        """Permutation p is small on planted data, large on random data."""
+        X, Y, _, _, _ = self._planted_rank1()
+        model = PLS(X, Y, seed=42)
+        model.fit()
+        model.permutation_test(n_perms=499)
+        # A no-op permutation (perm_order = arange) would leave the observed
+        # singular value in the null every time, forcing p = 1.0 here.
+        assert model.p_values[0] < 0.05
+
+        rng = np.random.default_rng(0)
+        from plsdo.io import zscore_columns
+
+        Xr = zscore_columns(rng.standard_normal((60, 4)))
+        Yr = zscore_columns(rng.standard_normal((60, 3)))
+        model_r = PLS(Xr, Yr, seed=7)
+        model_r.fit()
+        model_r.permutation_test(n_perms=499)
+        assert model_r.p_values[0] > 0.05
+
+    def test_bootstrap_ratios_reliable_on_planted_features(self):
+        """Dominant planted features are reliable (|BSR| > 1.96) on both sides
+        with the correct relative-sign structure."""
+        X, Y, _, _, _ = self._planted_rank1()
+        model = PLS(X, Y, seed=42)
+        model.fit()
+        model.bootstrap(n_bootstraps=500)
+
+        # Planted features 0 and 1 are reliable on both the X and Y sides.
+        assert abs(model.u_bootstrap_ratios[0, 0]) > 1.96
+        assert abs(model.u_bootstrap_ratios[1, 0]) > 1.96
+        assert abs(model.vt_bootstrap_ratios[0, 0]) > 1.96
+        assert abs(model.vt_bootstrap_ratios[0, 1]) > 1.96
+        # Noise features are not reliable.
+        assert abs(model.u_bootstrap_ratios[2, 0]) < 1.96
+        assert abs(model.u_bootstrap_ratios[3, 0]) < 1.96
+        # The two planted features carry opposite signs (matching a = [+, -]),
+        # a global-sign-invariant structural check.
+        assert np.sign(model.u_bootstrap_ratios[0, 0]) != np.sign(
+            model.u_bootstrap_ratios[1, 0]
+        )
+        assert np.sign(model.vt_bootstrap_ratios[0, 0]) != np.sign(
+            model.vt_bootstrap_ratios[0, 1]
+        )
+
+    def test_procrustes_keeps_degenerate_loadings_reliable(self):
+        """With near-degenerate components, Procrustes alignment keeps all
+        planted features reliable on the top LV. Without it (Q = I) the
+        rotating components inflate the standard errors and the ratios drop."""
+        X, Y = self._near_degenerate()
+        model = PLS(X, Y, seed=42)
+        model.fit()
+        model.bootstrap(n_bootstraps=500)
+
+        # All four planted features (two per component) are reliable on LV1.
+        assert np.all(np.abs(model.u_bootstrap_ratios[:4, 0]) > 1.96)

From 82f883c99d8013dcffcdb026d60cc5d3ea3e318d Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 20:53:50 +0100
Subject: [PATCH 03/15] test: add PLS invariant checks

SVD reconstruction (U diag(s) Vt == xcorr), orthonormality of the
singular vectors, and bootstrap SE shrinking as injected SNR rises.
Mutation-checked: scaling u in _decompose breaks reconstruction and
orthonormality.
---
 tests/test_core.py | 54 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 54 insertions(+)

diff --git a/tests/test_core.py b/tests/test_core.py
index ecfdc08..d9ed43c 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -429,3 +429,57 @@ def test_procrustes_keeps_degenerate_loadings_reliable(self):
 
         # All four planted features (two per component) are reliable on LV1.
         assert np.all(np.abs(model.u_bootstrap_ratios[:4, 0]) > 1.96)
+
+
+class TestInvariants:
+    """Cheap properties any correct PLS decomposition must satisfy. No
+    reference values needed; these guard against a future ``_decompose`` (e.g.
+    a SparsePLS override) silently breaking the SVD contract.
+    """
+
+    @staticmethod
+    def _planted(eps, seed=0):
+        from plsdo.io import zscore_columns
+
+        n = 60
+        a = np.array([1.0, -1.0, 0, 0])
+        b = np.array([1.0, -1.0, 0])
+        rng = np.random.default_rng(seed)
+        t = rng.standard_normal(n)
+        X = np.outer(t, a) + eps * rng.standard_normal((n, len(a)))
+        Y = np.outer(t, b) + eps * rng.standard_normal((n, len(b)))
+        return zscore_columns(X), zscore_columns(Y)
+
+    def test_svd_reconstructs_cross_covariance(self, x_array, y_array):
+        from plsdo.io import zscore_columns
+
+        model = PLS(zscore_columns(x_array), zscore_columns(y_array))
+        model.fit()
+        reconstructed = model.u @ np.diag(model.s) @ model.vt
+        np.testing.assert_allclose(reconstructed, model.xcorr, atol=1e-12)
+
+    def test_singular_vectors_orthonormal(self, x_array, y_array):
+        from plsdo.io import zscore_columns
+
+        model = PLS(zscore_columns(x_array), zscore_columns(y_array))
+        model.fit()
+        k = model.s.shape[0]
+        np.testing.assert_allclose(model.u.T @ model.u, np.eye(k), atol=1e-12)
+        np.testing.assert_allclose(model.vt @ model.vt.T, np.eye(k), atol=1e-12)
+
+    def test_bootstrap_se_shrinks_as_snr_rises(self):
+        """Stronger planted signal ⇒ smaller bootstrap SE / larger BSR on the
+        dominant feature."""
+        X_hi, Y_hi = self._planted(eps=0.1)
+        X_lo, Y_lo = self._planted(eps=0.8)
+
+        m_hi = PLS(X_hi, Y_hi, seed=42)
+        m_hi.fit()
+        m_hi.bootstrap(n_bootstraps=400)
+
+        m_lo = PLS(X_lo, Y_lo, seed=42)
+        m_lo.fit()
+        m_lo.bootstrap(n_bootstraps=400)
+
+        assert m_hi.u_se[0, 0] < m_lo.u_se[0, 0]
+        assert abs(m_hi.u_bootstrap_ratios[0, 0]) > abs(m_lo.u_bootstrap_ratios[0, 0])

From 94ed76b95bb1f04c1a9b6ef228922d8e959d850c Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 20:55:35 +0100
Subject: [PATCH 04/15] enh: pin component signs deterministically in fit()

A PLS component's global sign is arbitrary and np.linalg.svd can
return a different one across BLAS builds, flipping loadings, scores,
and bootstrap ratios machine-to-machine. Flip each component so its
largest-magnitude X loading is positive, making all outputs
reproducible across machines. Guarded by a new test; not a semantic
change (sign carries no scientific meaning).
---
 plsdo/core.py      | 17 +++++++++++++++++
 tests/test_core.py | 14 ++++++++++++++
 2 files changed, 31 insertions(+)

diff --git a/plsdo/core.py b/plsdo/core.py
index 56b7e37..91ad9ef 100644
--- a/plsdo/core.py
+++ b/plsdo/core.py
@@ -49,6 +49,7 @@ def fit(self):
         """Run PLS: cross-covariance, SVD, loadings, and subject scores."""
         self.xcorr = self.X.T @ self.Y / (self.n_subjects - 1)
         self._decompose()
+        self._fix_component_signs()
         self.u_loadings = self.u * self.s[np.newaxis, :]
         self.vt_loadings = self.s[:, np.newaxis] * self.vt
         self.x_scores = self.X @ self.u
@@ -63,6 +64,22 @@ def _decompose(self):
         """
         self.u, self.s, self.vt = np.linalg.svd(self.xcorr, full_matrices=False)
 
+    def _fix_component_signs(self):
+        """Pin each component's arbitrary global sign deterministically.
+
+        A PLS component's sign is not scientifically meaningful, but
+        ``np.linalg.svd`` can return a different one across BLAS builds, which
+        would flip loadings, scores, and bootstrap ratios machine-to-machine.
+        Flip each component so its largest-magnitude X loading is positive,
+        making all package outputs reproducible across machines. ``u`` and
+        ``vt`` share a component's sign, so both are flipped together.
+        """
+        max_idx = np.argmax(np.abs(self.u), axis=0)
+        signs = np.sign(self.u[max_idx, np.arange(self.u.shape[1])])
+        signs[signs == 0] = 1.0
+        self.u = self.u * signs
+        self.vt = self.vt * signs[:, np.newaxis]
+
     def _check_fitted(self):
         """Raise if fit() has not been called."""
         if not self._fitted:
diff --git a/tests/test_core.py b/tests/test_core.py
index d9ed43c..27a66ad 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -60,6 +60,20 @@ def test_scores_are_projections(self, x_array, y_array):
         np.testing.assert_allclose(model.x_scores, X @ model.u)
         np.testing.assert_allclose(model.y_scores, Y @ model.vt.T)
 
+    def test_sign_convention_largest_loading_positive(self, x_array, y_array):
+        """Each component is sign-fixed so its largest-magnitude X loading is
+        positive. A PLS component's global sign is arbitrary and can flip
+        across BLAS builds; pinning it makes outputs reproducible across
+        machines."""
+        from plsdo.io import zscore_columns
+
+        model = PLS(zscore_columns(x_array), zscore_columns(y_array))
+        model.fit()
+
+        for i in range(model.s.shape[0]):
+            col = model.u_loadings[:, i]
+            assert col[np.argmax(np.abs(col))] > 0
+
     def test_singular_values_descending(self, x_array, y_array):
         from plsdo.io import zscore_columns
 

From 299d10459c5da7e2ce7378cd02298101283454e4 Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 21:00:01 +0100
Subject: [PATCH 05/15] test: snapshot pipeline outputs as a drift guard

Lock the result-defining numeric outputs of both run_pipeline methods
and cross_validate_pipeline on the committed synthetic data. Tiered
tolerances: tight for deterministic outputs (singular values, loadings,
scores), absolute for permutation p-values, looser for bootstrap
ratios. Regenerate deliberately by deleting tests/data/regression/.
Mutation-checked: the n-1->n divisor breaks the snapshot.
---
 .../regression/correlational_p_values.npy     | Bin 0 -> 160 bytes
 .../correlational_singular_values.npy         | Bin 0 -> 160 bytes
 .../correlational_subject_scores.npy          | Bin 0 -> 320 bytes
 .../correlational_x_bootstrap_ratios.npy      | Bin 0 -> 288 bytes
 .../regression/correlational_x_loadings.npy   | Bin 0 -> 288 bytes
 .../correlational_y_bootstrap_ratios.npy      | Bin 0 -> 256 bytes
 .../regression/correlational_y_loadings.npy   | Bin 0 -> 256 bytes
 tests/data/regression/cv_mean_accuracy.npy    | Bin 0 -> 136 bytes
 .../data/regression/cv_mean_null_accuracy.npy | Bin 0 -> 136 bytes
 .../regression/discriminatory_p_values.npy    | Bin 0 -> 152 bytes
 .../discriminatory_singular_values.npy        | Bin 0 -> 152 bytes
 .../discriminatory_subject_scores.npy         | Bin 0 -> 320 bytes
 .../discriminatory_x_bootstrap_ratios.npy     | Bin 0 -> 200 bytes
 .../regression/discriminatory_x_loadings.npy  | Bin 0 -> 200 bytes
 .../discriminatory_y_bootstrap_ratios.npy     | Bin 0 -> 224 bytes
 .../regression/discriminatory_y_loadings.npy  | Bin 0 -> 224 bytes
 tests/test_regression.py                      | 154 ++++++++++++++++++
 17 files changed, 154 insertions(+)
 create mode 100644 tests/data/regression/correlational_p_values.npy
 create mode 100644 tests/data/regression/correlational_singular_values.npy
 create mode 100644 tests/data/regression/correlational_subject_scores.npy
 create mode 100644 tests/data/regression/correlational_x_bootstrap_ratios.npy
 create mode 100644 tests/data/regression/correlational_x_loadings.npy
 create mode 100644 tests/data/regression/correlational_y_bootstrap_ratios.npy
 create mode 100644 tests/data/regression/correlational_y_loadings.npy
 create mode 100644 tests/data/regression/cv_mean_accuracy.npy
 create mode 100644 tests/data/regression/cv_mean_null_accuracy.npy
 create mode 100644 tests/data/regression/discriminatory_p_values.npy
 create mode 100644 tests/data/regression/discriminatory_singular_values.npy
 create mode 100644 tests/data/regression/discriminatory_subject_scores.npy
 create mode 100644 tests/data/regression/discriminatory_x_bootstrap_ratios.npy
 create mode 100644 tests/data/regression/discriminatory_x_loadings.npy
 create mode 100644 tests/data/regression/discriminatory_y_bootstrap_ratios.npy
 create mode 100644 tests/data/regression/discriminatory_y_loadings.npy
 create mode 100644 tests/test_regression.py

diff --git a/tests/data/regression/correlational_p_values.npy b/tests/data/regression/correlational_p_values.npy
new file mode 100644
index 0000000000000000000000000000000000000000..eb26c44e41906ab858dfd6a5358a95c1f302d892
GIT binary patch
literal 160
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
zXCxM+0{I$-ItnJ5ItsN4WCO0)^IwmjPAIXjGgZ69JoT~tyrMGh{ZfbRpS}?T3E2Yx
D)H^EJ

literal 0
HcmV?d00001

diff --git a/tests/data/regression/correlational_singular_values.npy b/tests/data/regression/correlational_singular_values.npy
new file mode 100644
index 0000000000000000000000000000000000000000..6cb3e56a090e07bd69a50c875ebdf873820ce9d5
GIT binary patch
literal 160
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
zXCxM+0{I$-ItnJ5ItsN4WCN~;I&mff+QJUqPilj1dCjmlH)iJ#Fs`<53RS2ysGDd5
E0D&SVegFUf

literal 0
HcmV?d00001

diff --git a/tests/data/regression/correlational_subject_scores.npy b/tests/data/regression/correlational_subject_scores.npy
new file mode 100644
index 0000000000000000000000000000000000000000..0d5bd87e5628c2186dcceb0138502e5098abfe04
GIT binary patch
literal 320
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u
zBo?Fsxf+H>ItoUbItsN4WCJd5^HtN0%~%hF9njrZ@|5|&<P|DUcn-53V9MChkpG0`
zz?X}<>p73GIe0ya`+0gYqeJ7bD{Gg_b2|t&O)m1{WOZPBl&RhM@Z|nVsR!wCWoPXd
zeP5#=<M(9$jKdoruUyt+|48JMsCp~&fj^z`JFG+*4jiA~QY`$H?LfhKjpqhu84nce
zr)quaV0JLyCuF01=#TyB{|B~R+rZ(#_-gj=zZV!CHvf*|Tv4=U|L*S@uX#USw_mYs
T|J&`n_xC%`uzj^7>4-f5!VP?I

literal 0
HcmV?d00001

diff --git a/tests/data/regression/correlational_x_bootstrap_ratios.npy b/tests/data/regression/correlational_x_bootstrap_ratios.npy
new file mode 100644
index 0000000000000000000000000000000000000000..f256075eb48176f21c37f91fce5e5f318e6340ee
GIT binary patch
literal 288
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u
zBo?Fsxf-TA3MQI53bhL41Fkgj485I<-VQ=W8F?9tk`H7{w(H-WkmK-U?z8_qJChGM
zAB+^zc^u>r7Ni<&8S}*cO3U`UmM<^cuY5L3v99^u{)QC28%8_k+W(R0t1fllZU3hC
zjE9^1sr|V(<(CGgzp+33@;P_+rHl5GGY@(v_+7BysJ@cj!{ePjV|n6|-IG7t%oIEO
j@?NgQ-p}iWj;=p`XivaduS=ICRBTUtvD4ew!eI*lf3|9s

literal 0
HcmV?d00001

diff --git a/tests/data/regression/correlational_x_loadings.npy b/tests/data/regression/correlational_x_loadings.npy
new file mode 100644
index 0000000000000000000000000000000000000000..4e25fbbc7f2836911353c2a3048817cfa3b0f2de
GIT binary patch
literal 288
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u
zBo?Fsxf-TA3MQI53bhL411=+`2}Q3TF*wZXw$R#P!FWJOWP`ShKcj=gdX>xqXT}3}
zP0XLZ<79MD`)t~7|7DW>@^vv04Yif_^Zp2*bT{qZZ?|A`+VKyr_FpreUwIIlXur=-
z;Pe^sko`utPDF<`<=b<v@0``(=5Mch<<&J$J%9UYCugTAJ}j}<@1Ig+E)rwo)gEri
hp47O<c*5DT)U{fBHVTW#Y<k>dbL7E=i<`U(Z2&PLV?F=?

literal 0
HcmV?d00001

diff --git a/tests/data/regression/correlational_y_bootstrap_ratios.npy b/tests/data/regression/correlational_y_bootstrap_ratios.npy
new file mode 100644
index 0000000000000000000000000000000000000000..ec133852d13d0cedbcdee4695bd6cd03299de912
GIT binary patch
literal 256
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u
zBo?Fsxf&)q3MQI53bhL41FpCKPx&6#&p8ktGJj^)lWd1qvhx@_SGgYO4VF$<Sdry0
zO+4QJP4<obQIfYU|8<_Uf1>95K|=SJeO<MFwrTVUduQFhOPYG$+8>f?)zmup!rr1h
z;?alB<Mx}k+`Ra|=Cb|SEjsI-&1BkpU&wWtsi?SZ+Sjs)O@E*6kx*6IurNV!F90l+
BTe<)M

literal 0
HcmV?d00001

diff --git a/tests/data/regression/correlational_y_loadings.npy b/tests/data/regression/correlational_y_loadings.npy
new file mode 100644
index 0000000000000000000000000000000000000000..c328d242e001f695037f5550146e19a1b979cde5
GIT binary patch
literal 256
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u
zBo?Fsxf&)q3MQI53bhL41Flc5yCzMJVLp(pP}sF5jM?E;^G#d3*-QugpQ%;<jAC|J
z<kuJC(BHD(!L!tjy`aQi*<tpx+{aVwgBWCHEk9CZe>@>cZ}sG2`}7NsgLV}b+CRK&
zz_&&q!2W~zCT<V@2z%f23%)bX&D-NwoY`q`y~(EF9LEf$_x^isbFzmDr1k6p0AZw5
AK>z>%

literal 0
HcmV?d00001

diff --git a/tests/data/regression/cv_mean_accuracy.npy b/tests/data/regression/cv_mean_accuracy.npy
new file mode 100644
index 0000000000000000000000000000000000000000..651688ccb8c89d5540e2cb292522fcd7e5843cf9
GIT binary patch
literal 136
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
aXCxM+0{I$-I+{8PwF(pfE<rGOX%7H1CmiYk

literal 0
HcmV?d00001

diff --git a/tests/data/regression/cv_mean_null_accuracy.npy b/tests/data/regression/cv_mean_null_accuracy.npy
new file mode 100644
index 0000000000000000000000000000000000000000..67ce234f49ebaa361daa88380dad3bf2463b0579
GIT binary patch
literal 136
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
bXCxM+0{I$-I+{8PwF(pfuJ#TPxL^+eLPQ^$

literal 0
HcmV?d00001

diff --git a/tests/data/regression/discriminatory_p_values.npy b/tests/data/regression/discriminatory_p_values.npy
new file mode 100644
index 0000000000000000000000000000000000000000..0a6c4afde6b6a79d9ed6553e9d43d6500f5bfff3
GIT binary patch
literal 152
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
tXCxM+0{I$-Its>`ItsN4WCO0)^IwmjPAIXT4yKRTCvE?*+Hm4JdjO~#Dd+$I

literal 0
HcmV?d00001

diff --git a/tests/data/regression/discriminatory_singular_values.npy b/tests/data/regression/discriminatory_singular_values.npy
new file mode 100644
index 0000000000000000000000000000000000000000..a96868ea7bf042e6d5c4740434215a4b04ec653e
GIT binary patch
literal 152
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
wXCxM+0{I$-Its>`ItsN4WCJc?&&7dnmVCBnY1N(A?z7Q8sry^%vQ4XO0Cip`9smFU

literal 0
HcmV?d00001

diff --git a/tests/data/regression/discriminatory_subject_scores.npy b/tests/data/regression/discriminatory_subject_scores.npy
new file mode 100644
index 0000000000000000000000000000000000000000..806b44b55e4bf716be4f11e6c9eee308d6bf0797
GIT binary patch
literal 320
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u
zBo?Fsxf+H>ItoUbItsN4WCN~sojx_)bDr(TMT^K!mRrB%nLRFg-YKnK7x#euxag?A
zclYQtGasnqZhS5+#c;slz11Y!H*5!TTJ(dKoMSwY#Vg&!-No!+?<(K&<i;O+^9r>?
zOO|su)K19}Z@$XtU@}?zmpQ}Q{p~V~R@bk%ZXeF5yQ=raz5Q~lV%Ah}9I*!g4*zvj

literal 0
HcmV?d00001

diff --git a/tests/data/regression/discriminatory_x_bootstrap_ratios.npy b/tests/data/regression/discriminatory_x_bootstrap_ratios.npy
new file mode 100644
index 0000000000000000000000000000000000000000..ec46f1233c1987c1dc6388df67896553a2bc7d26
GIT binary patch
literal 200
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u
zBo?Fsxf;eg3dWi`3bhL411^JgOEx{s5;<_`$}z12GesQgW~s0S<`(ZS-z|8{{L=sZ
vU$;xHJDvV{zmip_|AY^}><^}Go&Rqqx9!T6V!G$JS!}(hYO-E^_sa$VXF)(^

literal 0
HcmV?d00001

diff --git a/tests/data/regression/discriminatory_x_loadings.npy b/tests/data/regression/discriminatory_x_loadings.npy
new file mode 100644
index 0000000000000000000000000000000000000000..c516cc4272c1ebc56f21d7107f33a8dae3575461
GIT binary patch
literal 200
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u
zBo?Fsxf;eg3dWi`3bhL41Fl62_Zt2<^m>2q^XKW`&b+pNb)`Eb=TGqdP0UNL{_~r&
kU-m@a`zC`q`@cQye00EmjlG-4mB^o$=h=8cX)`Dd0ERh73;+NC

literal 0
HcmV?d00001

diff --git a/tests/data/regression/discriminatory_y_bootstrap_ratios.npy b/tests/data/regression/discriminatory_y_bootstrap_ratios.npy
new file mode 100644
index 0000000000000000000000000000000000000000..247136ea8bcab316047fccf27a4b3ffcbc3b4a80
GIT binary patch
literal 224
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u
zBo?Fsxf&)q3dWi`3bhL411|l;dknOhR1duObh)(VfU3hnwtGn-6^aMKT|VwVw?frH
zm-oR<@2-FLXFvRom^=TP{SBe)ou4JY?$7hl4BdS2%Kl{Lvq4i9{IfC2Nx5Hh@~2Hj
Upyz~~ls`7i&U$~QKa;Qp0HMK4y#N3J

literal 0
HcmV?d00001

diff --git a/tests/data/regression/discriminatory_y_loadings.npy b/tests/data/regression/discriminatory_y_loadings.npy
new file mode 100644
index 0000000000000000000000000000000000000000..fa0a503788706a190546c745fa7066f6c273adc8
GIT binary patch
literal 224
zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u
zBo?Fsxf&)q3dWi`3bhL41Fp1Gj!=frkN2OQvLI<A(-Zq5$0Z3W^B?WMyug6@<=4md
z%ckTqlrLRw-}5qt@6*R}`?}A6TZ=!e*#A8Fa?&E9(*34y8JSPco^7*HN!{pqNSDpm
UFtd&3kEYwKl`ec)6}HF*01j_VzW@LL

literal 0
HcmV?d00001

diff --git a/tests/test_regression.py b/tests/test_regression.py
new file mode 100644
index 0000000..4cfcc87
--- /dev/null
+++ b/tests/test_regression.py
@@ -0,0 +1,154 @@
+"""Regression snapshot guard against silent drift in pipeline outputs.
+
+This locks the result-defining numbers produced end-to-end by the pipelines on
+the committed synthetic data, so a future change that silently alters results
+(e.g. the 1/(n-1) divisor, the z-scoring ddof, the bootstrap-ratio definition,
+or the sign convention) is caught. It is a *drift* guard, not a correctness
+oracle — correctness rests on the known-answer tests and the cited literature.
+
+Tolerances are tiered by how reproducible each quantity is across machines:
+
+* deterministic outputs (singular values, loadings, subject scores) are pinned
+  tightly — they depend only on the SVD and the fixed sign convention;
+* permutation p-values use an absolute tolerance, since they are counts and a
+  near-tied null singular value can shift the count by one across BLAS builds;
+* bootstrap ratios use a looser tolerance, as they accumulate per-resample SVD
+  differences. A genuine code regression moves these far more than the
+  tolerance, so the guard still bites.
+
+Regenerating the snapshot is a deliberate, reviewed step (e.g. a major
+numpy/scipy bump): delete ``tests/data/regression/`` and run the suite once to
+recreate it, then commit. Never regenerate to paper over an unexplained diff.
+"""
+
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from plsdo.pipeline import cross_validate_pipeline, run_pipeline
+
+DATA_DIR = Path(__file__).parent / "data"
+REF_DIR = DATA_DIR / "regression"
+
+TIGHT = dict(rtol=1e-6, atol=1e-9)
+PVAL = dict(rtol=0.0, atol=0.02)
+LOOSE = dict(rtol=1e-2, atol=1e-2)
+SCALAR = dict(rtol=0.0, atol=0.05)
+
+
+def _load(path: Path, index_col=None) -> np.ndarray:
+    return pd.read_csv(path, index_col=index_col).to_numpy(dtype=float)
+
+
+def _snapshot(name: str, arr: np.ndarray, tol: dict) -> None:
+    """Compare ``arr`` against the committed snapshot, or create it if absent.
+
+    On a fresh ``tests/data/regression/`` (deliberate regeneration) the
+    snapshot is written and the assertion is skipped for that quantity; on all
+    later runs the committed snapshot is asserted against.
+    """
+    ref_path = REF_DIR / f"{name}.npy"
+    arr = np.asarray(arr)
+    if not ref_path.exists():
+        REF_DIR.mkdir(parents=True, exist_ok=True)
+        np.save(ref_path, arr)
+        return
+    np.testing.assert_allclose(
+        arr, np.load(ref_path), err_msg=f"regression drift in '{name}'", **tol
+    )
+
+
+def _run_pls(method: str, out: Path) -> Path:
+    kwargs = dict(
+        method=method,
+        y_path=DATA_DIR / "behaviour.csv",
+        demographics_path=DATA_DIR / "demographics.csv",
+        output_dir=out,
+        group_col="group",
+        subject_id="subject_id",
+        n_perms=200,
+        n_bootstraps=200,
+        seed=42,
+        img_format="png",
+        dpi=72,
+    )
+    if method == "correlational":
+        kwargs["x_path"] = DATA_DIR / "brain.csv"
+    run_pipeline(**kwargs)
+    return out / "data"
+
+
+@pytest.fixture(scope="module")
+def correlational_data(tmp_path_factory):
+    return _run_pls("correlational", tmp_path_factory.mktemp("corr"))
+
+
+@pytest.fixture(scope="module")
+def discriminatory_data(tmp_path_factory):
+    return _run_pls("discriminatory", tmp_path_factory.mktemp("disc"))
+
+
+@pytest.fixture(scope="module")
+def cv_data(tmp_path_factory):
+    out = tmp_path_factory.mktemp("cv")
+    cross_validate_pipeline(
+        y_path=DATA_DIR / "behaviour.csv",
+        demographics_path=DATA_DIR / "demographics.csv",
+        output_dir=out,
+        group_col="group",
+        subject_id="subject_id",
+        n_folds=3,
+        n_repeats=5,
+        n_permutations=20,
+        seed=42,
+        img_format="png",
+        dpi=72,
+    )
+    return out / "data"
+
+
+class TestPipelineSnapshot:
+    @pytest.mark.parametrize("method", ["correlational", "discriminatory"])
+    def test_deterministic_outputs(self, method, request):
+        data = request.getfixturevalue(f"{method}_data")
+        _snapshot(f"{method}_singular_values", _load(data / "singular_values.csv"), TIGHT)
+        _snapshot(
+            f"{method}_x_loadings", _load(data / "x_loadings.csv", index_col=0), TIGHT
+        )
+        _snapshot(
+            f"{method}_y_loadings", _load(data / "y_loadings.csv", index_col=0), TIGHT
+        )
+        _snapshot(
+            f"{method}_subject_scores",
+            _load(data / "subject_scores.csv", index_col=0),
+            TIGHT,
+        )
+
+    @pytest.mark.parametrize("method", ["correlational", "discriminatory"])
+    def test_permutation_pvalues(self, method, request):
+        data = request.getfixturevalue(f"{method}_data")
+        _snapshot(f"{method}_p_values", _load(data / "p_values.csv"), PVAL)
+
+    @pytest.mark.parametrize("method", ["correlational", "discriminatory"])
+    def test_bootstrap_ratios(self, method, request):
+        data = request.getfixturevalue(f"{method}_data")
+        _snapshot(
+            f"{method}_x_bootstrap_ratios",
+            _load(data / "x_bootstrap_ratios.csv", index_col=0),
+            LOOSE,
+        )
+        _snapshot(
+            f"{method}_y_bootstrap_ratios",
+            _load(data / "y_bootstrap_ratios.csv", index_col=0),
+            LOOSE,
+        )
+
+    def test_cross_validation(self, cv_data):
+        fold = pd.read_csv(cv_data / "cv_fold_results.csv")
+        null = pd.read_csv(cv_data / "cv_permutation_accuracies.csv")
+        _snapshot("cv_mean_accuracy", np.array([fold["accuracy"].mean()]), SCALAR)
+        _snapshot(
+            "cv_mean_null_accuracy", np.array([null["null_accuracy"].mean()]), SCALAR
+        )

From 2610f67d6f09c417a40ff859b1516d5930a21ce3 Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 21:02:04 +0100
Subject: [PATCH 06/15] enh: give a clear error when scikit-learn is missing

Wrap the sklearn imports in cross_validate.py so a missing optional
dependency raises a helpful 'install plsdo[cv]' message instead of a
raw ImportError. import plsdo stays sklearn-free. Guarded by a test
that simulates sklearn absence.
---
 plsdo/cross_validate.py      | 22 ++++++++++++++--------
 tests/test_cross_validate.py | 24 ++++++++++++++++++++++++
 2 files changed, 38 insertions(+), 8 deletions(-)

diff --git a/plsdo/cross_validate.py b/plsdo/cross_validate.py
index 3d0c30a..df394c4 100644
--- a/plsdo/cross_validate.py
+++ b/plsdo/cross_validate.py
@@ -2,14 +2,20 @@
 
 import numpy as np
 import pandas as pd
-from sklearn.cross_decomposition import PLSRegression
-from sklearn.metrics import (
-    accuracy_score,
-    balanced_accuracy_score,
-    confusion_matrix,
-)
-from sklearn.model_selection import RepeatedStratifiedKFold
-from sklearn.preprocessing import StandardScaler
+try:
+    from sklearn.cross_decomposition import PLSRegression
+    from sklearn.metrics import (
+        accuracy_score,
+        balanced_accuracy_score,
+        confusion_matrix,
+    )
+    from sklearn.model_selection import RepeatedStratifiedKFold
+    from sklearn.preprocessing import StandardScaler
+except ImportError as exc:
+    raise ImportError(
+        "scikit-learn is required for cross-validation but is not installed. "
+        "Install the optional dependency with: pip install 'plsdo[cv]'"
+    ) from exc
 
 from plsdo.io import corrected_pvalue
 
diff --git a/tests/test_cross_validate.py b/tests/test_cross_validate.py
index 53a2fbd..269fb44 100644
--- a/tests/test_cross_validate.py
+++ b/tests/test_cross_validate.py
@@ -117,3 +117,27 @@ def test_predicted_labels_within_group_range(self):
 
         results = run_cv(X, labels, n_splits=5, n_repeats=2, n_components=2, seed=42)
         assert set(results["pred_labels"]).issubset({0, 1, 2})
+
+
+class TestSklearnImportGuard:
+    def test_missing_sklearn_raises_helpful_error(self, monkeypatch):
+        """Importing cross_validate without scikit-learn points at plsdo[cv]."""
+        import builtins
+        import importlib
+        import sys
+
+        for mod in list(sys.modules):
+            if mod == "sklearn" or mod.startswith("sklearn.") or mod == "plsdo.cross_validate":
+                monkeypatch.delitem(sys.modules, mod, raising=False)
+
+        real_import = builtins.__import__
+
+        def fake_import(name, *args, **kwargs):
+            if name == "sklearn" or name.startswith("sklearn."):
+                raise ImportError("simulated missing scikit-learn")
+            return real_import(name, *args, **kwargs)
+
+        monkeypatch.setattr(builtins, "__import__", fake_import)
+
+        with pytest.raises(ImportError, match=r"plsdo\[cv\]"):
+            importlib.import_module("plsdo.cross_validate")

From 4befb4c90436193d3556bc42420a746fd7932b86 Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 21:16:39 +0100
Subject: [PATCH 07/15] enh: --bsr-threshold also controls LV survival

filter_lvs was hardcoded to 1.96 while --bsr-threshold only affected
which loadings were plotted, so the flag silently failed to change
which LVs survived. Pass the flag through so one knob means one thing;
default 1.96 keeps existing behaviour unchanged. Update the help and
docstring accordingly.
---
 plsdo/cli.py           |  5 +++--
 plsdo/pipeline.py      | 11 +++++++----
 tests/test_pipeline.py | 29 +++++++++++++++++++++++++++++
 3 files changed, 39 insertions(+), 6 deletions(-)

diff --git a/plsdo/cli.py b/plsdo/cli.py
index 7eecb02..7a2e2dd 100644
--- a/plsdo/cli.py
+++ b/plsdo/cli.py
@@ -81,8 +81,9 @@ def pls_main(argv=None):
         default=1.96,
         type=float,
         help=(
-            "Plot loading bars only for features with |bootstrap ratio| > "
-            "THRESHOLD (default: 1.96). Does not affect CSV outputs."
+            "Bootstrap-ratio reliability threshold (default: 1.96). Controls "
+            "both which latent variables survive and which loading bars are "
+            "plotted. Per-component CSV outputs are written regardless."
         ),
     )
     run_common.add_argument(
diff --git a/plsdo/pipeline.py b/plsdo/pipeline.py
index 8f20ffa..d75f85e 100644
--- a/plsdo/pipeline.py
+++ b/plsdo/pipeline.py
@@ -131,9 +131,12 @@ def run_pipeline(
     all_plots : bool
         If True, generate additional diagnostic plots.
     bsr_threshold : float
-        Plot loading bars only for features with |bootstrap ratio|
-        exceeding this threshold. Default 1.96 (≈ 95% CI under the
-        standard-normal approximation). CSV outputs are unaffected.
+        Bootstrap-ratio reliability threshold. Default 1.96 (≈ 95% CI under
+        the standard-normal approximation). Controls both which latent
+        variables survive ``filter_lvs`` (a surviving LV needs at least one
+        feature with |bootstrap ratio| > threshold on each side) and which
+        loading bars are plotted. The loading and bootstrap-ratio CSVs are
+        written for every component regardless.
     verbose_feature_limit : int, optional
         Maximum number of features before verbose plots (except scree)
         are skipped. Defaults to ``VERBOSE_FEATURE_LIMIT`` (100).
@@ -211,7 +214,7 @@ def run_pipeline(
     model.fit()
     model.permutation_test(n_perms=n_perms)
     model.bootstrap(n_bootstraps=n_bootstraps)
-    model.filter_lvs()
+    model.filter_lvs(bsr_threshold=bsr_threshold)
 
     # --- Save data CSVs ---
     _save_csv(
diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py
index 3e0826b..8cbf6ff 100644
--- a/tests/test_pipeline.py
+++ b/tests/test_pipeline.py
@@ -516,3 +516,32 @@ def test_facet_col_wrap_inert_with_facet_role_warns(
             calls = self._capture_boxstrip_calls(config, monkeypatch, tmp_path)
         assert calls[0]["col_wrap"] is None
         assert "facet_col_wrap" in caplog.text
+
+
+class TestBsrThresholdControlsSurvival:
+    """The --bsr-threshold flag controls LV survival, not only plotting."""
+
+    def _run(self, out, bsr_threshold):
+        run_pipeline(
+            method="correlational",
+            x_path=DATA_DIR / "brain.csv",
+            y_path=DATA_DIR / "behaviour.csv",
+            demographics_path=DATA_DIR / "demographics.csv",
+            output_dir=out,
+            group_col="group",
+            subject_id="subject_id",
+            n_perms=200,
+            n_bootstraps=200,
+            seed=42,
+            img_format="png",
+            dpi=72,
+            bsr_threshold=bsr_threshold,
+        )
+        return out / "data" / "subject_scores.csv"
+
+    def test_high_threshold_drops_all_lvs(self, tmp_path):
+        # At the default threshold a latent variable survives and scores are
+        # written; an unreachably high threshold makes no feature reliable, so
+        # no LV survives and the scores file is not written.
+        assert self._run(tmp_path / "default", 1.96).exists()
+        assert not self._run(tmp_path / "high", 1e6).exists()

From a66e4189551ebd6c5b6726379459735848535cde Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 21:19:11 +0100
Subject: [PATCH 08/15] enh: record dependency versions and fix set_ticklabels
 warning

- _write_log now records numpy/scipy/scikit-learn versions for
  reproducibility (sklearn reported as 'not installed' when absent).
- filter_lvs reuses self.significant_lvs instead of recomputing
  p_values < 0.05 (behaviour identical).
- pin the tick locator before relabelling in plot_scores_scatter so
  runs are free of the matplotlib set_ticklabels warning.
---
 plsdo/core.py          |  2 +-
 plsdo/pipeline.py      | 22 ++++++++++++++++++++++
 plsdo/plotting.py      |  3 +++
 tests/test_pipeline.py | 10 ++++++++++
 4 files changed, 36 insertions(+), 1 deletion(-)

diff --git a/plsdo/core.py b/plsdo/core.py
index 91ad9ef..480d1c3 100644
--- a/plsdo/core.py
+++ b/plsdo/core.py
@@ -186,7 +186,7 @@ def filter_lvs(self, bsr_threshold: float = 1.96) -> None:
         if not self._bootstrapped:
             raise RuntimeError("Call .bootstrap() before .filter_lvs().")
 
-        significant = self.p_values < 0.05
+        significant = self.significant_lvs
 
         # Check if any feature exceeds threshold on X side
         x_reliable = np.any(np.abs(self.u_bootstrap_ratios) > bsr_threshold, axis=0)
diff --git a/plsdo/pipeline.py b/plsdo/pipeline.py
index d75f85e..618b87b 100644
--- a/plsdo/pipeline.py
+++ b/plsdo/pipeline.py
@@ -39,6 +39,25 @@
 logger = logging.getLogger("plsdo")
 
 
+def _dependency_versions() -> dict[str, str]:
+    """Versions of the numerical dependencies, for reproducibility.
+
+    scikit-learn is optional (only the cross-validate path needs it), so it is
+    reported as "not installed" when absent rather than failing.
+    """
+    import numpy
+    import scipy
+
+    versions = {"numpy": numpy.__version__, "scipy": scipy.__version__}
+    try:
+        import sklearn
+
+        versions["scikit-learn"] = sklearn.__version__
+    except ImportError:
+        versions["scikit-learn"] = "not installed"
+    return versions
+
+
 def _write_log(output_dir: Path, params: dict, notes: list[str] | None = None) -> None:
     """Write a log.txt with run parameters and optional trailing notes."""
     log_path = output_dir / "log.txt"
@@ -46,6 +65,9 @@ def _write_log(output_dir: Path, params: dict, notes: list[str] | None = None) -
         f.write("PLS analysis log\n")
         f.write(f"Version: {__version__}\n")
         f.write(f"Timestamp: {datetime.now().isoformat()}\n")
+        f.write("\nLibrary versions:\n")
+        for name, version in _dependency_versions().items():
+            f.write(f"  {name}: {version}\n")
         f.write("\nParameters:\n")
         for k, v in params.items():
             f.write(f"  {k}: {v}\n")
diff --git a/plsdo/plotting.py b/plsdo/plotting.py
index 0579fda..908a6e6 100644
--- a/plsdo/plotting.py
+++ b/plsdo/plotting.py
@@ -366,6 +366,9 @@ def _box_strip_facet(
     )
     if rotate_xticklabels:
         for ax in g.axes.flat:
+            # Pin the tick locator before relabelling: set_xticklabels alone
+            # warns when the number of ticks is not fixed first.
+            ax.set_xticks(ax.get_xticks())
             ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha="right")
     g.add_legend()
     _finalise(g, out_path, dpi)
diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py
index 8cbf6ff..54f51f8 100644
--- a/tests/test_pipeline.py
+++ b/tests/test_pipeline.py
@@ -244,6 +244,16 @@ def test_log_written_with_version_and_params(self, run_out):
         assert __version__ in log
         assert f"method: {method}" in log
 
+    def test_log_records_dependency_versions(self, run_out):
+        import numpy
+        import scipy
+
+        _method, out = run_out
+        log = (out / "log.txt").read_text()
+        assert f"numpy: {numpy.__version__}" in log
+        assert f"scipy: {scipy.__version__}" in log
+        assert "scikit-learn:" in log
+
     def test_core_figures_produced(self, run_out):
         _method, out = run_out
         figs = out / "figures"

From 8864f3aa2fcce8c59a2a40902b86a8fd0bd9e4b9 Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 21:22:21 +0100
Subject: [PATCH 09/15] test: thin redundant tests (leanness)

- delete test_singular_values_descending (a property of numpy's svd)
- delete the duplicate verbose-limit guard test (default-limit test
  already covers the firing path)
- reduce the three subcommand-alias tests to fast stubbed-dispatch
  checks of method/func instead of full pipeline runs

Kept the argparse interface and io/align guards: cheap, deliberate
interface checks not worth the coverage risk to remove.
---
 tests/test_cli.py      | 59 +++++++++++++++++++++---------------------
 tests/test_core.py     | 11 --------
 tests/test_pipeline.py | 31 ----------------------
 3 files changed, 29 insertions(+), 72 deletions(-)

diff --git a/tests/test_cli.py b/tests/test_cli.py
index 321b350..c976a9e 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -68,7 +68,14 @@ def test_discriminatory_without_group_col_errors(self, data_dir, tmp_path, capsy
         captured = capsys.readouterr()
         assert "requires --group-col" in captured.err.lower()
 
-    def test_corr_alias_runs(self, data_dir, tmp_path):
+    def test_corr_alias_dispatches_correlational(self, data_dir, tmp_path, monkeypatch):
+        # The alias only needs to resolve to the correlational subcommand;
+        # the full run is exercised elsewhere. Stub the pipeline to keep this a
+        # fast dispatch check.
+        captured = {}
+        monkeypatch.setattr(
+            "plsdo.pipeline.run_pipeline", lambda **kw: captured.update(kw)
+        )
         pls_main(
             [
                 "corr",
@@ -80,18 +87,17 @@ def test_corr_alias_runs(self, data_dir, tmp_path):
                 str(data_dir / "demographics.csv"),
                 "--output",
                 str(tmp_path / "out"),
-                "--n-perms",
-                "10",
-                "--n-bootstraps",
-                "10",
-                "--subject-id",
-                "subject_id",
             ]
         )
-        assert (tmp_path / "out" / "data").exists()
+        assert captured["method"] == "correlational"
 
-    def test_discrim_alias_runs(self, data_dir, tmp_path):
-        out = tmp_path / "out_discrim"
+    def test_discrim_alias_dispatches_discriminatory(
+        self, data_dir, tmp_path, monkeypatch
+    ):
+        captured = {}
+        monkeypatch.setattr(
+            "plsdo.pipeline.run_pipeline", lambda **kw: captured.update(kw)
+        )
         pls_main(
             [
                 "discrim",
@@ -101,17 +107,11 @@ def test_discrim_alias_runs(self, data_dir, tmp_path):
                 str(data_dir / "demographics.csv"),
                 "--group-col",
                 "group",
-                "--subject-id",
-                "subject_id",
                 "--output",
-                str(out),
-                "--n-perms",
-                "10",
-                "--n-bootstraps",
-                "10",
+                str(tmp_path / "out_discrim"),
             ]
         )
-        assert (out / "data").exists()
+        assert captured["method"] == "discriminatory"
 
     def test_group_col_and_groups_mutually_exclusive(
         self,
@@ -312,8 +312,14 @@ def test_runs_successfully(self, data_dir, tmp_path):
         assert (out / "data").exists()
         assert (out / "log.txt").exists()
 
-    def test_cv_alias_runs(self, data_dir, tmp_path):
-        out = tmp_path / "cv_alias"
+    def test_cv_alias_dispatches_cross_validate(self, data_dir, tmp_path, monkeypatch):
+        # The alias only needs to resolve to the cross-validate subcommand; the
+        # full run is exercised by test_runs_successfully. Stub the pipeline.
+        called = {}
+        monkeypatch.setattr(
+            "plsdo.pipeline.cross_validate_pipeline",
+            lambda **kw: called.update(kw, dispatched=True),
+        )
         pls_main(
             [
                 "cv",
@@ -323,19 +329,12 @@ def test_cv_alias_runs(self, data_dir, tmp_path):
                 str(data_dir / "demographics.csv"),
                 "--group-col",
                 "group",
-                "--subject-id",
-                "subject_id",
                 "--output",
-                str(out),
-                "--n-folds",
-                "3",
-                "--n-repeats",
-                "2",
-                "--n-permutations",
-                "10",
+                str(tmp_path / "cv_alias"),
             ]
         )
-        assert (out / "data").exists()
+        assert called["dispatched"] is True
+        assert called["group_col"] == "group"
 
     def test_accepts_groups_yaml(self, data_dir, tmp_path):
         out = tmp_path / "cv_yaml"
diff --git a/tests/test_core.py b/tests/test_core.py
index 27a66ad..7176f92 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -74,17 +74,6 @@ def test_sign_convention_largest_loading_positive(self, x_array, y_array):
             col = model.u_loadings[:, i]
             assert col[np.argmax(np.abs(col))] > 0
 
-    def test_singular_values_descending(self, x_array, y_array):
-        from plsdo.io import zscore_columns
-
-        X = zscore_columns(x_array)
-        Y = zscore_columns(y_array)
-        model = PLS(X, Y)
-        model.fit()
-
-        assert all(model.s[i] >= model.s[i + 1] for i in range(len(model.s) - 1))
-
-
 class TestPermutationTest:
     def _fitted_model(self, x_array, y_array):
         from plsdo.io import zscore_columns
diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py
index 54f51f8..69f30a9 100644
--- a/tests/test_pipeline.py
+++ b/tests/test_pipeline.py
@@ -111,37 +111,6 @@ def test_explicit_override_bypasses_guard(self, figures_dir, caplog):
         assert len(produced) > 1
         assert "scree.svg" in produced
 
-    def test_explicit_lower_limit_fires_guard(self, figures_dir, caplog):
-        """60 features with limit=50 — guard fires."""
-        n_features = 60
-        model = _make_mock_model(n_x=n_features, n_y=n_features)
-
-        with caplog.at_level(logging.WARNING, logger="plsdo"):
-            _plot_verbose(
-                model=model,
-                method="correlational",
-                X=np.zeros((10, n_features)),
-                Y=np.zeros((10, n_features)),
-                x_feature_names=[f"x{i}" for i in range(n_features)],
-                x_display_names=[f"x{i}" for i in range(n_features)],
-                y_feature_names=[f"y{i}" for i in range(n_features)],
-                x_colours=None,
-                y_colours=None,
-                final_lv_indices=np.array([0]),
-                final_lv_names=["LV1"],
-                config=None,
-                demo_aligned=None,
-                figures_dir=figures_dir,
-                ext="svg",
-                dpi=72,
-                verbose_feature_limit=50,
-            )
-
-        assert any("Skipping verbose plots" in msg for msg in caplog.messages)
-        produced = sorted(p.name for p in figures_dir.iterdir())
-        assert produced == ["scree.svg"]
-
-
 class TestMultiIndexSubjectScores:
     """Integration: compound subject ID produces a two-level index in CSV."""
 

From ffcc53bc30c13c26e89869c0792f41cb6e8d6c29 Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 21:24:01 +0100
Subject: [PATCH 10/15] chore: silence seaborn's matplotlib bxp deprecation
 warning

seaborn 0.13.2 (the latest release) calls bxp(vert=...), which
matplotlib >= 3.10 deprecates; the fix is only on seaborn's unreleased
main branch, so there is no version to bump to. Ignore this specific
third-party PendingDeprecationWarning until a seaborn release lands.
Test runs are now warning-clean.
---
 pyproject.toml | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/pyproject.toml b/pyproject.toml
index fd94d99..cb78099 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -61,6 +61,12 @@ path = "plsdo/__init__.py"
 testpaths = ["tests"]
 filterwarnings = [
     "error::UserWarning:seaborn",
+    # seaborn 0.13.2 (the latest release) calls matplotlib's bxp() with the
+    # vert= argument, which matplotlib >= 3.10 deprecates in favour of
+    # orientation=. The fix exists only on seaborn's unreleased main branch, so
+    # there is no version to bump to; ignore this third-party warning until a
+    # seaborn release lands, then remove this line.
+    "ignore:vert.*deprecated:PendingDeprecationWarning",
 ]
 
 [tool.coverage.run]

From c2293a4cf756cc647c41fa10220d32862a87e8cc Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 21:25:14 +0100
Subject: [PATCH 11/15] test: add SVD-path and CV edge-case coverage

- single X feature (n_x=1) runs end-to-end with one component
- many-group discriminatory design gives min(n_groups, n_y) components
- CV with fewer subjects than folds raises a clear ValueError
---
 tests/test_core.py           | 43 ++++++++++++++++++++++++++++++++++++
 tests/test_cross_validate.py |  9 ++++++++
 2 files changed, 52 insertions(+)

diff --git a/tests/test_core.py b/tests/test_core.py
index 7176f92..ba30898 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -434,6 +434,49 @@ def test_procrustes_keeps_degenerate_loadings_reliable(self):
         assert np.all(np.abs(model.u_bootstrap_ratios[:4, 0]) > 1.96)
 
 
+class TestEdgeCases:
+    """SVD-path edge cases: extreme feature counts and the discriminatory
+    many-group design. Assert shapes hold so the engine degrades gracefully."""
+
+    def test_single_x_feature(self):
+        """A single X feature yields one component and runs end-to-end."""
+        from plsdo.io import zscore_columns
+
+        rng = np.random.default_rng(0)
+        X = zscore_columns(rng.standard_normal((20, 1)))
+        Y = zscore_columns(rng.standard_normal((20, 4)))
+        model = PLS(X, Y, seed=1)
+        model.fit()
+        model.permutation_test(n_perms=50)
+        model.bootstrap(n_bootstraps=50)
+        model.filter_lvs()
+
+        assert model.s.shape == (1,)
+        assert model.u.shape == (1, 1)
+        assert model.vt.shape == (1, 4)
+        assert model.u_bootstrap_ratios.shape == (1, 1)
+        assert model.final_lvs.shape == (1,)
+
+    def test_many_group_discriminatory(self):
+        """Dummy-coded X with five groups gives min(n_groups, n_y) components."""
+        from plsdo.io import zscore_columns
+
+        labels = np.repeat(np.arange(5), 5)  # 5 groups, 25 subjects
+        X = np.eye(5)[labels].astype(float)  # discriminatory design (not z-scored)
+        Y = zscore_columns(np.random.default_rng(0).standard_normal((25, 3)))
+        model = PLS(X, Y, seed=1, zscore_x=False)
+        model.fit()
+        model.permutation_test(n_perms=30)
+        model.bootstrap(n_bootstraps=30)
+        model.filter_lvs()
+
+        n_components = min(5, 3)
+        assert model.s.shape == (n_components,)
+        assert model.u.shape == (5, n_components)
+        assert model.vt.shape == (n_components, 3)
+        assert model.final_lvs.shape == (n_components,)
+
+
 class TestInvariants:
     """Cheap properties any correct PLS decomposition must satisfy. No
     reference values needed; these guard against a future ``_decompose`` (e.g.
diff --git a/tests/test_cross_validate.py b/tests/test_cross_validate.py
index 269fb44..bf4854a 100644
--- a/tests/test_cross_validate.py
+++ b/tests/test_cross_validate.py
@@ -119,6 +119,15 @@ def test_predicted_labels_within_group_range(self):
         assert set(results["pred_labels"]).issubset({0, 1, 2})
 
 
+class TestCVEdgeCases:
+    def test_too_few_subjects_raises_clear_error(self):
+        """Fewer subjects than folds must fail loudly, not silently."""
+        X = np.random.default_rng(0).standard_normal((4, 5))
+        labels = np.array([0, 0, 1, 1])
+        with pytest.raises(ValueError, match="n_splits"):
+            run_cv(X, labels, n_splits=5, n_repeats=1, n_components=1, seed=0)
+
+
 class TestSklearnImportGuard:
     def test_missing_sklearn_raises_helpful_error(self, monkeypatch):
         """Importing cross_validate without scikit-learn points at plsdo[cv]."""

From 76cf96e7b030efafa56720a4db5718e218041fc6 Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 21:26:45 +0100
Subject: [PATCH 12/15] test: lock the deliberate CV X/Y-flip

Comment the flip at both decision points and add a test asserting
cross_validate_pipeline passes the continuous Y-matrix as run_cv's
predictor and the group codes as the target (opposite to the
discriminatory run_pipeline convention).
---
 plsdo/cross_validate.py      |  4 ++++
 plsdo/pipeline.py            |  3 +++
 tests/test_cross_validate.py | 41 ++++++++++++++++++++++++++++++++++++
 3 files changed, 48 insertions(+)

diff --git a/plsdo/cross_validate.py b/plsdo/cross_validate.py
index df394c4..8d6e53c 100644
--- a/plsdo/cross_validate.py
+++ b/plsdo/cross_validate.py
@@ -51,6 +51,10 @@ def run_cv(
         mean_accuracy, mean_balanced_accuracy, fold_results (DataFrame),
         true_labels, pred_labels, confusion_matrix
     """
+    # Deliberate X/Y flip, opposite to `plsdo discriminatory`: the continuous
+    # data ``X`` is the *predictor* and the dummy-coded groups are the
+    # *target*, so PLSRegression.predict yields predicted group scores that
+    # argmax into class labels.
     n_groups = len(np.unique(labels))
     Y_dummy = np.eye(n_groups)[labels]
 
diff --git a/plsdo/pipeline.py b/plsdo/pipeline.py
index 618b87b..aec07cd 100644
--- a/plsdo/pipeline.py
+++ b/plsdo/pipeline.py
@@ -569,6 +569,9 @@ def cross_validate_pipeline(
 
     # --- Run CV ---
     logger.info("Running %d-fold CV with %d repeats...", n_folds, n_repeats)
+    # X/Y flip: the continuous Y-matrix is the CV *predictor* and the
+    # demographic groups are the classification *target* (opposite to the
+    # discriminatory run_pipeline convention).
     cv_result = run_cv(
         Y,
         labels,
diff --git a/tests/test_cross_validate.py b/tests/test_cross_validate.py
index bf4854a..29da171 100644
--- a/tests/test_cross_validate.py
+++ b/tests/test_cross_validate.py
@@ -1,3 +1,5 @@
+from pathlib import Path
+
 import numpy as np
 import pytest
 from plsdo.cross_validate import run_cv, permutation_test_cv
@@ -119,6 +121,45 @@ def test_predicted_labels_within_group_range(self):
         assert set(results["pred_labels"]).issubset({0, 1, 2})
 
 
+class TestCVFlip:
+    def test_pipeline_passes_continuous_data_as_predictor(self, monkeypatch, tmp_path):
+        """Lock the deliberate X/Y flip: cross_validate_pipeline must pass the
+        continuous Y-matrix as run_cv's predictor and the group codes as the
+        target — the opposite of the discriminatory run_pipeline convention."""
+        import plsdo.cross_validate as cv_mod
+        from plsdo.pipeline import cross_validate_pipeline
+
+        data_dir = Path(__file__).parent / "data"
+        real_run_cv = cv_mod.run_cv
+        calls = []
+
+        def spy(X, labels, **kwargs):
+            calls.append((np.asarray(X), np.asarray(labels)))
+            return real_run_cv(X, labels, **kwargs)
+
+        monkeypatch.setattr("plsdo.cross_validate.run_cv", spy)
+
+        cross_validate_pipeline(
+            y_path=data_dir / "behaviour.csv",
+            demographics_path=data_dir / "demographics.csv",
+            output_dir=tmp_path / "out",
+            group_col="group",
+            subject_id="subject_id",
+            n_folds=3,
+            n_repeats=2,
+            n_permutations=5,
+            seed=42,
+            img_format="png",
+            dpi=72,
+        )
+
+        predictor, target = calls[0]
+        # behaviour.csv has four continuous features → predictor is the data.
+        assert predictor.shape[1] == 4
+        # The target is the integer group codes (3 groups: A, B, C).
+        assert set(np.unique(target)) == {0, 1, 2}
+
+
 class TestCVEdgeCases:
     def test_too_few_subjects_raises_clear_error(self):
         """Fewer subjects than folds must fail loudly, not silently."""

From e2123dcf6c995568a859541e16b15002d60064d5 Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 21:28:45 +0100
Subject: [PATCH 13/15] test: guard that the additive multifactor trailing LV
 is inert
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

An additive K-factor dummy design is rank-deficient by K-1. Assert the
degenerate trailing LV has ~zero singular value, is non-significant,
is dropped by filter_lvs, and has ~zero loadings — verifying the design
is harmless without resorting to contrast coding.
---
 tests/test_core.py | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/tests/test_core.py b/tests/test_core.py
index ba30898..b931944 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -457,6 +457,29 @@ def test_single_x_feature(self):
         assert model.u_bootstrap_ratios.shape == (1, 1)
         assert model.final_lvs.shape == (1,)
 
+    def test_additive_multifactor_trailing_lv_is_inert(self):
+        """An additive K-factor dummy design is rank-deficient by K-1. The
+        degenerate trailing latent variable must be harmless: ~zero singular
+        value, non-significant, dropped by filter_lvs, with ~zero loadings.
+        (Guards the design without switching to contrast coding.)"""
+        from plsdo.io import zscore_columns
+
+        a = np.repeat(np.arange(3), 4)  # factor A, 3 levels, 12 subjects
+        b = np.tile([0, 1], 6)  # factor B, 2 levels
+        X = np.column_stack([np.eye(3)[a], np.eye(2)[b]]).astype(float)
+        Y = zscore_columns(np.random.default_rng(0).standard_normal((12, 4)))
+
+        model = PLS(X, Y, seed=42, zscore_x=False)
+        model.fit()
+        model.permutation_test(n_perms=100)
+        model.bootstrap(n_bootstraps=100)
+        model.filter_lvs()
+
+        assert model.s[-1] < 1e-8
+        assert not model.significant_lvs[-1]
+        assert not model.final_lvs[-1]
+        assert np.abs(model.u_loadings[:, -1]).max() < 1e-6
+
     def test_many_group_discriminatory(self):
         """Dummy-coded X with five groups gives min(n_groups, n_y) components."""
         from plsdo.io import zscore_columns

From 7ff7fa6984b9b79e2d66f7f72774339766d49a90 Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Thu, 4 Jun 2026 21:31:10 +0100
Subject: [PATCH 14/15] docs: clarify the seaborn warning-ignore is an
 indefinite wait

The bxp vert= fix is merged on seaborn main (PR #3820) but unreleased;
seaborn has not shipped a release since 0.13.2 (Jan 2024), so note
there is no firm timeline rather than implying an imminent bump.
---
 pyproject.toml | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index cb78099..4ef17f2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -61,11 +61,12 @@ path = "plsdo/__init__.py"
 testpaths = ["tests"]
 filterwarnings = [
     "error::UserWarning:seaborn",
-    # seaborn 0.13.2 (the latest release) calls matplotlib's bxp() with the
-    # vert= argument, which matplotlib >= 3.10 deprecates in favour of
-    # orientation=. The fix exists only on seaborn's unreleased main branch, so
-    # there is no version to bump to; ignore this third-party warning until a
-    # seaborn release lands, then remove this line.
+    # seaborn 0.13.2 (the latest release, Jan 2024) calls matplotlib's bxp()
+    # with vert=, which matplotlib >= 3.10 deprecates in favour of orientation=.
+    # The fix is merged on seaborn's main branch (PR #3820) but unreleased, and
+    # seaborn has not cut a release in over two years — so there is no version
+    # to bump to and no firm timeline. Ignore this third-party warning; remove
+    # this line if/when seaborn next releases.
     "ignore:vert.*deprecated:PendingDeprecationWarning",
 ]
 

From 98f5c539d8844b2cd762d558728299891dfc96f0 Mon Sep 17 00:00:00 2001
From: Eilidh MacNicol <eilidhmacnicol@gmail.com>
Date: Fri, 5 Jun 2026 12:52:12 +0100
Subject: [PATCH 15/15] test: pin snapshot p-values and bootstrap ratios only
 for significant LVs

Non-significant latent variables are discarded by filter_lvs and their permutation p-values and bootstrap ratios are not reproducible across BLAS builds: a near-zero or weak singular value puts the observed statistic in a dense pile of near-equal null values (flipping permutation counts), and a near-zero loading makes the ratio noise/noise. This drifted up to 0.27 on the Linux CI runners while passing on macOS. Restrict these two snapshots to the significant LVs; the deterministic snapshots (singular values, loadings, scores) still cover all LVs and catch the result-defining drift.
---
 tests/test_regression.py | 68 ++++++++++++++++++++++++++++++++--------
 1 file changed, 55 insertions(+), 13 deletions(-)

diff --git a/tests/test_regression.py b/tests/test_regression.py
index 4cfcc87..235ab53 100644
--- a/tests/test_regression.py
+++ b/tests/test_regression.py
@@ -9,12 +9,16 @@
 Tolerances are tiered by how reproducible each quantity is across machines:
 
 * deterministic outputs (singular values, loadings, subject scores) are pinned
-  tightly — they depend only on the SVD and the fixed sign convention;
-* permutation p-values use an absolute tolerance, since they are counts and a
-  near-tied null singular value can shift the count by one across BLAS builds;
-* bootstrap ratios use a looser tolerance, as they accumulate per-resample SVD
-  differences. A genuine code regression moves these far more than the
-  tolerance, so the guard still bites.
+  tightly over all latent variables — they depend only on the SVD and the fixed
+  sign convention, and are reproducible across BLAS builds;
+* permutation p-values and bootstrap ratios are pinned only for the *significant*
+  latent variables (see ``_significant_lvs``); for non-significant LVs these
+  quantities are both meaningless (discarded by ``filter_lvs``) and unstable
+  across BLAS builds, so pinning them gives false failures, not drift detection.
+  p-values use an absolute tolerance (they are counts); bootstrap ratios use a
+  looser tolerance, as they accumulate per-resample SVD differences. A genuine
+  code regression moves these far more than the tolerance, so the guard still
+  bites.
 
 Regenerating the snapshot is a deliberate, reviewed step (e.g. a major
 numpy/scipy bump): delete ``tests/data/regression/`` and run the suite once to
@@ -42,12 +46,16 @@ def _load(path: Path, index_col=None) -> np.ndarray:
     return pd.read_csv(path, index_col=index_col).to_numpy(dtype=float)
 
 
-def _snapshot(name: str, arr: np.ndarray, tol: dict) -> None:
+def _snapshot(name: str, arr: np.ndarray, tol: dict, cols=None) -> None:
     """Compare ``arr`` against the committed snapshot, or create it if absent.
 
     On a fresh ``tests/data/regression/`` (deliberate regeneration) the
     snapshot is written and the assertion is skipped for that quantity; on all
-    later runs the committed snapshot is asserted against.
+    later runs the committed snapshot is asserted against. The full array is
+    always saved; ``cols`` restricts only the *comparison* to a subset of
+    latent-variable columns (the last axis), so the snapshot file stays
+    complete and transparent while the assertion ignores quantities that are
+    not reproducible across BLAS builds (see ``_significant_lvs``).
     """
     ref_path = REF_DIR / f"{name}.npy"
     arr = np.asarray(arr)
@@ -55,9 +63,33 @@ def _snapshot(name: str, arr: np.ndarray, tol: dict) -> None:
         REF_DIR.mkdir(parents=True, exist_ok=True)
         np.save(ref_path, arr)
         return
-    np.testing.assert_allclose(
-        arr, np.load(ref_path), err_msg=f"regression drift in '{name}'", **tol
-    )
+    ref = np.load(ref_path)
+    if cols is not None:
+        arr, ref = arr[..., cols], ref[..., cols]
+    np.testing.assert_allclose(arr, ref, err_msg=f"regression drift in '{name}'", **tol)
+
+
+def _significant_lvs(method: str):
+    """Column mask of the latent variables significant in the reference snapshot.
+
+    Permutation p-values and bootstrap ratios are pinned only for these LVs.
+    A non-significant LV is discarded by ``filter_lvs`` and never affects a
+    result, and its p-value and bootstrap ratios are also unstable across BLAS
+    builds: a near-zero or weak singular value puts the observed statistic in a
+    dense pile of near-equal null values, so sub-precision SVD differences flip
+    many permutation counts, and a near-zero loading makes the bootstrap ratio a
+    ratio of numerical noise. Non-significant LVs remain guarded by the
+    deterministic singular-value, loading and score snapshots, which are
+    reproducible across machines and catch the result-defining drift (the
+    1/(n-1) divisor, the z-scoring ddof, the sign convention).
+
+    Returns ``None`` during a deliberate regeneration (no committed p-values
+    yet), which leaves the comparison unrestricted.
+    """
+    p_path = REF_DIR / f"{method}_p_values.npy"
+    if not p_path.exists():
+        return None
+    return np.load(p_path).ravel() < 0.05
 
 
 def _run_pls(method: str, out: Path) -> Path:
@@ -113,7 +145,9 @@ class TestPipelineSnapshot:
     @pytest.mark.parametrize("method", ["correlational", "discriminatory"])
     def test_deterministic_outputs(self, method, request):
         data = request.getfixturevalue(f"{method}_data")
-        _snapshot(f"{method}_singular_values", _load(data / "singular_values.csv"), TIGHT)
+        _snapshot(
+            f"{method}_singular_values", _load(data / "singular_values.csv"), TIGHT
+        )
         _snapshot(
             f"{method}_x_loadings", _load(data / "x_loadings.csv", index_col=0), TIGHT
         )
@@ -129,20 +163,28 @@ def test_deterministic_outputs(self, method, request):
     @pytest.mark.parametrize("method", ["correlational", "discriminatory"])
     def test_permutation_pvalues(self, method, request):
         data = request.getfixturevalue(f"{method}_data")
-        _snapshot(f"{method}_p_values", _load(data / "p_values.csv"), PVAL)
+        _snapshot(
+            f"{method}_p_values",
+            _load(data / "p_values.csv"),
+            PVAL,
+            cols=_significant_lvs(method),
+        )
 
     @pytest.mark.parametrize("method", ["correlational", "discriminatory"])
     def test_bootstrap_ratios(self, method, request):
         data = request.getfixturevalue(f"{method}_data")
+        sig = _significant_lvs(method)
         _snapshot(
             f"{method}_x_bootstrap_ratios",
             _load(data / "x_bootstrap_ratios.csv", index_col=0),
             LOOSE,
+            cols=sig,
         )
         _snapshot(
             f"{method}_y_bootstrap_ratios",
             _load(data / "y_bootstrap_ratios.csv", index_col=0),
             LOOSE,
+            cols=sig,
         )
 
     def test_cross_validation(self, cv_data):