Skip to content

Share one spectrum across a centroid's fits, in both optimisers #46

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
brmather merged 2 commits into from
Aug 3, 2026
Merged

Share one spectrum across a centroid's fits, in both optimisers #46

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4a321f7
Carry a spectrum's provenance on the spectrum itself
brmather Aug 3, 2026
69e22df
Give Tanaka the spectrum sharing the Bouligand side had
brmather Aug 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 16 additions & 2 deletions CLAUDE.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 # 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
```
Expand Down Expand Up @@ -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
Expand Down
115 changes: 115 additions & 0 deletions pycurious/grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)))

Expand Down Expand Up @@ -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.
Expand Down
123 changes: 15 additions & 108 deletions pycurious/optimise_bouligand.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -215,12 +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`,
# which uses `_spectrum_key` to catch it being handed to a different one
self.last_spectrum = None
self._spectrum_key = None

self.max_processors = kwargs.pop("max_processors", cpu_count())

def _max_thickness(self):
Expand Down Expand Up @@ -361,9 +355,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

Expand Down Expand Up @@ -415,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.
Expand All @@ -433,94 +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.

`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.
"""
if spectrum is None:
spectrum = self._spectrum(
window, xc, yc, taper, process_subgrid, dof_factor, **kwargs
)
key = (window, xc, yc, taper, process_subgrid, dof_factor)
else:
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,
)

self.last_spectrum = spectrum
self._spectrum_key = key
return spectrum

def _bound_arrays(self, free=None):
"""
`self.bounds` as a pair of arrays, with `None` meaning infinite.
Expand Down Expand Up @@ -845,7 +755,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.
"""
Expand Down Expand Up @@ -1047,10 +958,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`

Expand Down Expand Up @@ -1290,9 +1199,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,)
Expand Down Expand Up @@ -1508,9 +1416,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,)
Expand Down
Loading
Loading