From 484896691bd3121c1c40538d343b3fcd1562a185 Mon Sep 17 00:00:00 2001 From: Ben Mather Date: Mon, 27 Jul 2026 10:08:12 +1000 Subject: [PATCH 1/4] Fit with least_squares instead of L-BFGS-B Every Bouligand fit now goes through one private method, `_fit`, which is `scipy.optimize.least_squares` (trust-region reflective) on the existing vector-valued `residuals` rather than `minimize` on the scalar `min_func`. The problem is a sum of squares and `residuals` already exposed the vector, so TRF can use structure L-BFGS-B cannot see. The difference is not marginal. Measured in CPU time -- wall clock is useless on a shared box -- over synthetic windows, L-BFGS-B spent 3.1 s on the four-parameter fit and 16.5 s on a `dz` profile against 13 ms and 49 ms here. Its cost did not track the number of function evaluations at all (34x spread at equal counts), so the time was going into its own machinery, not the forward model. The objective was 1-3% of a profile. Answers are unchanged, which is the reason to do it: - fitted parameters agree to four figures on 8/8 synthetic seeds, and where they differ least_squares finds the lower misfit - the four `dz` intervals pinned in the new test are the values L-BFGS-B produced - over cached EMAG2 spectra in a downstream project, eight vertices at 2000 km reproduce stock exactly, and 0 of 108 profile intervals across three window sizes and three targets excluded their point estimate `_fit` also supplies the two exact Jacobian columns -- `zt` and `C` enter the forward model linearly, so d r/d zt = -2k/sigma and d r/d C = 1/sigma -- and differences only `beta` and `dz`, with a small cache so the Jacobian reuses the residual `least_squares` just evaluated. Both are needed: the columns alone, recomputing the base, gave 1.32x where together they give 1.51x. `beta` enters through the *order* of a Bessel function and has no closed-form derivative; `dz` has one, but it costs an extra `kv` call, exactly what the difference it would replace costs. `_covariance` and `_jacobian` are deliberately left on central differences. They are 0.16 ms, and the calibrated sigmas depend on them. One behavioural difference: `least_squares` refuses `lb == ub`, which L-BFGS-B accepted by pinning the parameter. `_fit` widens a collapsed bound to the narrowest numerically distinct box, which pins it just as effectively and stays close enough to the edge for `_warn_on_bounds` to fire. Two tests drove the fit by monkeypatching `min_func`, which nothing calls now. The sensitivity one patches `_fit` and counts fits instead, so it still blows up inside the simulation loop where the prior is copied; the Metropolis one is unaffected, since `log_posterior` still uses `min_func`, but its docstring explained a discrepancy that no longer exists. Suite: 100 -> 107 tests, 293 s -> 12 s. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01BLGk2RXegPShhDRHcH7FtL --- CLAUDE.md | 26 +++- pycurious/optimise_bouligand.py | 221 +++++++++++++++++++++++++++----- tests/test_bouligand.py | 123 ++++++++++++++++-- 3 files changed, 324 insertions(+), 46 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 27a9070..7e41462 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 # 108 tests, ~25 s +pytest -m "not slow" # 105 tests, ~14 s -- skips the calibration tests that # fit a few hundred realisations pytest tests/test_tanaka.py -q ``` @@ -78,6 +78,17 @@ 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 spends seconds of CPU on this +four-parameter problem — 3.1 s for a fit and 16.5 s for a `dz` profile against +13 ms and 49 ms, measured in CPU time over synthetic windows — and its cost +does not track the number of function evaluations, so the time is going into +its own machinery rather than the forward model. `_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. Benchmark in +`time.process_time`, never wall clock — this is a shared machine. + **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 +213,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..dedc145 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 @@ -65,6 +65,30 @@ # Relative step for the finite-difference Jacobian behind the covariance. _JACOBIAN_STEP = 1.0e-6 +# 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 = { + 1: lambda kh, sigma: -2.0 * kh / sigma, # zt + 3: lambda kh, sigma: 1.0 / sigma, # C +} + # 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 +140,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 +409,154 @@ def _spectrum(self, window, xc, yc, taper, process_subgrid, dof_factor, **kwargs **kwargs ) + 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 the structure that L-BFGS-B + cannot see. The difference is not marginal: measured over eight + synthetic windows, L-BFGS-B spent 3.1 s of CPU on the four-parameter + fit and 16.5 s on a `dz` profile, against 13 ms and 49 ms here, for + the same parameters to four figures and the same intervals. Its + cost did not track the number of function evaluations at all, so + the time was going into its own machinery rather than the forward + model. + """ + 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) + + # map from the free coordinates to the full parameter vector, so the + # chain rule for a constrained fit is one matrix rather than a special + # case in every column + transform = np.zeros((len(_PARAMETERS), len(free))) + for column, index in enumerate(free): + transform[index, column] = 1.0 + if curie: + # dz = CPD - zt, so the dz row differentiates through zt + transform[fixed[0], free.index(1)] = -1.0 + + def expand(y): + x = np.empty(len(_PARAMETERS)) + x[free] = y + if fixed is not None: + # zt is already in place, which the Curie constraint needs + x[fixed[0]] = fixed[1] - x[1] if curie else fixed[1] + return x + + # `least_squares` evaluates the residual at a point and then asks for + # the Jacobian there, so one entry is enough to spare the Jacobian its + # base evaluation. A handful are kept rather than exactly one so that a + # step which is rejected and retried does not thrash it. + cache = {} + + def residuals_at(x): + key = x.tobytes() + r = cache.get(key) + if r is None: + if len(cache) > 3: + cache.clear() + r = cache[key] = self.residuals(x, *args, prior=prior) + return r + + # rows `residuals` appends for the priors, in the order it appends them + prior_rows = [ + index + for index, key in enumerate(_PARAMETERS) + if prior.get(key) is not None + ] + # a column of the full Jacobian is only worth forming if the transform + # carries it through -- the held coordinate of a plain profile does not + needed = [i for i in range(len(_PARAMETERS)) if transform[i].any()] + + def jacobian(y): + x = expand(y) + J = np.zeros((nk + len(prior_rows), len(_PARAMETERS))) + base = None + for index in needed: + analytic = _ANALYTIC_COLUMNS.get(index) + if analytic is not None: + J[:nk, index] = analytic(kh, sigma) + continue + if base is None: + base = residuals_at(x) + step = _JACOBIAN_STEP * max(abs(x[index]), 1.0) + xp = x.copy() + xp[index] += step + # spectral rows only; the prior rows are exact and are filled + # below for every column, so differencing them here would + # count them twice + J[:nk, index] = (residuals_at(xp)[:nk] - base[:nk]) / step + for row, index in enumerate(prior_rows): + J[nk + row, index] = 1.0 / prior[_PARAMETERS[index]][1] + return J @ transform + + lower = np.array( + [ + self.bounds[i][0] if self.bounds[i][0] is not None else -np.inf + for i in free + ] + ) + upper = np.array( + [ + self.bounds[i][1] if self.bounds[i][1] is not None else np.inf + for i in 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. + collapsed = np.isfinite(lower) & np.isfinite(upper) & (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 +692,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 +815,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 +926,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 @@ -1045,11 +1211,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 +1370,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 +1381,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..a4fc035 100644 --- a/tests/test_bouligand.py +++ b/tests/test_bouligand.py @@ -221,6 +221,101 @@ 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. These are the numbers, not a tolerance band: the swap was + adopted on the evidence that it reproduces them. + + A better inner optimiser is not automatically safe here. It finds lower + constrained minima, which re-anchors the deviance, and on a multimodal + window that moves the reported interval -- see the test above. dz was + checked across seeds for exactly that reason. + """ + 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=2e-3) + + +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 + + 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 index in (1, 3): + np.testing.assert_allclose( + _ANALYTIC_COLUMNS[index](k, sigma), J[: k.size, index], rtol=1e-6 + ) + + +def test_fit_cost_is_min_func(bouligand): + """ + `_profiled_misfit` returns `res.cost` where it used to return the value of + `min_func`, and the deviance -- so every profile interval -- is differences + of those. The two are the same quantity only for as long as `residuals` and + `min_func` agree. + """ + grid, xc, yc = bouligand + grid.reset_priors() + k, Phi, sigma = grid.window_spectrum(WINDOW, xc, yc, power=2.0) + args = (k, Phi, sigma) + res = grid._fit(np.array([3.0, 1.0, 10.0, 5.0]), args) + assert res.cost == pytest.approx(grid.min_func(res.x, *args)) + + +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 +330,27 @@ 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 calls["n"] > 2 assert grid.prior == before grid.reset_priors() @@ -415,12 +515,13 @@ 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 moves nothing but the acceptance + ratio and the two chains should now agree to the bit. They are compared + loosely anyway: the point of the test is that the chain is not frozen by a + constant, and pinning it exactly would make the test fail for reasons that + have nothing to do with that. Reintroducing exp() would not miss this + tolerance narrowly; it would freeze the chain outright. """ grid, xc, yc = bouligand grid.reset_priors() From 07daea54794d1a55d26315b0224d0159680f1917 Mon Sep 17 00:00:00 2001 From: Ben Mather Date: Mon, 27 Jul 2026 10:31:04 +1000 Subject: [PATCH 2/4] Test the Jacobian directly, and correct the performance claim Adversarial review found three things. **The performance numbers were an artefact.** The 3.1 s / 16.5 s figures came from an unpinned BLAS: L-BFGS-B calls a threaded one whose workers spin-wait, and `time.process_time` charges every spinning thread, so on a contended machine the same comparison reads 20x to 400x depending on load. Pinned to one thread it is **3x** -- 12.5 ms against 5.0 ms for a fit, 102 ms against 33 ms for a `dz` profile. The durable number, which no threading can distort, is evaluations of the forward model per vertex: 2012 against 359. CLAUDE.md said to prefer `process_time` over wall clock; that is the advice that produced the wrong number, and it now says to pin the BLAS instead. The change is still worth making. It is 3x and 5.6x fewer model evaluations, not two orders of magnitude. **The Curie chain rule had no test.** Deleting `transform[fixed[0], free.index(1)] = -1.0` or flipping its sign passed the entire suite, because trust-region reflective converges from a wrong Jacobian too -- just by a different route, for a few extra evaluations. The old test asserted `curie.cost == min_func(...)`, which tests `expand()` and holds for any Jacobian at all. Replaced with a test that captures the callable `_fit` hands scipy and differences the residual it hands it alongside, in all three configurations and with priors active. Both mutations now fail it. **An inverted bound was silently accepted.** The widening predicate was `upper <= lower`, so `bounds[2] = (30.0, 10.0)` -- a typo, and `bounds` is documented as reassignable -- was rewritten to `[30, 30+3e-11]` and pinned `dz` at 30 with no error, where L-BFGS-B raised. Now `==` only, and an inverted bound reaches `least_squares` and raises. Also replaced `test_fit_cost_is_min_func`, which was a tautology -- both sides are half the sum of squares of the same vector, and it survived a mutation returning `2 * res.cost` from `_profiled_misfit`, the exact failure its docstring claimed to guard. It now checks that profiling a parameter at its own fitted value recovers the unconstrained misfit, which pins the scale. Suite 107 -> 111. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01BLGk2RXegPShhDRHcH7FtL --- CLAUDE.md | 26 ++++--- pycurious/optimise_bouligand.py | 36 +++++++--- tests/test_bouligand.py | 123 ++++++++++++++++++++++++++++---- 3 files changed, 154 insertions(+), 31 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 7e41462..7ef9e51 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -79,15 +79,23 @@ depth. gradients the fits produce. **Every Bouligand fit goes through `_fit`, which is `least_squares`, not -`minimize`.** Do not put L-BFGS-B back. It spends seconds of CPU on this -four-parameter problem — 3.1 s for a fit and 16.5 s for a `dz` profile against -13 ms and 49 ms, measured in CPU time over synthetic windows — and its cost -does not track the number of function evaluations, so the time is going into -its own machinery rather than the forward model. `_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. Benchmark in -`time.process_time`, never wall clock — this is a shared machine. +`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 diff --git a/pycurious/optimise_bouligand.py b/pycurious/optimise_bouligand.py index dedc145..bb03085 100644 --- a/pycurious/optimise_bouligand.py +++ b/pycurious/optimise_bouligand.py @@ -440,14 +440,27 @@ def _fit(self, y0, args, prior=None, free=None, fixed=None, curie=False): 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 the structure that L-BFGS-B - cannot see. The difference is not marginal: measured over eight - synthetic windows, L-BFGS-B spent 3.1 s of CPU on the four-parameter - fit and 16.5 s on a `dz` profile, against 13 ms and 49 ms here, for - the same parameters to four figures and the same intervals. Its - cost did not track the number of function evaluations at all, so - the time was going into its own machinery rather than the forward - model. + 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))) @@ -541,7 +554,12 @@ def jacobian(y): # 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. - collapsed = np.isfinite(lower) & np.isfinite(upper) & (upper <= lower) + # + # 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. + collapsed = np.isfinite(lower) & np.isfinite(upper) & (upper == lower) if collapsed.any(): upper[collapsed] = lower[collapsed] + _DEGENERATE_BOUND * np.maximum( np.abs(lower[collapsed]), 1.0 diff --git a/tests/test_bouligand.py b/tests/test_bouligand.py index a4fc035..00b064c 100644 --- a/tests/test_bouligand.py +++ b/tests/test_bouligand.py @@ -272,19 +272,118 @@ def test_analytic_jacobian_columns_match_the_finite_difference(bouligand): ) -def test_fit_cost_is_min_func(bouligand): +def test_profiled_misfit_is_on_the_same_scale_as_min_func(bouligand): """ - `_profiled_misfit` returns `res.cost` where it used to return the value of - `min_func`, and the deviance -- so every profile interval -- is differences - of those. The two are the same quantity only for as long as `residuals` and - `min_func` agree. + `_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) - res = grid._fit(np.array([3.0, 1.0, 10.0, 5.0]), args) - assert res.cost == pytest.approx(grid.min_func(res.x, *args)) + 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): + """Run `_fit` and hand back the `(fun, jac, y)` it passed to scipy.""" + from pycurious import optimise_bouligand as mod + + grabbed = {} + real = mod.least_squares + + def spy(fun, y0, jac=None, **kw): + grabbed.update(fun=fun, jac=jac, y0=np.asarray(y0, dtype=float)) + return real(fun, y0, jac=jac, **kw) + + monkeypatch.setattr(mod, "least_squares", spy) + grid._fit(*fit_args, **fit_kwargs) + return grabbed + + +@pytest.mark.parametrize( + "free, fixed, curie", + [ + (None, None, False), # unconstrained + ([0, 1, 3], (2, 25.0), False), # dz held, as profile("dz") does + ([0, 1, 3], (2, 30.0), True), # dz = CPD - zt, as profile("CPD") does + ], +) +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] + + got = _captured_jacobian(grid, monkeypatch, y0, args, free=free, + fixed=fixed, curie=curie) + fun, jac, y = got["fun"], got["jac"], got["y0"] + + analytic = jac(y) + numeric = np.empty_like(analytic) + for i in range(y.size): + h = 1.0e-6 * 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): @@ -516,12 +615,10 @@ def test_metropolis_hastings_is_invariant_to_a_constant_misfit(bouligand): proposal was rejected. The chain starts from a `least_squares` fit of `residuals`, which does not - go through `min_func`, so the offset moves nothing but the acceptance - ratio and the two chains should now agree to the bit. They are compared - loosely anyway: the point of the test is that the chain is not frozen by a - constant, and pinning it exactly would make the test fail for reasons that - have nothing to do with that. Reintroducing exp() would not miss this - tolerance narrowly; it would freeze the chain outright. + 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() From aef207dc633c7af46219d9aca1d0e2cd93c8302b Mon Sep 17 00:00:00 2001 From: Ben Mather Date: Mon, 27 Jul 2026 11:15:33 +1000 Subject: [PATCH 3/4] Loosen the pinned dz intervals to a tolerance the quantity supports The four pinned intervals failed on macOS at rtol=2e-3, by up to 1.3%, with identical code -- the only CI leg that failed. An interval endpoint is where `brentq` crosses the threshold on a deviance curve that is nearly flat there, which is the whole reason dz needs a profile rather than a sigma. A last-ulp difference in `kv`, or in the FFT behind the synthetic, moves the crossing far more than it moves the fit. Pinning to 0.2% was asserting on the platform's libm, not on this package. 5%, which is ~4x the observed cross-platform spread. Still load-bearing: a sign-flipped zt column, zeroed finite-difference columns, and a `_profiled_misfit` off by a factor of two each fail all four cases at the new tolerance. Those move an endpoint by tens of percent. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01BLGk2RXegPShhDRHcH7FtL --- tests/test_bouligand.py | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/tests/test_bouligand.py b/tests/test_bouligand.py index 00b064c..260f492 100644 --- a/tests/test_bouligand.py +++ b/tests/test_bouligand.py @@ -234,20 +234,29 @@ 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. These are the numbers, not a tolerance band: the swap was - adopted on the evidence that it reproduces them. - - A better inner optimiser is not automatically safe here. It finds lower - constrained minima, which re-anchors the deviance, and on a multimodal - window that moves the reported interval -- see the test above. dz was - checked across seeds for exactly that reason. + 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=2e-3) + np.testing.assert_allclose([lower, upper], interval, rtol=5e-2) def test_analytic_jacobian_columns_match_the_finite_difference(bouligand): From 3bbed8ee51e26e15bb87c925cbad2d7a0e37d7b7 Mon Sep 17 00:00:00 2001 From: Ben Mather Date: Mon, 27 Jul 2026 11:32:28 +1000 Subject: [PATCH 4/4] Fold the Jacobian's constant part out, and stop stating things twice Four parallel cleanup reviews. Applied: **The affine map had two implementations.** `transform` said the `dz` row differentiates through `zt`; `expand` separately hand-coded `x[fixed[0]] = value - x[1]`. Same relation, two representations, kept in step by review -- and an inconsistency between them is close to invisible, because the fit converges from a wrong Jacobian anyway, just by a longer route. `expand` is now `transform @ y + offset`, so the Jacobian is the derivative of the expansion by construction. (The matmul stays: measured at 0.98 us against 1.55 us for column indexing, and indexing cannot express the chain-rule term.) **The Jacobian rebuilt its constant part on every call.** The analytic columns depend only on `kh` and `sigma`, and the prior rows only on the prior -- neither on `x`. Hoisting them out drops the lazy `base = None` guard and a four-line comment about not double-counting prior rows, because the differenced columns no longer touch those rows at all. **`needed` probed the matrix it had just built** -- four `ndarray.any()` calls on 3-element rows, measured at 4.4-6.2 us per fit, over half the per-call setup. It is `free` minus the closed-form columns, plus the held one under the Curie constraint, which is also what it means. **The residual cache kept four entries for a one-deep pattern.** Instrumented over six seeds: `residuals` calls are exactly `nfev + 2*njev`, so the base lookup hits every time and nothing else ever does. The comment justified the extra entries by rejected steps being retried, which trust-region reflective does not do -- a shrunken region gives a different point. One entry, same hit rate. (The cache itself earns its keep: ~0.1 ms of hashing against 9-27 ms of forward model per profile.) **`_bound_arrays`** -- `_fit` and `metropolis_hastings` had each written the `None`-means-infinite unpacking, and had already drifted on dtype. **`_ANALYTIC_COLUMNS` keyed by name** rather than by position, so `_PARAMETERS` order is encoded once. Same for `free.index(_ZT)`. Tests: the Jacobian spy no longer runs a fit it discards; the FD reference step uses `_JACOBIAN_STEP` instead of a literal, which would otherwise let a retuned constant pass a test that no longer measures the scheme; parametrised ids; and `assert calls["n"] > 2` is gone, since `exploding` raises on exactly that condition and the assertion could not fail. Mutation-checked after the refactor: chain-rule dropped and sign-flipped, `zt` column sign-flipped, prior rows zeroed, `_profiled_misfit` doubled, and the new `differenced` list missing its Curie entry all still fail. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01BLGk2RXegPShhDRHcH7FtL --- CLAUDE.md | 4 +- pycurious/optimise_bouligand.py | 142 ++++++++++++++++++-------------- tests/test_bouligand.py | 35 ++++---- 3 files changed, 101 insertions(+), 80 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 7ef9e51..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 # 108 tests, ~25 s -pytest -m "not slow" # 105 tests, ~14 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 ``` diff --git a/pycurious/optimise_bouligand.py b/pycurious/optimise_bouligand.py index bb03085..afd017e 100644 --- a/pycurious/optimise_bouligand.py +++ b/pycurious/optimise_bouligand.py @@ -62,9 +62,13 @@ # 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 @@ -85,8 +89,8 @@ # `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 = { - 1: lambda kh, sigma: -2.0 * kh / sigma, # zt - 3: lambda kh, sigma: 1.0 / sigma, # C + "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 @@ -409,6 +413,28 @@ 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. @@ -470,83 +496,76 @@ def _fit(self, y0, args, prior=None, free=None, fixed=None, curie=False): kh, sigma = args[0], args[2] nk = np.size(kh) - # map from the free coordinates to the full parameter vector, so the - # chain rule for a constrained fit is one matrix rather than a special - # case in every column + # 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 - if curie: - # dz = CPD - zt, so the dz row differentiates through zt - transform[fixed[0], free.index(1)] = -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): - x = np.empty(len(_PARAMETERS)) - x[free] = y - if fixed is not None: - # zt is already in place, which the Curie constraint needs - x[fixed[0]] = fixed[1] - x[1] if curie else fixed[1] - return x + return transform @ y + offset # `least_squares` evaluates the residual at a point and then asks for - # the Jacobian there, so one entry is enough to spare the Jacobian its - # base evaluation. A handful are kept rather than exactly one so that a - # step which is rejected and retried does not thrash it. - cache = {} + # 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() - r = cache.get(key) - if r is None: - if len(cache) > 3: - cache.clear() - r = cache[key] = self.residuals(x, *args, prior=prior) - return r - - # rows `residuals` appends for the priors, in the order it appends them + 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 ] - # a column of the full Jacobian is only worth forming if the transform - # carries it through -- the held coordinate of a plain profile does not - needed = [i for i in range(len(_PARAMETERS)) if transform[i].any()] + 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 = np.zeros((nk + len(prior_rows), len(_PARAMETERS))) - base = None - for index in needed: - analytic = _ANALYTIC_COLUMNS.get(index) - if analytic is not None: - J[:nk, index] = analytic(kh, sigma) - continue - if base is None: - base = residuals_at(x) + 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 exact and are filled - # below for every column, so differencing them here would - # count them twice + # 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 - for row, index in enumerate(prior_rows): - J[nk + row, index] = 1.0 / prior[_PARAMETERS[index]][1] return J @ transform - lower = np.array( - [ - self.bounds[i][0] if self.bounds[i][0] is not None else -np.inf - for i in free - ] - ) - upper = np.array( - [ - self.bounds[i][1] if self.bounds[i][1] is not None else np.inf - for i in free - ] - ) + 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 @@ -559,7 +578,9 @@ def jacobian(y): # `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. - collapsed = np.isfinite(lower) & np.isfinite(upper) & (upper == lower) + # 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 @@ -1200,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): diff --git a/tests/test_bouligand.py b/tests/test_bouligand.py index 260f492..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 @@ -266,7 +267,7 @@ def test_analytic_jacobian_columns_match_the_finite_difference(bouligand): 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 + from pycurious.optimise_bouligand import _ANALYTIC_COLUMNS, _PARAMETERS grid, xc, yc = bouligand grid.reset_priors() @@ -275,9 +276,10 @@ def test_analytic_jacobian_columns_match_the_finite_difference(bouligand): args = (k, Phi, sigma) J = grid._jacobian(x, grid.residuals(x, *args), args) - for index in (1, 3): + for name, column in _ANALYTIC_COLUMNS.items(): + index = _PARAMETERS.index(name) np.testing.assert_allclose( - _ANALYTIC_COLUMNS[index](k, sigma), J[: k.size, index], rtol=1e-6 + column(k, sigma), J[: k.size, index], rtol=1e-6, err_msg=name ) @@ -308,28 +310,33 @@ def test_profiled_misfit_is_on_the_same_scale_as_min_func(bouligand): def _captured_jacobian(grid, monkeypatch, *fit_args, **fit_kwargs): - """Run `_fit` and hand back the `(fun, jac, y)` it passed to scipy.""" + """ + 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 = {} - real = mod.least_squares def spy(fun, y0, jac=None, **kw): grabbed.update(fun=fun, jac=jac, y0=np.asarray(y0, dtype=float)) - return real(fun, y0, jac=jac, **kw) monkeypatch.setattr(mod, "least_squares", spy) grid._fit(*fit_args, **fit_kwargs) - return grabbed + return grabbed["fun"], grabbed["jac"], grabbed["y0"] @pytest.mark.parametrize( "free, fixed, curie", [ - (None, None, False), # unconstrained - ([0, 1, 3], (2, 25.0), False), # dz held, as profile("dz") does - ([0, 1, 3], (2, 30.0), True), # dz = CPD - zt, as profile("CPD") does + (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): @@ -354,14 +361,13 @@ def test_fit_jacobian_matches_finite_differences(bouligand, monkeypatch, free, 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] - got = _captured_jacobian(grid, monkeypatch, y0, args, free=free, - fixed=fixed, curie=curie) - fun, jac, y = got["fun"], got["jac"], got["y0"] + 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 = 1.0e-6 * max(abs(y[i]), 1.0) + h = _JACOBIAN_STEP * max(abs(y[i]), 1.0) yp, ym = y.copy(), y.copy() yp[i] += h ym[i] -= h @@ -458,7 +464,6 @@ def exploding(*args, **kwargs): grid.sensitivity(300e3, xc, yc, 4, taper=np.hanning, seed=1) finally: del grid._fit - assert calls["n"] > 2 assert grid.prior == before grid.reset_priors()