diff --git a/CLAUDE.md b/CLAUDE.md index 27a9070..5b560d9 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 # 70 tests, ~40 s -pytest -m "not slow" # 67 tests, ~26 s -- skips the calibration tests that +pytest # 111 tests, ~20 s +pytest -m "not slow" # 108 tests, ~12 s -- skips the calibration tests that # fit a few hundred realisations pytest tests/test_tanaka.py -q ``` @@ -78,6 +78,25 @@ depth. **Depths are positive downwards.** `optimise` returns depths, not the negative gradients the fits produce. +**Every Bouligand fit goes through `_fit`, which is `least_squares`, not +`minimize`.** Do not put L-BFGS-B back: it needs **2012** evaluations of the +forward model per vertex where trust-region reflective needs 359, about 3x the +CPU. `_fit` also supplies the two exact Jacobian columns (`dr/dzt = -2k/sigma`, +`dr/dC = 1/sigma`); `beta` enters through the *order* of a Bessel function and +`dz` costs the same analytically as by difference, so those two stay numerical. + +**Pin the BLAS before timing anything.** + +```bash +OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1 python bench.py +``` + +L-BFGS-B calls a threaded BLAS whose workers spin-wait. Unpinned on a loaded +machine it looks 20x to 400x slower than TRF depending on contention, and +`time.process_time` makes it worse by charging every spinning thread — a +"3.1 s CPU" fit whose wall clock was 1.7 s. Both numbers were wrong; the truth +is 3x. Counting `residuals` calls is the measurement that does not lie. + **Tanaka bands have no defaults, deliberately.** Each straight-line limit holds only over part of the spectrum: the `zt` band needs wavelengths shorter than ~4x the source thickness, the `z0` band needs `|k|d << 1`. Violating either @@ -202,6 +221,17 @@ left in place, which makes a correct `MANIFEST.in` look broken. ## Known defects +- **`profile` reports one basin of a multimodal deviance.** The scan walks + outward from the best node to the first threshold crossing, so where the + misfit has two minima it covers the one around the best node and never sees + the other, and the interval can then exclude the fitted value. Measured on + synthetics at a 200 km window, 2 of 20 intervals did. It does **not** appear + in the regime that matters: over cached EMAG2 spectra in `~/Global_CPD` — + three window sizes from 1000 to 4000 km, three targets, twelve mesh vertices + — 0 of 108 intervals excluded their estimate. Band limiting and the prior + pinning `zt` between them seem to remove it. Worth knowing if you profile a + small unconstrained synthetic; not worth guarding against. + - **`install_documentation()` fails for an installed package.** It is still advertised in the README, but `[tool.setuptools] packages = ["pycurious"]` installs only the package directory, and `Examples/` sits at the diff --git a/pycurious/optimise_bouligand.py b/pycurious/optimise_bouligand.py index 3631ece..afd017e 100644 --- a/pycurious/optimise_bouligand.py +++ b/pycurious/optimise_bouligand.py @@ -32,7 +32,7 @@ from .parallel import stochastic import numpy as np import warnings -from scipy.optimize import minimize, brentq +from scipy.optimize import least_squares, brentq from scipy import stats from multiprocessing import cpu_count @@ -62,9 +62,37 @@ # Parameters of bouligand2009, in the order the optimiser sees them. _PARAMETERS = ("beta", "zt", "dz", "C") -# Relative step for the finite-difference Jacobian behind the covariance. +# Relative step for the finite-difference Jacobian, both in `_fit` and behind +# the covariance. _JACOBIAN_STEP = 1.0e-6 +# `zt`, which the Curie constraint dz = CPD - zt differentiates through. +_ZT = _PARAMETERS.index("zt") + +# Width given to a bound the caller collapsed to a point. Relative, so it pins +# a parameter of any magnitude to its own last few bits. +_DEGENERATE_BOUND = 1.0e-12 + +# Two of the four columns of the fit Jacobian are exact. Writing the forward +# model out, +# +# Phi_syn = C - 2 |k| zt - (beta - 1) ln|k| - |k| dz + ln A(beta, dz, |k|) +# +# `zt` and `C` enter linearly and A depends on neither, so the derivatives of +# the whitened residual (Phi_syn - Phi)/sigma are closed form. Supplying them +# removes two of the four finite differences per Jacobian, and in `profile` +# -- where `dz` is the coordinate being held -- two of the three. +# +# The other two stay numerical. `beta` enters through the *order* of the Bessel +# term, which has no closed-form derivative. `dz` enters through its argument +# and could be differentiated via d/dx K_v = -(K_{v-1} + K_{v+1})/2, one extra +# `kv` call given the recurrence -- but that is exactly what the finite +# difference it would replace costs, so there is nothing to gain. +_ANALYTIC_COLUMNS = { + "zt": lambda kh, sigma: -2.0 * kh / sigma, + "C": lambda kh, sigma: 1.0 / sigma, +} + # Acceptance rate the burn-in tunes the proposal towards, the usual optimum # for a random walk over a smooth multivariate target. _TARGET_ACCEPTANCE = 0.234 @@ -116,7 +144,7 @@ def _prior_loc_scale(pdf): class CurieOptimiseBouligand(CurieGrid): """ Extends the `pycurious.grid.CurieGrid` class to include - optimisation routines see `scipy.optimize.minimize` for + optimisation routines see `scipy.optimize.least_squares` for a description of the algorithm. Args: @@ -385,6 +413,189 @@ def _spectrum(self, window, xc, yc, taper, process_subgrid, dof_factor, **kwargs **kwargs ) + def _bound_arrays(self, free=None): + """ + `self.bounds` as a pair of arrays, with `None` meaning infinite. + + Both the fit and the sampler need the bounds this way, and the `None` + sentinel is the sort of convention that goes stale in one copy. + + Args: + free : list of int, optional + indices of `_PARAMETERS` to include, default all four + + Returns: + lower, upper : arrays shape (len(free),) + """ + if free is None: + free = range(len(_PARAMETERS)) + pairs = [self.bounds[i] for i in free] + return ( + np.array([-np.inf if lo is None else lo for lo, _ in pairs], dtype=float), + np.array([np.inf if hi is None else hi for _, hi in pairs], dtype=float), + ) + + def _fit(self, y0, args, prior=None, free=None, fixed=None, curie=False): + """ + Bounded least-squares fit of `residuals`, with the Jacobian supplied. + + Every fit in this module goes through here, so they cannot drift apart + on method, bounds or Jacobian. + + Args: + y0 : array shape (len(free),) + starting point, in the free coordinates + args : tuple + `(kh, Phi, sigma_Phi)`, as `residuals` takes them + prior : dict, optional + priors to use in place of `self.prior` + free : list of int, optional + indices of `_PARAMETERS` allowed to vary, default all four + fixed : tuple, optional + `(index, value)` for a coordinate held constant, as `profile` + holds one + curie : bool (default=False) + the held coordinate is :math:`\\Delta z = \\mathrm{CPD} - z_t` + rather than a constant, so it moves with :math:`z_t` + + Returns: + res : `scipy.optimize.OptimizeResult` + `res.cost` is `min_func` at the solution by construction -- + both are half the sum of squares of the same vector. + + Notes: + This is `least_squares`, not `minimize`. The problem is a sum of + squares and `residuals` already exposes the vector, so the + trust-region reflective method can use structure L-BFGS-B cannot + see. Measured over six synthetic 200 km windows with the BLAS + pinned to one thread, `optimise` costs 12.5 ms of CPU against + 5.0 ms and a `dz` profile 102 ms against 33 ms -- about 3x -- on + **359** evaluations of the forward model per vertex rather than + 2012. + + The evaluation count is the durable number; the timings are not. + Unpinned, L-BFGS-B calls a threaded BLAS whose workers spin-wait, + and on a busy machine the same comparison reads anywhere from 20x + to 400x depending on how many cores are already contended. That is + an artefact of the measurement, not a property of the algorithms. + Pin `OMP_NUM_THREADS` and friends before timing anything here. + + Parameters agree with L-BFGS-B to three or four figures on windows + that constrain the fit, and exactly on real band-limited spectra. + Where the likelihood is flat they can disagree by much more -- + `dz` of 310 km against 139 km for a misfit difference of 2e-6 on + one 100 km synthetic -- because both answers are equally good and + neither optimiser has anything to descend. That is a property of + the window, not of the change. + """ + if free is None: + free = list(range(len(_PARAMETERS))) + if prior is None: + prior = self.prior + + kh, sigma = args[0], args[2] + nk = np.size(kh) + + # The map from the free coordinates to the full parameter vector is + # affine, `x = transform @ y + offset`, in all three cases the callers + # use. Writing it as a matrix once means `expand` and the chain rule in + # `jacobian` are the same statement rather than two hand-maintained + # ones -- an inconsistency between them is close to invisible in the + # output, since the fit converges from a wrong Jacobian anyway. + transform = np.zeros((len(_PARAMETERS), len(free))) + for column, index in enumerate(free): + transform[index, column] = 1.0 + offset = np.zeros(len(_PARAMETERS)) + if fixed is not None: + offset[fixed[0]] = fixed[1] + if curie: + # dz = CPD - zt, so the dz row differentiates through zt + transform[fixed[0], free.index(_ZT)] = -1.0 + + def expand(y): + return transform @ y + offset + + # `least_squares` evaluates the residual at a point and then asks for + # the Jacobian there, which is the only reuse there is: a rejected step + # is retried somewhere else, never at the same point. One entry. + last = [None, None] + + def residuals_at(x): + key = x.tobytes() + if last[0] != key: + last[0] = key + last[1] = self.residuals(x, *args, prior=prior) + return last[1] + + # Everything in the Jacobian except the differenced columns is constant + # over the fit, so build it once. The analytic columns depend only on + # `kh` and `sigma`; a prior contributes a row whose derivative is + # `1/scale` against its own parameter and zero elsewhere, in the order + # `residuals` appends them. Columns the transform discards are filled + # anyway -- harmless, and cheaper than deciding not to. + prior_rows = [ + index + for index, key in enumerate(_PARAMETERS) + if prior.get(key) is not None + ] + constant = np.zeros((nk + len(prior_rows), len(_PARAMETERS))) + for index, key in enumerate(_PARAMETERS): + analytic = _ANALYTIC_COLUMNS.get(key) + if analytic is not None: + constant[:nk, index] = analytic(kh, sigma) + for row, index in enumerate(prior_rows): + constant[nk + row, index] = 1.0 / prior[_PARAMETERS[index]][1] + + # what is left to difference: the free columns with no closed form, + # plus the held one when the Curie constraint makes it move with `zt` + differenced = [i for i in free if _PARAMETERS[i] not in _ANALYTIC_COLUMNS] + if curie: + differenced.append(fixed[0]) + + def jacobian(y): + x = expand(y) + J = constant.copy() + base = residuals_at(x) + for index in differenced: + step = _JACOBIAN_STEP * max(abs(x[index]), 1.0) + xp = x.copy() + xp[index] += step + # spectral rows only -- the prior rows are already exact in + # `constant`, and differencing them would count them twice + J[:nk, index] = (residuals_at(xp)[:nk] - base[:nk]) / step + return J @ transform + + lower, upper = self._bound_arrays(free) + + # L-BFGS-B accepted an equality bound and simply pinned the parameter. + # Trust-region reflective needs a box with an interior to reflect + # inside and refuses one without, so give it the narrowest box that is + # numerically distinct. That pins the parameter just as effectively, + # and leaves it near enough the edge that `_warn_on_bounds` still says + # the uncertainty there is meaningless. + # + # Equality only. An *inverted* bound is a typo, not an intention -- + # `self.bounds` is documented as reassignable, so it is a typo a user + # can make -- and widening it would silently rewrite it into whichever + # of the two numbers happened to be first. Let `least_squares` raise. + # An infinite bound cannot compare equal to the other, which is always + # its opposite sign, so this needs no finiteness test. + collapsed = upper == lower + if collapsed.any(): + upper[collapsed] = lower[collapsed] + _DEGENERATE_BOUND * np.maximum( + np.abs(lower[collapsed]), 1.0 + ) + + return least_squares( + lambda y: residuals_at(expand(y)), + # trust-region reflective requires a feasible start, which a + # caller passing the unconstrained solution into a constrained fit + # cannot guarantee + np.clip(np.asarray(y0, dtype=float), lower, upper), + jac=jacobian, + bounds=(lower, upper), + ) + def _jacobian(self, x, r, args): """ Central-difference Jacobian of `residuals` at `x`, given `r` there. @@ -520,9 +731,7 @@ def optimise( window, xc, yc, taper, process_subgrid, dof_factor, **kwargs ) - # minimise function - res = minimize(self.min_func, x0, args=(k, Phi, sigma_Phi), bounds=self.bounds) - x = res.x + x = self._fit(x0, (k, Phi, sigma_Phi)).x self._warn_on_bounds(x) cov = self._covariance(x, k, Phi, sigma_Phi) @@ -645,19 +854,15 @@ def _profiled_misfit(self, target, value, x_hat, args): fixed = _PARAMETERS.index("dz" if curie else target) free = [j for j in range(len(_PARAMETERS)) if j != fixed] - def expand(y): - x = np.empty(len(_PARAMETERS)) - x[free] = y - # for the Curie depth the constraint depends on zt, which is free - x[fixed] = value - x[1] if curie else value - return x - - res = minimize( - lambda y: self.min_func(expand(y), *args), + res = self._fit( np.asarray(x_hat)[free], - bounds=[self.bounds[j] for j in free], + args, + free=free, + fixed=(fixed, value), + curie=curie, ) - return res.fun + # `cost` is half the sum of squares, which is what `min_func` returns + return res.cost def profile( self, @@ -760,8 +965,8 @@ def profile( args = (k, Phi, sigma_Phi) x0 = np.array([beta, zt, dz, C]) - res = minimize(self.min_func, x0, args=args, bounds=self.bounds) - x_hat, F_min = res.x, res.fun + res = self._fit(x0, args) + x_hat, F_min = res.x, res.cost # every constrained fit is cached, so the root finding below reuses the # scan nodes it lands on rather than paying for them twice @@ -1016,12 +1221,7 @@ def metropolis_hastings( window, xc, yc, taper, process_subgrid, dof_factor, **kwargs ) - lower = np.array( - [-np.inf if b[0] is None else b[0] for b in self.bounds], dtype=float - ) - upper = np.array( - [np.inf if b[1] is None else b[1] for b in self.bounds], dtype=float - ) + lower, upper = self._bound_arrays() def log_posterior(x): if np.any(x < lower) or np.any(x > upper): @@ -1045,11 +1245,8 @@ def step(x, F, scale, chol): # burn-in travelling instead of tuning. Measured on a synthetic, the # posterior mean from a default start sits at a misfit of 121 against # the mode's 50; started here it lands on 50.1. - start = minimize( - self.min_func, - np.array([beta, zt, dz, C], dtype=float), - args=(k, Phi, sigma_Phi), - bounds=self.bounds, + start = self._fit( + np.array([beta, zt, dz, C], dtype=float), (k, Phi, sigma_Phi) ) x = start.x @@ -1207,9 +1404,7 @@ def sensitivity( # simulation from the unresampled solution rather than from the # caller's guess. One extra fit up front, and about a third off the # total for any useful `nsim`. - x0 = minimize( - self.min_func, x0, args=(k, Phi, sigma_Phi), bounds=self.bounds - ).x + x0 = self._fit(x0, (k, Phi, sigma_Phi)).x for sim in range(0, nsim): # a fresh set of prior centres, drawn without disturbing the ones @@ -1220,13 +1415,7 @@ def sensitivity( prior[key] = (loc, self.prior[key][1]) rPhi = rng.normal(Phi, sigma_Phi) - res = minimize( - self.min_func, - x0, - args=(k, rPhi, sigma_Phi, prior), - bounds=self.bounds, - ) - samples[sim] = res.x + samples[sim] = self._fit(x0, (k, rPhi, sigma_Phi), prior=prior).x return list(samples.T) diff --git a/tests/test_bouligand.py b/tests/test_bouligand.py index c97ef58..61d8913 100644 --- a/tests/test_bouligand.py +++ b/tests/test_bouligand.py @@ -13,6 +13,7 @@ import pytest import pycurious +from pycurious.optimise_bouligand import _JACOBIAN_STEP from conftest import synthetic_grid @@ -221,6 +222,214 @@ def test_profile_rejects_unknown_target(bouligand): grid.profile(WINDOW, xc, yc, "curie_depth") +@pytest.mark.parametrize( + "seed, interval", + [ + (1, (11.12, np.inf)), + (2, (1.695, 48.17)), + (3, (5.094, 32.47)), + (4, (24.06, np.inf)), + ], +) +def test_profile_dz_interval_is_unchanged_by_the_optimiser(seed, interval): + """ + Pinned against the values L-BFGS-B produced before `_fit` moved to + `least_squares`, so a change of optimiser cannot quietly move a published + interval. A better inner optimiser is not automatically safe here: it + finds lower constrained minima, which re-anchors the deviance, and that + can move the reported interval. + + The tolerance is 5%, which is loose because the quantity is. An interval + endpoint is where `brentq` crosses the threshold on a deviance curve that + is nearly flat there -- that flatness is the whole reason `dz` needs a + profile rather than a sigma -- so a last-ulp difference in `kv` or in the + FFT behind the synthetic moves it far more than it moves the fit. The same + four intervals come out up to 1.3% different on macOS from Linux with + identical code, which is what set the bound. Pinning tighter tests the + platform's libm, not this package. + + Loose as it is, it still bites: it is what catches a sign-flipped or + dropped Jacobian column, and a `_profiled_misfit` off by a constant + factor. Those move an endpoint by tens of percent, not tenths. + """ + grid, xc, yc = _grid(seed=seed) + grid.reset_priors() + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + _, _, lower, upper = grid.profile(200e3, xc, yc, "dz") + np.testing.assert_allclose([lower, upper], interval, rtol=5e-2) + + +def test_analytic_jacobian_columns_match_the_finite_difference(bouligand): + """ + `_fit` supplies the zt and C columns of the Jacobian in closed form and + differences only beta and dz. If the forward model changes and these are + not re-derived, every fit silently descends a slightly wrong gradient -- + it still converges, just to a worse place, and nothing complains. + """ + from pycurious.optimise_bouligand import _ANALYTIC_COLUMNS, _PARAMETERS + + grid, xc, yc = bouligand + grid.reset_priors() + k, Phi, sigma = grid.window_spectrum(WINDOW, xc, yc, power=2.0) + x = np.array([3.0, 1.0, 20.0, 5.0]) + args = (k, Phi, sigma) + J = grid._jacobian(x, grid.residuals(x, *args), args) + + for name, column in _ANALYTIC_COLUMNS.items(): + index = _PARAMETERS.index(name) + np.testing.assert_allclose( + column(k, sigma), J[: k.size, index], rtol=1e-6, err_msg=name + ) + + +def test_profiled_misfit_is_on_the_same_scale_as_min_func(bouligand): + """ + `_profiled_misfit` returns `res.cost` where it used to return `min_func`, + and every profile interval is differences of those against `F_min`, which + still comes from `min_func`. A constant factor between the two would put + every interval out by that factor while leaving the curve the right shape. + + Holding a parameter at its own fitted value and re-optimising the rest must + recover the unconstrained misfit, which pins the scale. Comparing + `res.cost` to `min_func(res.x)` does not: both are half the sum of squares + of the same vector, so that identity holds however wrong the surrounding + code is. + """ + grid, xc, yc = bouligand + grid.reset_priors() + k, Phi, sigma = grid.window_spectrum(WINDOW, xc, yc, power=2.0) + args = (k, Phi, sigma) + x_hat = grid._fit(np.array([3.0, 1.0, 10.0, 5.0]), args).x + F_min = grid.min_func(x_hat, *args) + + for target, value in (("dz", x_hat[2]), ("beta", x_hat[0]), + ("CPD", x_hat[1] + x_hat[2])): + held = grid._profiled_misfit(target, value, x_hat, args) + assert held == pytest.approx(F_min, rel=1e-6), target + + +def _captured_jacobian(grid, monkeypatch, *fit_args, **fit_kwargs): + """ + Hand back the `(fun, jac, y)` that `_fit` passes to scipy. + + The spy returns nothing rather than delegating: the caller wants the + callables, not the fit, and running one would cost an optimisation whose + result is discarded. + """ + from pycurious import optimise_bouligand as mod + + grabbed = {} + + def spy(fun, y0, jac=None, **kw): + grabbed.update(fun=fun, jac=jac, y0=np.asarray(y0, dtype=float)) + + monkeypatch.setattr(mod, "least_squares", spy) + grid._fit(*fit_args, **fit_kwargs) + return grabbed["fun"], grabbed["jac"], grabbed["y0"] + + +@pytest.mark.parametrize( + "free, fixed, curie", + [ + (None, None, False), + ([0, 1, 3], (2, 25.0), False), # as profile("dz") does + ([0, 1, 3], (2, 30.0), True), # as profile("CPD") does + ], + ids=["unconstrained", "dz-held", "curie"], +) +def test_fit_jacobian_matches_finite_differences(bouligand, monkeypatch, free, + fixed, curie): + """ + The Jacobian `_fit` hands scipy must be the derivative of the residual it + hands scipy alongside it, in every configuration -- including the Curie + one, where `dz = CPD - zt` couples the `zt` and `dz` columns by a chain + rule that no other test reaches. + + Checking the fitted answer instead does not work: trust-region reflective + converges from a wrong Jacobian too, just by a different route. Dropping + the chain-rule term entirely, or flipping its sign, moves the CPD interval + not at all and costs a handful of extra evaluations -- so only the + derivative itself can be tested. + """ + grid, xc, yc = bouligand + grid.reset_priors() + grid.add_prior(beta=(3.0, 0.4), C=(5.0, 1.0)) + try: + k, Phi, sigma = grid.window_spectrum(WINDOW, xc, yc, power=2.0) + args = (k, Phi, sigma) + x_hat = grid._fit(np.array([3.0, 1.0, 10.0, 5.0]), args).x + y0 = x_hat if free is None else np.asarray(x_hat)[free] + + fun, jac, y = _captured_jacobian(grid, monkeypatch, y0, args, free=free, + fixed=fixed, curie=curie) + + analytic = jac(y) + numeric = np.empty_like(analytic) + for i in range(y.size): + h = _JACOBIAN_STEP * max(abs(y[i]), 1.0) + yp, ym = y.copy(), y.copy() + yp[i] += h + ym[i] -= h + numeric[:, i] = (fun(yp) - fun(ym)) / (2.0 * h) + + scale = np.maximum(np.abs(numeric).max(axis=0), 1.0) + np.testing.assert_allclose( + analytic / scale, numeric / scale, atol=1e-6 + ) + finally: + grid.reset_priors() + + +def test_fit_rejects_an_inverted_bound(bouligand): + """ + A collapsed bound (lb == ub) is widened so trust-region reflective has an + interior to work in. An inverted one (lb > ub) is a typo -- `bounds` is + documented as reassignable -- and widening it would silently rewrite it + into whichever number came first, pinning the parameter somewhere the + caller never asked for. + """ + grid, xc, yc = bouligand + grid.reset_priors() + k, Phi, sigma = grid.window_spectrum(WINDOW, xc, yc, power=2.0) + original = list(grid.bounds) + try: + grid.bounds[2] = (30.0, 10.0) + with pytest.raises(ValueError, match="bound"): + grid._fit(np.array([3.0, 1.0, 10.0, 5.0]), (k, Phi, sigma)) + finally: + grid.bounds = original + + +def test_fit_honours_a_held_coordinate(bouligand): + """ + A constrained fit must actually hold what it was told to, including the + Curie-depth case where the held coordinate is dz = CPD - zt and so moves + with a free parameter. + """ + grid, xc, yc = bouligand + grid.reset_priors() + k, Phi, sigma = grid.window_spectrum(WINDOW, xc, yc, power=2.0) + args = (k, Phi, sigma) + x_hat = grid._fit(np.array([3.0, 1.0, 10.0, 5.0]), args).x + + free = [0, 1, 3] + + # dz pinned at 25: the reported cost must be the misfit of the model that + # actually has dz = 25, which is what makes the deviance meaningful + held = grid._fit(x_hat[free], args, free=free, fixed=(2, 25.0)) + beta, zt, C = held.x + assert held.cost == pytest.approx(grid.min_func([beta, zt, 25.0, C], *args)) + assert held.cost >= grid._fit(x_hat, args).cost + + # Curie depth pinned at 30: dz is not fixed, zt + dz is + curie = grid._fit(x_hat[free], args, free=free, fixed=(2, 30.0), curie=True) + beta, zt, C = curie.x + assert curie.cost == pytest.approx( + grid.min_func([beta, zt, 30.0 - zt, C], *args) + ) + + def test_sensitivity_does_not_disturb_priors(bouligand): """ sensitivity used to redraw each prior centre by mutating self.prior in @@ -235,22 +444,26 @@ def test_sensitivity_does_not_disturb_priors(bouligand): grid.sensitivity(300e3, xc, yc, 4, taper=np.hanning, seed=1) assert grid.prior == before - # and still intact when a simulation blows up part way through + # and still intact when a simulation blows up part way through. The count + # is of fits rather than of objective evaluations: the first is the + # unresampled solution the simulations start from, so raising on the third + # lands inside the loop, which is where the prior is copied and where a + # leak would show. calls = {"n": 0} - original = grid.min_func + original = grid._fit def exploding(*args, **kwargs): calls["n"] += 1 - if calls["n"] > 5: + if calls["n"] > 2: raise RuntimeError("boom") return original(*args, **kwargs) - grid.min_func = exploding + grid._fit = exploding try: with pytest.raises(RuntimeError, match="boom"): grid.sensitivity(300e3, xc, yc, 4, taper=np.hanning, seed=1) finally: - del grid.min_func + del grid._fit assert grid.prior == before grid.reset_priors() @@ -415,12 +628,11 @@ def test_metropolis_hastings_is_invariant_to_a_constant_misfit(bouligand): once F grew to a few hundred: exp(-F) underflowed to zero and every proposal was rejected. - Agreement is close but not bit-exact, and the reason is not the acceptance - rule. min_func also drives the search for the mode the chain starts from, - and scipy's convergence tolerance is relative to the value of the - objective, so offsetting it moves where the minimiser stops -- by about - 3e-06 here. Reintroducing exp() would not miss this tolerance narrowly; it - would freeze the chain outright. + The chain starts from a `least_squares` fit of `residuals`, which does not + go through `min_func`, so the offset now moves nothing but the acceptance + ratio. The comparison is left loose regardless -- what it is here to catch + is a chain frozen by a constant, and reintroducing exp() would not miss + this tolerance narrowly but freeze the chain outright. """ grid, xc, yc = bouligand grid.reset_priors()