Task 17: platform-independent whole-fit check (CI fix)

docs/changelog.rst

@@ -36,8 +36,8 @@ Added:
   ``||X||**2 - <F^T X, W>``. A cost evaluation is about 3 times and an
   accepted solver step 1.4 times faster on development_fmri (K = 50 and 100);
   whole fits 1.1 to 1.5 times (see :doc:`performance`).
-  Checked against the earlier formulas over 492 generated scenarios and 36
-  saved synthetic fits (see :doc:`validation`).
+  Checked against the earlier formulas over 492 generated scenarios and 35
+  well-conditioned synthetic fits (see :doc:`validation`).
 - An ``'ols'`` fit whose width shrank (with a small ``lower_ratio``) until a
   factor was 0 at every voxel raised ``LinAlgError`` in its final weights;
   ``ridge_weights`` now falls back to the minimum-norm weights the varpro

docs/validation.rst

@@ -185,16 +185,22 @@ keyed on the widths only failed 32 of its 40 cases, one keyed on the centers
 only 32, and a cost cached on the point alone (ignoring prior changes) 27.
 
 **Whole fits, well-conditioned synthetic problems**
-(``tests/test_fit_equivalence.py``). 36 fits with the defaults (2D grids 20 x
-20, 30 x 30, 25 x 15 and 3D grids 10³, 12³, 14 x 10 x 8; K = 2, 4, 6; noise 0.1
-and 0.5) against the saved fits of ``3298278``
-(``tests/data/fits_3298278.npz``, written by ``scripts/make_fit_references.py``).
-35 reference fits recover the true centers to 0.2 lengths; one ends at a local
-optimum in both versions. Tolerances are those of the stopping rule
-(``threshold`` = 1e-3: centers in lengths ``sqrt(max_sigma / 2)``, widths and
-weights relative) and 1e-8 for the cost.
-
-.. list-table:: Whole-fit differences, 36 fits
+(``tests/test_fit_equivalence.py``). 36 candidate problems (2D grids 20 x 20,
+30 x 30, 25 x 15 and 3D grids 10³, 12³, 14 x 10 x 8; K = 2, 4, 6; noise 0.1
+and 0.5), fitted with the defaults by the current code and by the same
+``TFA`` running on the 3298278 kernels, both on the machine that runs the test
+(saved outputs would carry one platform's BLAS rounding). A problem counts as
+well-conditioned when the 3298278 fit recovers the true centers to 0.2
+lengths ``sqrt(max_sigma / 2)`` and is stable: its fit on the data changed by
+a relative 1e-12 agrees with it within the tolerances below. Only those are
+compared; at least 30 must qualify. Tolerances are those of the stopping rule
+(``threshold`` = 1e-3: centers in lengths, widths and weights relative) and
+1e-8 for the cost. On macOS (Accelerate) 35 qualified (the one excluded
+ends at a local optimum, centers 1.75 lengths from the truth); on Linux
+(OpenBLAS, Python 3.10) all 36, with differences of at most 1.5e-9 (centers),
+6.2e-8 (widths), 2.8e-8 (weights) and 6.3e-11 (cost).
+
+.. list-table:: Whole-fit differences, 35 well-conditioned fits (macOS)
    :header-rows: 1
 
    * - quantity

tests/test_fit_equivalence.py

@@ -1,48 +1,80 @@
-"""Whole TFA fits against saved fits of commit 3298278 on well-conditioned problems.
-
-``tests/data/fits_3298278.npz`` holds the centers, widths, weights and cost of
-`FIT_CASES` fitted by the code of commit 3298278 (written by
-``scripts/make_fit_references.py``; provenance in the file). Each case is a
-synthetic problem whose reference fit recovers the true centers (checked here),
-so its optimum is well determined (one case, `LOCAL_OPTIMA`, ends at a local
-optimum in both versions and is compared all the same).
-
-A fit stops when no center moved by more than ``threshold`` (1e-3) times the
-length ``sqrt(max_sigma / 2)`` and no width by more than ``threshold``
-relatively in the last iteration, so two fits whose solver arithmetic differs by
-rounding can end that far apart. Those are the tolerances; the measured
-differences are printed (and recorded in docs/validation.rst).
+"""Whole TFA fits by the current kernels against the 3298278 kernels, on the same
+machine, on problems that are demonstrably well-conditioned.
+
+The oracle fit is ``htfa.tfa.TFA`` running on ``tests/reference_kernels_3298278.py``
+(`oracle_kernels`): ``htfa/tfa.py`` is unchanged since commit 3298278, so this
+is the 3298278 fit, computed with this machine's BLAS. Both fits therefore see
+the same platform rounding, and only the kernels differ.
+
+A case counts as well-conditioned when the oracle fit (a) recovers the true
+centers to `RECOVERY_LENGTHS` lengths ``sqrt(max_sigma / 2)`` and (b) is stable:
+the oracle fit on the data changed by a relative 1e-12 (`PERTURBATION`) agrees
+with it within the stopping rule's tolerances below. Only such cases are
+compared; the others (local optima, or optima that rounding moves) are counted
+and printed.
+
+A fit stops when no center moved by more than ``threshold`` (1e-3) lengths and
+no width by more than ``threshold`` relatively in the last iteration, so two
+fits whose arithmetic differs by rounding can end that far apart. Those are the
+tolerances for centers, widths and weights; the cost's is 1e-8.
 """
 
-import os
+from contextlib import contextmanager
 
 import numpy as np
-import pytest
+import reference_kernels_3298278 as ref
 
+import htfa.tfa as tfa_module
 from htfa._kernels import max_sigma
 from htfa._matching import match_factors
 from htfa.synthetic import make_tfa_data
 from htfa.tfa import TFA
 
 THRESHOLD = 1e-3  # TFA's default threshold
-REFERENCE = os.path.join(os.path.dirname(__file__), "data", "fits_3298278.npz")
+COST_RTOL = 1e-8
+RECOVERY_LENGTHS = 0.2
+PERTURBATION = 1e-12
+MIN_WELL_CONDITIONED = 30
 GRIDS = ((20, 20), (30, 30), (25, 15), (10, 10, 10), (12, 12, 12), (14, 10, 8))
 FIT_CASES = [
     dict(grid_shape=g, n_samples=40, K=K, noise=noise, random_state=100 + i)
     for i, (g, K, noise) in enumerate(
         (g, K, noise) for g in GRIDS for K in (2, 4, 6) for noise in (0.1, 0.5)
     )
 ]
-
-
-# Reference fits that ended at a local optimum (max center error 1.75
-# lengths); they are compared all the same, but are not well-conditioned.
-LOCAL_OPTIMA = {23}  # (10, 10, 10), K = 6, noise 0.5
-
-
-def fit_case(case: dict) -> dict[str, np.ndarray]:
+# The kernel names htfa.tfa imports.
+KERNEL_NAMES = (
+    "Prior",
+    "ProjectedProblem",
+    "ResidualProblem",
+    "bounds",
+    "factors",
+    "init_centers_widths",
+    "max_sigma",
+    "ridge_weights",
+    "solve_centers_widths",
+)
+
+
+@contextmanager
+def oracle_kernels():
+    """``htfa.tfa`` with the 3298278 kernels in place of the current ones."""
+    saved = {name: getattr(tfa_module, name) for name in KERNEL_NAMES}
+    try:
+        for name in KERNEL_NAMES:
+            setattr(tfa_module, name, getattr(ref, name))
+        yield
+    finally:
+        for name, value in saved.items():
+            setattr(tfa_module, name, value)
+
+
+def fit_case(case: dict, perturbation: float = 0.0) -> dict[str, np.ndarray]:
     """TFA with its defaults on the synthetic problem of ``case``."""
     X, R, truth = make_tfa_data(**case)
+    if perturbation:
+        rng = np.random.default_rng(case["random_state"] + 1)
+        X = X * (1 + perturbation * rng.standard_normal(X.shape))
     m = TFA(K=case["K"], random_state=case["random_state"]).fit(X, R)
     return dict(
         centers=m.centers_,
@@ -54,58 +86,56 @@ def fit_case(case: dict) -> dict[str, np.ndarray]:
     )
 
 
-@pytest.fixture(scope="module")
-def reference():
-    return np.load(REFERENCE)
+def _differences(a: dict, b: dict) -> dict[str, float]:
+    """Differences of fit ``b`` from fit ``a`` (factors matched by center)."""
+    perm = match_factors(a["centers"], b["centers"])
+    C, w, W = b["centers"][perm], b["widths"][perm], b["weights"][perm]
+    return dict(
+        centers=float(np.max(np.abs(C - a["centers"]))) / float(a["length"]),
+        widths=float(np.max(np.abs(w / a["widths"] - 1))),
+        weights=float(np.max(np.abs(W - a["weights"])) / np.max(np.abs(a["weights"]))),
+        cost=float(abs(b["cost"] - a["cost"]) / a["cost"]),
+    )
 
 
-def _differences(case_index: int, ref) -> dict[str, float]:
-    got = fit_case(FIT_CASES[case_index])
-    key = f"case{case_index}_"
-    C_ref, w_ref = ref[key + "centers"], ref[key + "widths"]
-    perm = match_factors(C_ref, got["centers"])
-    C, w, W = got["centers"][perm], got["widths"][perm], got["weights"][perm]
-    W_ref = ref[key + "weights"]
-    length = float(got["length"])
-    return dict(
-        centers=float(np.max(np.abs(C - C_ref))) / length,
-        widths=float(np.max(np.abs(w / w_ref - 1))),
-        weights=float(np.max(np.abs(W - W_ref)) / np.max(np.abs(W_ref))),
-        cost=float(abs(got["cost"] - ref[key + "cost"]) / ref[key + "cost"]),
-        recovery=float(
-            np.max(
-                np.abs(
-                    C_ref[match_factors(got["true_centers"], C_ref)]
-                    - got["true_centers"]
-                )
-            )
-        )
-        / length,
+def _within(d: dict[str, float]) -> bool:
+    return (
+        d["centers"] <= THRESHOLD
+        and d["widths"] <= THRESHOLD
+        and d["weights"] <= THRESHOLD
+        and d["cost"] <= COST_RTOL
     )
 
 
-def test_reference_provenance(reference):
-    """The reference file says where it came from."""
-    meta = str(reference["provenance"])
-    print(meta)
-    assert "3298278" in meta
-    assert len([k for k in reference.files if k.endswith("_centers")]) == len(FIT_CASES)
+def _recovery(fit: dict) -> float:
+    perm = match_factors(fit["true_centers"], fit["centers"])
+    err = np.max(np.abs(fit["centers"][perm] - fit["true_centers"]))
+    return float(err) / float(fit["length"])
 
 
-def test_fits_match_3298278_on_well_conditioned_problems(reference):
-    """36 fits (2D and 3D grids, K = 2, 4, 6, noise 0.1 and 0.5): centers within
+def test_fits_match_the_3298278_kernels_on_well_conditioned_problems():
+    """36 candidate problems (2D and 3D grids, K = 2, 4, 6, noise 0.1 and 0.5):
+    on every well-conditioned one (see the module docstring), centers within
     THRESHOLD lengths, widths and weights within THRESHOLD (relative), cost
-    within 1e-8 (relative) of the 3298278 fits. All reference fits but those in
-    LOCAL_OPTIMA recover the true centers to 0.2 lengths (a well-determined
-    optimum)."""
-    rows = [_differences(i, reference) for i in range(len(FIT_CASES))]
-    for name in ("centers", "widths", "weights", "cost", "recovery"):
-        v = np.array([r[name] for r in rows])
-        print(f"{name:9s} median {np.median(v):.2e} max {v.max():.2e}")
-    for i, r in enumerate(rows):
-        msg = f"case {i}: {FIT_CASES[i]}: {r}"
-        assert (r["recovery"] < 0.2) != (i in LOCAL_OPTIMA), msg
-        assert r["centers"] <= THRESHOLD, msg
-        assert r["widths"] <= THRESHOLD, msg
-        assert r["weights"] <= THRESHOLD, msg
-        assert r["cost"] <= 1e-8, msg
+    within COST_RTOL of the oracle fit. At least MIN_WELL_CONDITIONED of the
+    candidates must qualify."""
+    rows, excluded = [], []
+    for i, case in enumerate(FIT_CASES):
+        with oracle_kernels():
+            old = fit_case(case)
+            old_perturbed = fit_case(case, PERTURBATION)
+        recovery = _recovery(old)
+        stability = _differences(old, old_perturbed)
+        if recovery >= RECOVERY_LENGTHS or not _within(stability):
+            excluded.append((i, case, f"recovery {recovery:.3f}", stability))
+            continue
+        rows.append((i, case, _differences(old, fit_case(case))))
+    for i, case, rec, stab in excluded:
+        print(f"excluded case {i} {case}: {rec}, perturbed oracle {stab}")
+    print(f"{len(rows)} of {len(FIT_CASES)} cases well-conditioned")
+    for name in ("centers", "widths", "weights", "cost"):
+        v = np.array([d[name] for *_, d in rows])
+        print(f"{name:8s} median {np.median(v):.2e} max {v.max():.2e}")
+    assert len(rows) >= MIN_WELL_CONDITIONED
+    for i, case, d in rows:
+        assert _within(d), f"case {i}: {case}: {d}"

tests/test_tfa.py

@@ -816,10 +816,16 @@ def test_ols_fit_with_a_factor_that_vanishes_finishes():
     """With 'ols' and a small lower_ratio a width can shrink until its factor is
     0 at every voxel; the full-data weights then had a singular gram and
     ridge_weights raised LinAlgError (seeds 0, 3, 4, 18 of this setup). Now
-    the fit finishes, the vanished factor gets a weight of about 0 (the minimum-norm
-    weights the varpro objective uses; measured 8.6e-10 of the largest weight)
-    and cost_ is finite."""
-    from htfa._kernels import gram_solve
+    the fit finishes with the minimum-norm weights the varpro objective uses
+    and a finite cost_. The vanished factor's weight is 0 only up to rounding:
+    the eigensolver mixes the zero column's direction into the eigenvector of
+    the smallest kept eigenvalue by about eps times the kept condition number
+    (eigenvector perturbation with that eigenvalue as the gap), and that
+    eigenvector carries weights of size max|W|; summed over the K - 1 kept
+    directions, |w| <= (K - 1) eps cond max|W|. Measured: 8.6e-10 of the
+    largest weight (macOS, Accelerate) and 1.2e-8 (macOS runners with Python
+    3.10 and 3.11), with eps cond = 1.4e-8 here."""
+    from htfa._kernels import GRAM_RCOND, gram_solve
 
     rng = np.random.default_rng(3)
     R, X = rng.uniform(0, 100, (30, 3)), rng.normal(size=(30, 10))
@@ -835,7 +841,13 @@ def test_ols_fit_with_a_factor_that_vanishes_finishes():
         f"vanished factors {np.flatnonzero(vanished)}, widths {m.widths_}, "
         f"cost {m.cost_:.6g}, weight rows max {np.abs(m.weights_).max(1)}"
     )
+    vals = np.linalg.eigvalsh(F.T @ F)
+    kept = vals[vals > GRAM_RCOND * vals[-1]]
+    cond = float(vals[-1] / kept.min())
+    bound = (len(kept) - 1) * np.finfo(np.float64).eps * cond
+    ratio = np.abs(m.weights_[vanished]).max() / np.abs(m.weights_).max()
+    print(f"vanished weight / max weight {ratio:.2e}, bound {bound:.2e}")
     assert vanished.any()
-    assert np.abs(m.weights_[vanished]).max() <= 1e-8 * np.abs(m.weights_).max()
+    assert ratio <= bound
     np.testing.assert_allclose(m.weights_, W, rtol=1e-12, atol=0)
     assert np.isfinite(m.cost_)