From 4a321f774d4540d4db2662979811cca46e1c8509 Mon Sep 17 00:00:00 2001 From: Ben Mather Date: Mon, 3 Aug 2026 13:45:04 +1000 Subject: [PATCH 1/2] Carry a spectrum's provenance on the spectrum itself `_spectrum_key` was a second instance attribute that had to be kept in step with `last_spectrum` at every exit, and it only ever knew the provenance of the single most recent spectrum. Reuse survived only because `asarray` returns an already-float64 array unchanged, so the arrays could be matched by identity -- correct, but three lines of comment defending a coincidence. `_Spectrum` is a tuple subclass carrying its own `provenance`, so the provenance travels with the object instead of beside it. That removes `_spectrum_key`, the `ours` flag and the identity `zip`, and it is strictly more capable: a spectrum held across an intervening call no longer loses its label. It unpacks, indexes and pickles as a plain triple, so `spectrum=`, `last_spectrum` and Global_CPD's `fit_spectrum(*...)` are unchanged, and a triple built anywhere else still has no provenance and is taken as given. `taper` and `process_subgrid` are dropped from what is checked. Holding a callable on the instance kept its captured scope alive, and worse, made `pickle.dumps(grid.optimise)` raise for a closure or a lambda -- so one `optimise(..., process_subgrid=...)` sent `parallelise_routine` to serial for every later call on that instance, warning that it could not send `func`, which was perfectly picklable. Measured before this change: 48 centroids on 4 processors went 0.16 s -> 0.48 s and stayed there. The misuse worth catching is a spectrum from a different window or centroid, which `window`, `xc`, `yc` and `dof_factor` still catch. `parallelise_routine` now rejects `spectrum` outright, alongside `on_error` and `seed`. One spectrum cannot describe a list of centroids: forwarded, it gave every centroid the same answer -- a flat map that looks like a result -- and since the provenance warning fires per call, whether the user was told depended on how many processors were available. Also, `residuals` claimed `np.errstate` is thread safe where `catch_warnings` is not. It only became context-local in NumPy 2.0 and `pyproject.toml` allows `numpy>=1.20`, so on the supported floor it sets global state through `seterr` just the same. The two arguments that do hold on every version -- it does not swallow unrelated warnings, and it is cheaper to enter several hundred times per fit -- are enough on their own. Tidying alongside: the four identical `spectrum :` parameter blocks collapse to one description on `optimise` and three cross-references, the pattern already used for `dof_factor`; the `optimise` example passes its fitted parameters into `profile`, which is worth 5% of the pair and stops the bound names going nowhere; and the new tests use `monkeypatch`, `pytest.warns` and `recwarn` in place of hand-rolled restores and warning capture. CLAUDE.md's test counts were left behind by the previous two commits and are now measured rather than incremented. Suite 121 -> 123 tests. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01PLgEaDtuZw6QXDphSrzcdQ --- CLAUDE.md | 4 +- pycurious/optimise_bouligand.py | 147 +++++++++++++++++--------------- pycurious/parallel.py | 16 ++++ tests/test_optimise.py | 111 ++++++++++++------------ tests/test_routines.py | 24 ++++++ 5 files changed, 173 insertions(+), 129 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 5b560d9..d9d3192 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -15,8 +15,8 @@ returns a bare number without one is a pre-v2 remnant. ## Commands ```bash -pytest # 111 tests, ~20 s -pytest -m "not slow" # 108 tests, ~12 s -- skips the calibration tests that +pytest # 123 tests, ~13 s +pytest -m "not slow" # 120 tests, ~6 s -- skips the calibration tests that # fit a few hundred realisations pytest tests/test_tanaka.py -q ``` diff --git a/pycurious/optimise_bouligand.py b/pycurious/optimise_bouligand.py index ee22a32..f534e91 100644 --- a/pycurious/optimise_bouligand.py +++ b/pycurious/optimise_bouligand.py @@ -107,6 +107,40 @@ # forward model but a sum of two of them. _CPD = "CPD" +# The arguments a reused spectrum is checked against, in the order +# `_Spectrum.provenance` stores them. `taper` and `process_subgrid` are +# deliberately absent: they are callables, and holding one on the instance both +# keeps its captured scope alive and makes `grid.optimise` unpicklable, which +# drops `pycurious.parallel.CurieParallel.parallelise_routine` back to serial +# with a warning that blames the wrong thing. The misuse worth catching is a +# spectrum from a different window or centroid anyway. +_SPECTRUM_ARGS = ("window", "xc", "yc", "dof_factor") + + +class _Spectrum(tuple): + """ + A `(k, Phi, sigma_Phi)` triple that remembers what it was computed from. + + It is an ordinary tuple everywhere it matters -- it unpacks, indexes, zips + and pickles like the one `pycurious.grid.CurieGrid.window_spectrum` + returns -- so `spectrum=` still takes a plain triple and `last_spectrum` + still hands one back. Carrying the provenance on the spectrum rather than + beside it means it survives being passed from routine to routine, and that + a triple built anywhere else simply has none, which is exactly the + "nothing to check, take it as given" case. + """ + + def __new__(cls, arrays, provenance): + spectrum = super(_Spectrum, cls).__new__(cls, arrays) + spectrum.provenance = provenance + return spectrum + + def __getnewargs__(self): + # a tuple subclass with a two-argument `__new__` cannot be unpickled + # without this, and a grid carrying a `last_spectrum` is pickled every + # time `parallelise_routine` sends one to a worker + return (tuple(self), self.provenance) + def _prior_loc_scale(pdf): """ @@ -216,10 +250,8 @@ def __init__(self, grid, xmin, xmax, ymin, ymax, **kwargs): self.bounds = list(zip(lb, ub)) # the spectrum most recently computed or supplied, ready to be handed - # to another routine at the same centroid -- see `_resolve_spectrum`, - # which uses `_spectrum_key` to catch it being handed to a different one + # to another routine at the same centroid -- see `_resolve_spectrum` self.last_spectrum = None - self._spectrum_key = None self.max_processors = kwargs.pop("max_processors", cpu_count()) @@ -361,9 +393,8 @@ def residuals(self, x, kh, Phi, sigma_Phi, prior=None): This is `numpy.errstate` rather than `warnings.catch_warnings` deliberately. Both silence those, but `catch_warnings` swallows *every* warning raised in the block, including real ones from - elsewhere in the library, and it rewrites a global filter on each - of the several hundred evaluations a fit makes, which is not - thread safe. + elsewhere in the library, and it is the more expensive of the two + to enter several hundred times per fit. """ beta, zt, dz, C = x @@ -450,75 +481,54 @@ def _resolve_spectrum( reweight `sigma`, has already applied that to the array the caller is holding; applying it a second time would compound it. - `self.last_spectrum` is set either way, so a caller can hand what - `optimise` just used straight to `profile` without having to - reconstruct it -- and reconstructing it is easy to get wrong, since - `power` must be 2 and any `process_subgrid` must match. - Nothing is reused implicitly: a routine given no `spectrum` always computes one, and the library never reads `last_spectrum` itself. But a supplied spectrum makes `window`, `xc` and `yc` dead arguments, so handing over the one from a *different* window is accepted in silence - and answers a question the caller did not ask -- measured at a 128 km - spectrum passed to a 512 km call, dz came back 222.8 km against the - 21.7 km that window really gives. - - `_spectrum_key` guards the case that can be guarded. When the spectrum - handed back is the one this instance last computed -- which is what the - documented `spectrum=grid.last_spectrum` idiom passes -- the arguments - it was computed from are known, and disagreeing with them is a warning. - A spectrum from anywhere else has no provenance to check, so it is - taken at face value and the key is cleared rather than guessed at. + and answers a question the caller did not ask. `_Spectrum.provenance` + guards the case that can be guarded: a spectrum this library computed + knows the `_SPECTRUM_ARGS` it came from, and disagreeing with them is a + warning. One built anywhere else -- read from an archive, cast to + float32 -- has no provenance, so it is taken at face value rather than + guessed at. """ + provenance = (window, xc, yc, dof_factor) + if spectrum is None: - spectrum = self._spectrum( - window, xc, yc, taper, process_subgrid, dof_factor, **kwargs + spectrum = _Spectrum( + self._spectrum( + window, xc, yc, taper, process_subgrid, dof_factor, **kwargs + ), + provenance, ) - key = (window, xc, yc, taper, process_subgrid, dof_factor) else: + was = getattr(spectrum, "provenance", None) + k, Phi, sigma_Phi = (np.asarray(a, dtype=float) for a in spectrum) if not (k.shape == Phi.shape == sigma_Phi.shape): raise ValueError( "spectrum must be three arrays of the same shape, got " "{}, {} and {}".format(k.shape, Phi.shape, sigma_Phi.shape) ) - spectrum = (k, Phi, sigma_Phi) - - # `asarray` hands back the same object for an array that is already - # float64, so the arrays of `last_spectrum` survive the conversion - # by identity even though the tuple around them does not. - ours = self.last_spectrum is not None and all( - new is old for new, old in zip(spectrum, self.last_spectrum) - ) - key = self._spectrum_key if ours else None - if ours and key is not None: - mismatched = [ - name - for name, was, now in zip( - ("window", "xc", "yc", "taper", "process_subgrid", - "dof_factor"), - key, - (window, xc, yc, taper, process_subgrid, dof_factor), - ) - if was is not now and was != now - ] - if mismatched: - warnings.warn( - "the supplied spectrum was computed with a different " - "{}, and a supplied spectrum is used as given -- " - "{} of this call {} ignored, so the result describes " - "the window the spectrum came from, not the one asked " - "for here.".format( - ", ".join(mismatched), - ", ".join(mismatched), - "is" if len(mismatched) == 1 else "are", - ), - RuntimeWarning, - stacklevel=3, - ) + spectrum = _Spectrum((k, Phi, sigma_Phi), was) + + mismatched = [] if was is None else [ + name + for name, then, now in zip(_SPECTRUM_ARGS, was, provenance) + if then != now + ] + if mismatched: + names = ", ".join(mismatched) + warnings.warn( + "the supplied spectrum was computed with a different {0}; " + "a supplied spectrum is used as given, so the {0} of this " + "call is ignored and the result describes the window the " + "spectrum came from, not the one asked for here.".format(names), + RuntimeWarning, + stacklevel=3, + ) self.last_spectrum = spectrum - self._spectrum_key = key return spectrum def _bound_arrays(self, free=None): @@ -845,7 +855,8 @@ def optimise( >>> beta, zt, dz, C = grid.optimise(2000e3, xc, yc)[:4] >>> _, _, lo, hi = grid.profile( - ... 2000e3, xc, yc, "CPD", spectrum=grid.last_spectrum) + ... 2000e3, xc, yc, "CPD", spectrum=grid.last_spectrum, + ... beta=beta, zt=zt, dz=dz, C=C) which takes 41% off the pair at a 2049-cell window. """ @@ -1047,10 +1058,8 @@ def profile( dof_factor : float, optional see `pycurious.grid.CurieGrid.window_spectrum` spectrum : tuple (k, Phi, sigma_Phi), optional - a spectrum already in hand -- typically `last_spectrum` from - the `optimise` at this same centroid, which saves recomputing - it. See `optimise` for the idiom. `window`, `xc`, `yc`, - `taper`, `process_subgrid` and `dof_factor` are then unused. + a spectrum already in hand, typically `last_spectrum` from the + `optimise` at this same centroid -- see `optimise` kwargs : keyword arguments passed to `radial_spectrum` @@ -1290,9 +1299,8 @@ def metropolis_hastings( also return a dict of `acceptance`, `burnin_acceptance`, `x_scale` spectrum : tuple (k, Phi, sigma_Phi), optional - a spectrum already in hand, e.g. `last_spectrum` from the - `optimise` at this centroid. `window`, `xc`, `yc`, `taper`, - `process_subgrid` and `dof_factor` are then unused. + a spectrum already in hand, typically `last_spectrum` from + the `optimise` at this same centroid -- see `optimise` Returns: beta : ndarray shape (nsim,) @@ -1508,9 +1516,8 @@ def sensitivity( seed : int, optional seed for reproducibility spectrum : tuple (k, Phi, sigma_Phi), optional - a spectrum already in hand, e.g. `last_spectrum` from the - `optimise` at this centroid. `window`, `xc`, `yc`, `taper`, - `process_subgrid` and `dof_factor` are then unused. + a spectrum already in hand, typically `last_spectrum` from + the `optimise` at this same centroid -- see `optimise` Returns: beta : ndarray shape (nsim,) diff --git a/pycurious/parallel.py b/pycurious/parallel.py index 8e7e370..27f3dc9 100644 --- a/pycurious/parallel.py +++ b/pycurious/parallel.py @@ -135,6 +135,10 @@ def parallelise_routine(self, window, xc_list, yc_list, func, *args, **kwargs): were used. `func` must accept a `seed` keyword if this is supplied. + `spectrum` is rejected outright: one spectrum describes one + window at one centroid, so it cannot mean anything for a list + of them. + Returns: out : list of lists (depends on output of `func` - see notes) @@ -183,6 +187,18 @@ def parallelise_routine(self, window, xc_list, yc_list, func, *args, **kwargs): on_error = kwargs.pop("on_error", "raise") if on_error not in ("raise", "ignore"): raise ValueError("on_error must be 'raise' or 'ignore'") + + # forwarding this would hand every centroid the same spectrum, which + # returns the same answer at each of them -- a flat map that looks like + # a result. The per-centroid provenance warning only fires on the + # serial path, so whether the user is told would depend on nprocs. + if "spectrum" in kwargs: + raise ValueError( + "spectrum describes a single window at a single centroid, so " + "it cannot be shared across a list of them. Drop it and let " + "each centroid compute its own." + ) + seed = kwargs.pop("seed", None) if seed is not None and not getattr(func, "wants_seed", False): diff --git a/tests/test_optimise.py b/tests/test_optimise.py index 70aedac..af8ed73 100644 --- a/tests/test_optimise.py +++ b/tests/test_optimise.py @@ -1,9 +1,11 @@ +import warnings + import pytest import pycurious import numpy as np from scipy.optimize import minimize -from conftest import load_magnetic_anomaly +from conftest import load_magnetic_anomaly, synthetic_grid def test_optimisation_smoke(load_magnetic_anomaly): @@ -219,9 +221,7 @@ def test_valid_numbers(load_magnetic_anomaly): def _shared_spectrum_grid(): - from conftest import synthetic_grid - - return synthetic_grid(pycurious.CurieOptimiseBouligand, n=256, dx=2.0) + return synthetic_grid(pycurious.CurieOptimiseBouligand, n=256, dx=2.0)[:3] @pytest.mark.parametrize("target", ["dz", "CPD"]) @@ -233,7 +233,7 @@ def test_supplied_spectrum_reproduces_the_computed_one(target): its own code path is that the two cannot drift apart. Assert to the bit, so that they cannot. """ - grid, xc, yc, extent = _shared_spectrum_grid() + grid, xc, yc = _shared_spectrum_grid() window = 200e3 stock = grid.optimise(window, xc, yc) @@ -246,7 +246,7 @@ def test_supplied_spectrum_reproduces_the_computed_one(target): assert np.array_equal(lhs, rhs) -def test_supplied_spectrum_bypasses_the_spectrum_hook(): +def test_supplied_spectrum_bypasses_the_spectrum_hook(monkeypatch): """ A supplied spectrum must not be routed through `_spectrum`. @@ -254,7 +254,7 @@ def test_supplied_spectrum_bypasses_the_spectrum_hook(): Global_CPD workflow does both). Whatever they did was already applied to the array the caller is holding, so doing it again would compound it. """ - grid, xc, yc, extent = _shared_spectrum_grid() + grid, xc, yc = _shared_spectrum_grid() calls = [] original = grid._spectrum @@ -262,25 +262,20 @@ def counting(*args, **kwargs): calls.append(1) return original(*args, **kwargs) - grid._spectrum = counting - try: - grid.optimise(200e3, xc, yc) - assert len(calls) == 1 - grid.optimise(200e3, xc, yc, spectrum=grid.last_spectrum) - assert len(calls) == 1, "_spectrum was called for a supplied spectrum" - finally: - grid._spectrum = original + monkeypatch.setattr(grid, "_spectrum", counting) + grid.optimise(200e3, xc, yc) + assert len(calls) == 1 + grid.optimise(200e3, xc, yc, spectrum=grid.last_spectrum) + assert len(calls) == 1, "_spectrum was called for a supplied spectrum" def test_last_spectrum_tracks_both_paths(): - # a fresh instance, because `synthetic_grid` is cached and a grid that has - # already been fitted carries the spectrum from whichever test got there - # first - data, extent = pycurious.fractal_anomaly(n=128, dx=2.0, beta=3.0, zt=1.0, - dz=20.0, C=5.0, seed=1) - grid = pycurious.CurieOptimiseBouligand(data, *extent) - xc = 0.5 * (extent[0] + extent[1]) - yc = 0.5 * (extent[2] + extent[3]) + # `__wrapped__` is the uncached `synthetic_grid`. A fresh instance is the + # point of the test: the cached one carries the spectrum from whichever + # test reached it first, so it would never start at None + grid, xc, yc = synthetic_grid.__wrapped__( + pycurious.CurieOptimiseBouligand, n=128, dx=2.0 + )[:3] assert grid.last_spectrum is None grid.optimise(200e3, xc, yc) @@ -292,23 +287,40 @@ def test_last_spectrum_tracks_both_paths(): assert np.array_equal(lhs, rhs) +def test_a_used_grid_still_pickles(): + """ + `parallelise_routine` pickles the bound method, and so the instance behind + it, to reach a worker. Anything `optimise` leaves on the instance therefore + has to survive a round trip -- including the provenance riding on + `last_spectrum`, and including a `process_subgrid` the caller passed, which + must not be retained at all. + """ + import pickle + + grid, xc, yc = _shared_spectrum_grid() + grid.optimise(200e3, xc, yc, process_subgrid=lambda subgrid: subgrid) + + restored = pickle.loads(pickle.dumps(grid)) + for lhs, rhs in zip(restored.last_spectrum, grid.last_spectrum): + assert np.array_equal(lhs, rhs) + assert restored.last_spectrum.provenance == grid.last_spectrum.provenance + + def test_supplied_spectrum_rejects_mismatched_shapes(): - grid, xc, yc, extent = _shared_spectrum_grid() + grid, xc, yc = _shared_spectrum_grid() bad = (np.ones(5), np.ones(4), np.ones(5)) with pytest.raises(ValueError, match="same shape"): grid.optimise(200e3, xc, yc, spectrum=bad) -def test_residuals_do_not_swallow_unrelated_warnings(): +def test_residuals_do_not_swallow_unrelated_warnings(monkeypatch): """ `residuals` silences the forward model's floating-point errors, which it must -- the fit legitimately probes dz <= 0 -- but it used to do so with a blanket `catch_warnings`, which ate every other warning raised in the block as well. """ - import warnings - - grid, xc, yc, extent = _shared_spectrum_grid() + grid, xc, yc = _shared_spectrum_grid() k, Phi, sigma = grid.window_spectrum(200e3, xc, yc, power=2.0) # dz < 0 is what the Curie profile evaluates below zt, and what raises @@ -326,16 +338,9 @@ def noisy(*args, **kwargs): warnings.warn("a real warning", RuntimeWarning) return real(*args, **kwargs) - pycurious.optimise_bouligand.bouligand2009 = noisy - try: - with warnings.catch_warnings(record=True) as caught: - warnings.simplefilter("always") - grid.residuals(np.array([3.0, 1.0, 10.0, 5.0]), k, Phi, sigma) - assert [w for w in caught if "a real warning" in str(w.message)], ( - "residuals is still swallowing warnings it did not raise" - ) - finally: - pycurious.optimise_bouligand.bouligand2009 = real + monkeypatch.setattr(pycurious.optimise_bouligand, "bouligand2009", noisy) + with pytest.warns(RuntimeWarning, match="a real warning"): + grid.residuals(np.array([3.0, 1.0, 10.0, 5.0]), k, Phi, sigma) def test_supplied_spectrum_from_another_window_warns(): @@ -346,7 +351,7 @@ def test_supplied_spectrum_from_another_window_warns(): 128 km spectrum passed to a 512 km call returns dz = 222.8 km where that window really gives 21.7. """ - grid, xc, yc, extent = _shared_spectrum_grid() + grid, xc, yc = _shared_spectrum_grid() grid.optimise(100e3, xc, yc) ours = grid.last_spectrum @@ -359,36 +364,28 @@ def test_supplied_spectrum_from_another_window_warns(): grid.optimise(100e3, xc + 40e3, yc, spectrum=grid.last_spectrum) -def test_supplied_spectrum_at_matching_arguments_is_silent(): +def test_supplied_spectrum_at_matching_arguments_is_silent(recwarn): """The documented idiom must not warn, or the guard is worse than useless.""" - import warnings - - grid, xc, yc, extent = _shared_spectrum_grid() + grid, xc, yc = _shared_spectrum_grid() grid.optimise(200e3, xc, yc) - with warnings.catch_warnings(record=True) as caught: - warnings.simplefilter("always") - grid.profile(200e3, xc, yc, "dz", npoints=7, spectrum=grid.last_spectrum) - assert not [w for w in caught if "supplied spectrum" in str(w.message)] + grid.profile(200e3, xc, yc, "dz", npoints=7, spectrum=grid.last_spectrum) + assert not [w for w in recwarn if "supplied spectrum" in str(w.message)] -def test_spectrum_of_unknown_provenance_is_taken_at_face_value(): +def test_spectrum_of_unknown_provenance_is_taken_at_face_value(recwarn): """ - Only a spectrum this instance computed can be checked. One built by the + Only a spectrum this library built carries provenance. One assembled by the caller -- Global_CPD reads its archived spectra out of zarr as float32 -- has nothing to compare against, so it must be accepted without a warning rather than guessed at. """ - import warnings - - grid, xc, yc, extent = _shared_spectrum_grid() + grid, xc, yc = _shared_spectrum_grid() grid.optimise(100e3, xc, yc) foreign = tuple(np.asarray(a).astype(np.float32) for a in grid.last_spectrum) - with warnings.catch_warnings(record=True) as caught: - warnings.simplefilter("always") - grid.optimise(200e3, xc, yc, spectrum=foreign) - assert not [w for w in caught if "supplied spectrum" in str(w.message)] + grid.optimise(200e3, xc, yc, spectrum=foreign) + assert not [w for w in recwarn if "supplied spectrum" in str(w.message)] def test_provenance_survives_being_passed_along(): @@ -396,7 +393,7 @@ def test_provenance_survives_being_passed_along(): Reusing a spectrum must not relabel it with the arguments of whichever call reused it, or the guard would go blind after one hop. """ - grid, xc, yc, extent = _shared_spectrum_grid() + grid, xc, yc = _shared_spectrum_grid() grid.optimise(100e3, xc, yc) grid.profile(100e3, xc, yc, "dz", npoints=7, spectrum=grid.last_spectrum) diff --git a/tests/test_routines.py b/tests/test_routines.py index 874c8b9..6ad779a 100644 --- a/tests/test_routines.py +++ b/tests/test_routines.py @@ -59,3 +59,27 @@ def test_CurieOptimiseBouligand_routines(load_magnetic_anomaly): assert time_routine(cpd.optimise_routine, 0.5 * max_window, xc_list, yc_list) assert time_routine(cpd.metropolis_hastings, 0.5 * max_window, xc, yc, 100, 10) assert time_routine(cpd.sensitivity, 0.5 * max_window, xc, yc, 100) + + +def test_a_routine_refuses_a_shared_spectrum(load_magnetic_anomaly): + """ + One spectrum cannot describe a list of centroids. + + Forwarded, it would give every centroid the same answer -- a flat map that + looks like a result -- and the per-call provenance warning fires only on + the serial path, so whether the user was told would depend on how many + processors happened to be available. + """ + d = load_magnetic_anomaly["mag_data"] + xc = load_magnetic_anomaly["xc"] + yc = load_magnetic_anomaly["yc"] + xmin, xmax, ymin, ymax = load_magnetic_anomaly["extent"] + max_window = load_magnetic_anomaly["max_window"] + + cpd = pycurious.CurieOptimiseBouligand(d, xmin, xmax, ymin, ymax) + cpd.optimise(0.5 * max_window, xc, yc) + + with pytest.raises(ValueError, match="single window at a single centroid"): + cpd.optimise_routine( + 0.5 * max_window, [xc, xc], [yc, yc], spectrum=cpd.last_spectrum + ) From 69e22dfe35ed7f920bb7ad662d6b5e8247832686 Mon Sep 17 00:00:00 2001 From: Ben Mather Date: Mon, 3 Aug 2026 13:54:04 +1000 Subject: [PATCH 2/2] Give Tanaka the spectrum sharing the Bouligand side had `_resolve_spectrum`, `_Spectrum` and `last_spectrum` sat on `CurieOptimiseBouligand`, so `CurieOptimiseTanaka` recomputed a spectrum on every call. It has the stronger case of the two: two straight-line fits cost almost nothing beside computing the spectrum -- 24.8 ms of a 25.9 ms `optimise` at a 1025-cell window, 96% -- and the documented workflow is a band sweep, where `check_bands` needs `k`, then `optimise` needs all of it, then `sensitivity`, then a revised band needs it again. None of those depend on the bands. Measured on `optimise` / `check_bands` / `optimise`, 51.9 ms -> 26.6. The mechanism was already class-agnostic -- it never mentions `power`, which each `_spectrum` pins below the seam -- so it moves to `CurieGrid` beside `window_spectrum`. Two things were Bouligand-specific and both were parameters in disguise. A subclass defining `_spectrum` now also defines: - `_SPECTRUM_ARGS`, the names of its positional arguments in order, so the provenance is picked out by name rather than by position. Tanaka takes seven with `beta` inserted at position five, Bouligand six. - `_SPECTRUM_PROVENANCE`, the subset a reused spectrum is checked against. Tanaka's includes `beta`, which subtracts the fractal contribution and so changes the values of `Phi`: a spectrum computed at one beta is the wrong data for a fit at another, and is now a warning rather than a silent answer. - `_SPECTRUM_RETURNS`, which fixes the arity -- four arrays here against three -- and names them when it is wrong, which is worth having where the order of `(k, Phi, Phi_n, sigma)` is not obvious. A call site is unchanged apart from the name: it forwards exactly what it forwarded before, in `_spectrum`'s own order. `optimise` and `sensitivity` take `spectrum=`. `parallelise_routine` already rejects it and is inherited, so the Tanaka routines get that for free -- verified rather than assumed. The three names are a second list to keep in step with a signature, which is the drift this refactor exists to remove, so `test_grid.py` pins them against `inspect.signature` and asserts no callable reaches the provenance. Suite 123 -> 129 tests. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01PLgEaDtuZw6QXDphSrzcdQ --- CLAUDE.md | 18 ++++- pycurious/grid.py | 115 ++++++++++++++++++++++++++++++++ pycurious/optimise_bouligand.py | 110 ++---------------------------- pycurious/optimise_tanaka.py | 50 ++++++++++++-- tests/test_grid.py | 26 ++++++++ tests/test_tanaka.py | 78 ++++++++++++++++++++++ 6 files changed, 286 insertions(+), 111 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index d9d3192..af562f5 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -15,8 +15,8 @@ returns a bare number without one is a pre-v2 remnant. ## Commands ```bash -pytest # 123 tests, ~13 s -pytest -m "not slow" # 120 tests, ~6 s -- skips the calibration tests that +pytest # 129 tests, ~13 s +pytest -m "not slow" # 126 tests, ~6 s -- skips the calibration tests that # fit a few hundred realisations pytest tests/test_tanaka.py -q ``` @@ -118,6 +118,20 @@ annulus**. That is not what a fit needs. `window_spectrum` converts it to the **uncertainty of the annulus mean** and is what both optimisers call (through their private `_spectrum`). Prefer it over `radial_spectrum` for anything fitted. +Computing that spectrum is most of what a fit costs — 96% of a Tanaka +`optimise` at a 1025-cell window — so every fitting routine takes `spectrum=` +and sets `last_spectrum`, and one spectrum serves a whole sweep at a centroid. +`CurieGrid._resolve_spectrum` is the single seam: it calls the subclass's +`_spectrum` or takes the caller's, never both. A subclass supplying `_spectrum` +also supplies `_SPECTRUM_ARGS`, `_SPECTRUM_PROVENANCE` and `_SPECTRUM_RETURNS`, +which is what lets one implementation serve Bouligand's 3-array return and +Tanaka's 4-array one. Two rules that look arbitrary and are not: **callables +never go in the provenance** (holding a `process_subgrid` on the instance makes +the bound routine unpicklable, which silently drops `parallelise_routine` to +serial), and **anything that changes the spectrum's values does** — including +Tanaka's `beta`, which subtracts the fractal contribution. `parallelise_routine` +rejects `spectrum` outright; one spectrum cannot describe a list of centroids. + Two corrections, both measured by Monte Carlo rather than assumed: 1. **Within a bin**, cells are not independent — a real field is Hermitian, so diff --git a/pycurious/grid.py b/pycurious/grid.py index 7a83048..0e02dd4 100644 --- a/pycurious/grid.py +++ b/pycurious/grid.py @@ -190,6 +190,34 @@ def _dof_factor(taper, counts=None, dof_factor=None): return dof_inf * counts / np.maximum(counts - lost, 1.0) +class _Spectrum(tuple): + """ + A spectrum that remembers what it was computed from. + + It is an ordinary tuple everywhere it matters -- it unpacks, indexes, zips + and pickles like the one `CurieGrid.window_spectrum` returns -- so + `spectrum=` still takes a plain tuple and `last_spectrum` still hands one + back. Carrying the provenance on the spectrum rather than beside it means + it survives being passed from routine to routine, and that a tuple built + anywhere else simply has none, which is exactly the "nothing to check, take + it as given" case. + + The arity is whatever the subclass's `_spectrum` returns: three arrays on + the Bouligand side, four on the Tanaka one. + """ + + def __new__(cls, arrays, provenance): + spectrum = super(_Spectrum, cls).__new__(cls, arrays) + spectrum.provenance = provenance + return spectrum + + def __getnewargs__(self): + # a tuple subclass with a two-argument `__new__` cannot be unpickled + # without this, and a grid carrying a `last_spectrum` is pickled every + # time `parallelise_routine` sends one to a worker + return (tuple(self), self.provenance) + + class CurieGrid(CurieParallel): """ Accepts a 2D array and Cartesian coordinates specifying the @@ -252,6 +280,10 @@ def __init__(self, grid, xmin, xmax, ymin, ymax, **kwargs): self.nx, self.ny = nx, ny self.dx, self.dy = dx, dy + # the spectrum most recently computed or supplied, ready to be handed + # to another routine at the same centroid -- see `_resolve_spectrum` + self.last_spectrum = None + if not np.allclose(dx, dy, 1.0): raise ValueError("node spacing should be identical {}".format((dx, dy))) @@ -679,6 +711,89 @@ def process_subgrid(subgrid): return k, Phi, sigma + def _resolve_spectrum(self, spectrum, *args, **kwargs): + """ + The spectrum a fitting routine should use: the caller's, or a fresh one. + + Every fitting routine in both optimisers begins by turning a window + into a spectrum, and computing it is most of what a call at a large + window costs -- 30% of a Bouligand `optimise` plus `profile` pair at a + 1025-cell window, 41% at 2049, and 81% of a single Tanaka `optimise`, + whose two straight-line fits cost almost nothing beside it. Passing one + spectrum to several routines removes that entirely. + + `args` are `_spectrum`'s own positional arguments, in its own order, so + a routine forwards exactly what it would have forwarded anyway. A + subclass that defines `_spectrum` also defines the three names below; + they are what lets one implementation serve signatures that differ: + + - `_SPECTRUM_ARGS` -- names of those positional arguments, in order + - `_SPECTRUM_PROVENANCE` -- the subset a reused spectrum is checked + against. `taper` and `process_subgrid` are deliberately excluded: + they are callables, and holding one on the instance both keeps its + captured scope alive and makes the bound routine unpicklable, which + drops `pycurious.parallel.CurieParallel.parallelise_routine` back to + serial with a warning that blames the wrong thing. Anything that + changes the *values* of the spectrum belongs here -- including + Tanaka's `beta`, which subtracts the fractal contribution. + - `_SPECTRUM_RETURNS` -- names of the arrays `_spectrum` returns, which + fixes the arity and names them in the error when it is wrong + + A supplied spectrum bypasses `_spectrum` rather than being routed + through it. A subclass that overrides `_spectrum` to band limit, or to + reweight `sigma`, has already applied that to the array the caller is + holding; applying it a second time would compound it. + + Nothing is reused implicitly: a routine given no `spectrum` always + computes one, and the library never reads `last_spectrum` itself. But a + supplied spectrum makes `window`, `xc` and `yc` dead arguments, so + handing over the one from a *different* window is accepted in silence + and answers a question the caller did not ask. `_Spectrum.provenance` + guards the case that can be guarded: a spectrum this library computed + knows what it came from, and disagreeing with it is a warning. One + built anywhere else -- read from an archive, cast to float32 -- has no + provenance, so it is taken at face value rather than guessed at. + """ + named = dict(zip(self._SPECTRUM_ARGS, args)) + provenance = tuple(named[name] for name in self._SPECTRUM_PROVENANCE) + + if spectrum is None: + spectrum = _Spectrum(self._spectrum(*args, **kwargs), provenance) + else: + was = getattr(spectrum, "provenance", None) + + arrays = tuple(np.asarray(a, dtype=float) for a in spectrum) + names = self._SPECTRUM_RETURNS + if len(arrays) != len(names) or len({a.shape for a in arrays}) != 1: + raise ValueError( + "spectrum must be {} arrays of the same shape, ({}), got " + "{}".format( + len(names), + ", ".join(names), + ", ".join(str(a.shape) for a in arrays), + ) + ) + spectrum = _Spectrum(arrays, was) + + mismatched = [] if was is None else [ + name + for name, then, now in zip(self._SPECTRUM_PROVENANCE, was, provenance) + if then != now + ] + if mismatched: + names = ", ".join(mismatched) + warnings.warn( + "the supplied spectrum was computed with a different {0}; " + "a supplied spectrum is used as given, so the {0} of this " + "call is ignored and the result describes the window the " + "spectrum came from, not the one asked for here.".format(names), + RuntimeWarning, + stacklevel=3, + ) + + self.last_spectrum = spectrum + return spectrum + def reduce_to_pole(self, data, inc, dec, sinc=None, sdec=None): """ Reduce total field magnetic anomaly data to the pole. diff --git a/pycurious/optimise_bouligand.py b/pycurious/optimise_bouligand.py index f534e91..4b81989 100644 --- a/pycurious/optimise_bouligand.py +++ b/pycurious/optimise_bouligand.py @@ -107,40 +107,6 @@ # forward model but a sum of two of them. _CPD = "CPD" -# The arguments a reused spectrum is checked against, in the order -# `_Spectrum.provenance` stores them. `taper` and `process_subgrid` are -# deliberately absent: they are callables, and holding one on the instance both -# keeps its captured scope alive and makes `grid.optimise` unpicklable, which -# drops `pycurious.parallel.CurieParallel.parallelise_routine` back to serial -# with a warning that blames the wrong thing. The misuse worth catching is a -# spectrum from a different window or centroid anyway. -_SPECTRUM_ARGS = ("window", "xc", "yc", "dof_factor") - - -class _Spectrum(tuple): - """ - A `(k, Phi, sigma_Phi)` triple that remembers what it was computed from. - - It is an ordinary tuple everywhere it matters -- it unpacks, indexes, zips - and pickles like the one `pycurious.grid.CurieGrid.window_spectrum` - returns -- so `spectrum=` still takes a plain triple and `last_spectrum` - still hands one back. Carrying the provenance on the spectrum rather than - beside it means it survives being passed from routine to routine, and that - a triple built anywhere else simply has none, which is exactly the - "nothing to check, take it as given" case. - """ - - def __new__(cls, arrays, provenance): - spectrum = super(_Spectrum, cls).__new__(cls, arrays) - spectrum.provenance = provenance - return spectrum - - def __getnewargs__(self): - # a tuple subclass with a two-argument `__new__` cannot be unpickled - # without this, and a grid carrying a `last_spectrum` is pickled every - # time `parallelise_routine` sends one to a worker - return (tuple(self), self.provenance) - def _prior_loc_scale(pdf): """ @@ -249,10 +215,6 @@ def __init__(self, grid, xmin, xmax, ymin, ymax, **kwargs): ub = [None, None, self._max_thickness(), None] self.bounds = list(zip(lb, ub)) - # the spectrum most recently computed or supplied, ready to be handed - # to another routine at the same centroid -- see `_resolve_spectrum` - self.last_spectrum = None - self.max_processors = kwargs.pop("max_processors", cpu_count()) def _max_thickness(self): @@ -446,6 +408,11 @@ def min_func(self, x, kh, Phi, sigma_Phi, prior=None): """ return 0.5 * np.sum(self.residuals(x, kh, Phi, sigma_Phi, prior) ** 2) + # see `pycurious.grid.CurieGrid._resolve_spectrum` + _SPECTRUM_ARGS = ("window", "xc", "yc", "taper", "process_subgrid", "dof_factor") + _SPECTRUM_PROVENANCE = ("window", "xc", "yc", "dof_factor") + _SPECTRUM_RETURNS = ("k", "Phi", "sigma_Phi") + def _spectrum(self, window, xc, yc, taper, process_subgrid, dof_factor, **kwargs): """ Radial power spectrum of one window, weighted ready for fitting. @@ -464,73 +431,6 @@ def _spectrum(self, window, xc, yc, taper, process_subgrid, dof_factor, **kwargs **kwargs ) - def _resolve_spectrum( - self, spectrum, window, xc, yc, taper, process_subgrid, dof_factor, **kwargs - ): - """ - The spectrum a fitting routine should use: the caller's, or a fresh one. - - `optimise`, `profile`, `sensitivity` and `metropolis_hastings` all begin - by turning a window into a spectrum, and computing it is most of what a - call at a large window costs -- 30% of an `optimise` plus `profile` pair - at a 1025-cell window, 41% at 2049. Passing the same spectrum to both - removes that entirely. - - A supplied spectrum bypasses `_spectrum` rather than being routed - through it. A subclass that overrides `_spectrum` to band limit, or to - reweight `sigma`, has already applied that to the array the caller is - holding; applying it a second time would compound it. - - Nothing is reused implicitly: a routine given no `spectrum` always - computes one, and the library never reads `last_spectrum` itself. But a - supplied spectrum makes `window`, `xc` and `yc` dead arguments, so - handing over the one from a *different* window is accepted in silence - and answers a question the caller did not ask. `_Spectrum.provenance` - guards the case that can be guarded: a spectrum this library computed - knows the `_SPECTRUM_ARGS` it came from, and disagreeing with them is a - warning. One built anywhere else -- read from an archive, cast to - float32 -- has no provenance, so it is taken at face value rather than - guessed at. - """ - provenance = (window, xc, yc, dof_factor) - - if spectrum is None: - spectrum = _Spectrum( - self._spectrum( - window, xc, yc, taper, process_subgrid, dof_factor, **kwargs - ), - provenance, - ) - else: - was = getattr(spectrum, "provenance", None) - - k, Phi, sigma_Phi = (np.asarray(a, dtype=float) for a in spectrum) - if not (k.shape == Phi.shape == sigma_Phi.shape): - raise ValueError( - "spectrum must be three arrays of the same shape, got " - "{}, {} and {}".format(k.shape, Phi.shape, sigma_Phi.shape) - ) - spectrum = _Spectrum((k, Phi, sigma_Phi), was) - - mismatched = [] if was is None else [ - name - for name, then, now in zip(_SPECTRUM_ARGS, was, provenance) - if then != now - ] - if mismatched: - names = ", ".join(mismatched) - warnings.warn( - "the supplied spectrum was computed with a different {0}; " - "a supplied spectrum is used as given, so the {0} of this " - "call is ignored and the result describes the window the " - "spectrum came from, not the one asked for here.".format(names), - RuntimeWarning, - stacklevel=3, - ) - - self.last_spectrum = spectrum - return spectrum - def _bound_arrays(self, free=None): """ `self.bounds` as a pair of arrays, with `None` meaning infinite. diff --git a/pycurious/optimise_tanaka.py b/pycurious/optimise_tanaka.py index 1c20b3f..8e45b3b 100644 --- a/pycurious/optimise_tanaka.py +++ b/pycurious/optimise_tanaka.py @@ -151,6 +151,12 @@ def check_bands(self, k, zt_range, z0_range, thickness=None, verbose=True): Usage: >>> k, Phi, sigma_Phi = grid.radial_spectrum(subgrid, power=1) >>> grid.check_bands(k, (0.2, 0.6), (0.0, 0.05), thickness=20.0) + + After an `optimise` at the same centroid, take `k` from the + spectrum it already computed rather than building another:: + + >>> grid.check_bands(grid.last_spectrum[0], (0.2, 0.6), + ... (0.0, 0.05), thickness=20.0) """ k = np.asarray(k) _warn_if_cycles_per_km(zt_range, k) @@ -310,6 +316,15 @@ def _correlated_gradient_stdev(k, Phi, sigma, gradient, intercept, fallback): stdev = np.sqrt(np.diag(cov))[0] return stdev if np.isfinite(stdev) else fallback + # see `pycurious.grid.CurieGrid._resolve_spectrum`. `beta` is provenance + # because it subtracts the fractal contribution from `Phi`, so a spectrum + # computed at one `beta` is the wrong data for a fit at another. + _SPECTRUM_ARGS = ( + "window", "xc", "yc", "taper", "beta", "process_subgrid", "dof_factor" + ) + _SPECTRUM_PROVENANCE = ("window", "xc", "yc", "beta", "dof_factor") + _SPECTRUM_RETURNS = ("k", "Phi", "Phi_n", "sigma") + def _spectrum(self, window, xc, yc, taper, beta, process_subgrid, dof_factor, **kwargs): """ @@ -354,6 +369,7 @@ def optimise( process_subgrid=None, absolute_sigma=True, dof_factor=None, + spectrum=None, **kwargs ): """ @@ -388,6 +404,11 @@ def optimise( dof_factor : float, optional override the effective-degrees-of-freedom deflation applied to the spectral uncertainties (see Notes) + spectrum : tuple (k, Phi, Phi_n, sigma), optional + a spectrum already in hand, typically `last_spectrum` from an + earlier call at this same centroid. Skips computing one, and + `window`, `xc`, `yc`, `taper`, `beta`, `process_subgrid` and + `dof_factor` are then unused -- see Notes. kwargs : keyword arguments passed to `radial_spectrum` @@ -419,6 +440,22 @@ def optimise( only. They do not include the systematic error from the choice of band, which is usually larger -- see `sensitivity`. + Choosing those bands means fitting the same spectrum several times, + and the spectrum does not depend on them: `check_bands` needs `k`, + then `optimise` and `sensitivity` need all of it, then a revised + band needs it again. Two straight-line fits cost almost nothing + beside computing it -- 24.8 ms of a 25.9 ms `optimise` at a + 1025-cell window, 96% -- so compute it once and pass it along:: + + >>> grid.optimise(200e3, xc, yc, (0.2, 0.6), (0.0, 0.05)) + >>> spectrum = grid.last_spectrum + >>> grid.check_bands(spectrum[0], (0.2, 0.6), (0.0, 0.05)) + >>> grid.optimise(200e3, xc, yc, (0.25, 0.6), (0.0, 0.04), + ... spectrum=spectrum) + + which halves that sweep, 51.9 ms to 26.6. `sensitivity` takes the + same spectrum, and each further revision is then free. + `radial_spectrum` returns the scatter of the FFT cells within each annulus, whereas the fit needs the uncertainty of the annulus mean. That is the standard error, except that the cells are not @@ -449,8 +486,8 @@ def optimise( covariance by the reduced chi-squared afterwards, which the profile construction does not do. """ - k, Phi, Phi_n, sigma = self._spectrum( - window, xc, yc, taper, beta, process_subgrid, dof_factor, **kwargs + k, Phi, Phi_n, sigma = self._resolve_spectrum( + spectrum, window, xc, yc, taper, beta, process_subgrid, dof_factor, **kwargs ) _warn_if_cycles_per_km(zt_range, k) @@ -526,6 +563,7 @@ def sensitivity( absolute_sigma=True, dof_factor=None, seed=None, + spectrum=None, **kwargs ): """ @@ -554,6 +592,10 @@ def sensitivity( spectrum, which recovers the analytic covariance. seed : int, optional seed for reproducibility + spectrum : tuple (k, Phi, Phi_n, sigma), optional + a spectrum already in hand, typically `last_spectrum` from the + `optimise` at this same centroid -- see `optimise`. The bands + being sampled do not enter it, so one spectrum serves both. Returns: zt : 1D array shape (nsim,) @@ -580,8 +622,8 @@ def sensitivity( # the spectrum is computed once and resampled, as in # CurieOptimiseBouligand.sensitivity - k, Phi, Phi_n, sigma = self._spectrum( - window, xc, yc, taper, beta, process_subgrid, dof_factor, **kwargs + k, Phi, Phi_n, sigma = self._resolve_spectrum( + spectrum, window, xc, yc, taper, beta, process_subgrid, dof_factor, **kwargs ) _warn_if_cycles_per_km(zt_range, k) diff --git a/tests/test_grid.py b/tests/test_grid.py index 6ed8fa1..09d9496 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -441,3 +441,29 @@ def lstsq_detrend(data): detrended = grid.remove_trend_linear(line) assert np.all(np.isfinite(detrended)) np.testing.assert_allclose(detrended, lstsq_detrend(line), atol=1e-9) + + +@pytest.mark.parametrize( + "cls", [pycurious.CurieOptimiseBouligand, pycurious.CurieOptimiseTanaka] +) +def test_spectrum_arg_names_match_the_signature(cls): + """ + `_resolve_spectrum` names `_spectrum`'s positional arguments so it can pick + the provenance out of them by name. That is a second list to keep in step + with the signature, and a silent slip in it would put the wrong value under + the wrong name -- so pin it rather than trusting it. + """ + import inspect + + declared = list(cls._SPECTRUM_ARGS) + # drop `self`, and any keyword-only tail `_resolve_spectrum` never passes + actual = list(inspect.signature(cls._spectrum).parameters)[1:] + + assert actual[: len(declared)] == declared + assert set(cls._SPECTRUM_PROVENANCE) <= set(declared) + for callable_arg in ("taper", "process_subgrid"): + assert callable_arg not in cls._SPECTRUM_PROVENANCE, ( + "a callable in the provenance is retained on the instance, which " + "makes the bound routine unpicklable and drops parallelise_routine " + "to serial" + ) diff --git a/tests/test_tanaka.py b/tests/test_tanaka.py index 719c67e..026cd97 100644 --- a/tests/test_tanaka.py +++ b/tests/test_tanaka.py @@ -180,3 +180,81 @@ def test_CurieOptimiseTanaka_routines(load_magnetic_anomaly): CPD, sigma_CPD = grid.calculate_CPD(zt, z0, sigma_zt, sigma_z0) assert CPD.shape == (len(xc_list),) assert np.all(sigma_CPD > 0.0) + + +def test_supplied_spectrum_reproduces_the_computed_one(tanaka): + """ + `spectrum=` must be the same fit, not merely a similar one. + + Both routines route it through the same `_resolve_spectrum` the Bouligand + side uses, rather than getting their own code path, so the two cannot drift + apart. Assert to the bit, so that they cannot. + """ + grid, xc, yc = tanaka + window = 300e3 + + stock = grid.optimise(window, xc, yc, ZT_RANGE, Z0_RANGE) + shared = grid.optimise( + window, xc, yc, ZT_RANGE, Z0_RANGE, spectrum=grid.last_spectrum + ) + assert stock == shared + + a = grid.sensitivity(window, xc, yc, 20, ZT_RANGE, Z0_RANGE, seed=1) + b = grid.sensitivity( + window, xc, yc, 20, ZT_RANGE, Z0_RANGE, seed=1, spectrum=grid.last_spectrum + ) + for lhs, rhs in zip(a, b): + assert np.array_equal(lhs, rhs) + + +def test_last_spectrum_carries_all_four_arrays(tanaka): + """ + Tanaka's `_spectrum` returns `(k, Phi, Phi_n, sigma)`, not the triple the + Bouligand side returns. The shared machinery must not have assumed three. + """ + grid, xc, yc = tanaka + assert grid.last_spectrum is None + + grid.optimise(300e3, xc, yc, ZT_RANGE, Z0_RANGE) + assert len(grid.last_spectrum) == 4 + + bad = (np.ones(5), np.ones(5), np.ones(5)) + with pytest.raises(ValueError, match="4 arrays of the same shape"): + grid.optimise(300e3, xc, yc, ZT_RANGE, Z0_RANGE, spectrum=bad) + + +def test_a_spectrum_from_a_different_beta_warns(tanaka): + """ + `beta` subtracts the fractal contribution from `Phi`, so a spectrum + computed at one beta is the wrong data for a fit at another -- which makes + it provenance, exactly like the window and the centroid. + """ + grid, xc, yc = tanaka + grid.optimise(300e3, xc, yc, ZT_RANGE, Z0_RANGE, beta=3.0) + + with pytest.warns(RuntimeWarning, match="different beta"): + grid.optimise( + 300e3, xc, yc, ZT_RANGE, Z0_RANGE, beta=2.0, spectrum=grid.last_spectrum + ) + + grid.optimise(300e3, xc, yc, ZT_RANGE, Z0_RANGE, beta=3.0) + with pytest.warns(RuntimeWarning, match="different window"): + grid.optimise( + 200e3, xc, yc, ZT_RANGE, Z0_RANGE, beta=3.0, spectrum=grid.last_spectrum + ) + + +def test_bands_are_not_provenance(tanaka, recwarn): + """ + The spectrum does not depend on the fitting bands, so sweeping them while + reusing one spectrum -- the whole point of the idiom `optimise` documents + -- must not warn. + """ + grid, xc, yc = tanaka + grid.optimise(300e3, xc, yc, ZT_RANGE, Z0_RANGE) + + grid.check_bands(grid.last_spectrum[0], ZT_RANGE, Z0_RANGE, verbose=False) + grid.optimise( + 300e3, xc, yc, (1.3, 1.8), (0.0, 0.5), spectrum=grid.last_spectrum + ) + assert not [w for w in recwarn if "supplied spectrum" in str(w.message)]