diff --git a/docs/changelog.rst b/docs/changelog.rst index 6ee6938..d41cea0 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -36,8 +36,8 @@ Added: ``||X||**2 - ``. 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 diff --git a/docs/validation.rst b/docs/validation.rst index 0513dd5..7c6e6b2 100644 --- a/docs/validation.rst +++ b/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 diff --git a/scripts/make_fit_references.py b/scripts/make_fit_references.py deleted file mode 100644 index 8bfdf0c..0000000 --- a/scripts/make_fit_references.py +++ /dev/null @@ -1,35 +0,0 @@ -"""Write tests/data/fits_3298278.npz: the fits of tests/test_fit_equivalence.py's -FIT_CASES by the code that is importable as ``htfa`` (run it with the code of -commit 3298278 first on the path). - - git archive 3298278 htfa | tar -x -C /tmp/htfa_3298278 - PYTHONPATH=/tmp/htfa_3298278 python scripts/make_fit_references.py 3298278 -""" - -import datetime -import os -import platform -import sys - -import numpy as np - -here = os.path.dirname(os.path.abspath(__file__)) -sys.path.insert(0, os.path.join(here, "..", "tests")) -import test_fit_equivalence as fe # noqa: E402 - -commit = sys.argv[1] -arrays = {} -for i, case in enumerate(fe.FIT_CASES): - for k, v in fe.fit_case(case).items(): - if k in ("centers", "widths", "weights", "cost"): - arrays[f"case{i}_{k}"] = v -arrays["provenance"] = np.array( - f"commit {commit}; TFA defaults, FIT_CASES of tests/test_fit_equivalence.py; numpy " - f"{np.__version__}; {platform.platform()}; " - f"{datetime.datetime.now(datetime.timezone.utc).isoformat(timespec='seconds')}" -) -os.makedirs(os.path.join(here, "..", "tests", "data"), exist_ok=True) -np.savez_compressed( - os.path.join(here, "..", "tests", "data", f"fits_{commit}.npz"), **arrays -) -print("wrote", len(fe.FIT_CASES), "fits;", arrays["provenance"]) diff --git a/tests/data/fits_3298278.npz b/tests/data/fits_3298278.npz deleted file mode 100644 index df8c898..0000000 Binary files a/tests/data/fits_3298278.npz and /dev/null differ diff --git a/tests/test_fit_equivalence.py b/tests/test_fit_equivalence.py index 8b2f4d2..090147c 100644 --- a/tests/test_fit_equivalence.py +++ b/tests/test_fit_equivalence.py @@ -1,31 +1,40 @@ -"""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) @@ -33,16 +42,39 @@ (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}" diff --git a/tests/test_tfa.py b/tests/test_tfa.py index 96a3c85..0c7b82f 100644 --- a/tests/test_tfa.py +++ b/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_)