Publication-ready v2: GitHub Actions CI, Sphinx docs, PyPI publishing; drop Docker/Binder - #44
Conversation
Added/extended three forms of Tanaka method, with flexibility for with/without errors.
Code to generate synthetic magnetic anomaly data for Curie point depth analysis/comparison.
Collection of scripts from RDelhaye. Included are mag_grid2.npz (compressed Tellus aeromag data as of 18/05/2018), and tst_data.npz, a synthetic data set from the Bouligand_forward.py code.
- beta = 3 - dz = 10
Left room for Rob to add Tanaka optimisation strategy in Ex3-Optimisation-routine.ipynb
Ex2a essentially complete. Ex4 in progress.
…Testing Done to integrate Tanaka additions
... and a bunch of other modifications
Now to merge with master branch
Merge Testing to master
Fix Tanaka to fit the amplitude spectrum
Move project metadata from setup.py into pyproject.toml, leaving setup.py as a compatibility shim. Requires Python 3.9+. The Cython radon extension (src/*) is removed along with its build machinery; the azimuthal_spectrum method that depended on it is already gone. Optional dependencies are now grouped into download / mapping / examples / test extras rather than being all-or-nothing. MANIFEST.in and .gitignore are updated for the Examples/ move out of the package directory. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
Follow through on splitting CurieOptimise into CurieOptimiseBouligand and CurieOptimiseTanaka: update the exports, the tests and the package docstring. Drop the Cython import test and the azimuthal_spectrum call. install_documentation moves off the removed pkg_resources/distutils APIs onto importlib.resources and shutil, and now looks for Examples/ both inside the package and at the repository root so editable installs work. mapping.export_netcdf4 uses np.float64 rather than the np.float alias removed in numpy 2. The tests/test_mag_data.txt symlink is repointed at the relocated Examples/data. Note: tests/test_grid.py::test_Tanaka currently fails. ComputeTanaka lost its abs() in an earlier commit on this branch, so it returns -10.34. That is addressed in the Tanaka work that follows, not here. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
Document the three-class API, the Python 3.9+ requirement and the new optional-dependency extras. Notebook data paths follow Examples/ moving out of the package directory. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
Three ways the parallel path could hang indefinitely, all fixed: An exception inside func killed the worker without ever posting a result, so the parent blocked forever on q_out.get(). Failures are now captured and ferried back; on_error='raise' (default) reports which centroid failed, on_error='ignore' fills it with NaN and warns. A worker could also die before reaching that try/except -- most often an argument that pickles in the parent but cannot be reconstructed in the child, such as a taper or process_subgrid defined in a notebook cell. Such arguments are now detected up front and the routine falls back to serial with a warning. As a general backstop for any other cause of worker death, results are collected with a timeout and the workers are checked for liveness rather than waited on unconditionally. Feeding q_in could block the parent once the workers had stopped draining it, so it is no longer bounded at a single item. Also: seed= gives each centroid an independent, reproducible stream via SeedSequence, so stochastic routines no longer depend on the number of processors. Under fork every worker previously inherited the same RNG state and drew an identical sequence, which would put spatially coherent artefacts into a sensitivity map. An empty centroid list now raises instead of NameError, and the output shape is taken from the first successful result rather than whichever arrived last. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
_taper_spectrum used dk = 2*pi/((N-1)*dx), but the DFT fundamental is 2*pi/(N*dx). Every wavenumber was overstated by N/(N-1) and so every depth understated by (N-1)/N -- 0.5% at N=201, worse for small windows. This affected Bouligand and Tanaka alike. An injection test now recovers a known depth to 0.02 km, and the radial wavenumber grid is asserted equal to numpy's own fftfreq grid. _FFT_spectrum indexed the zero frequency as (nr-1)//2, which is the fftshift DC index only for odd nr. subgrid() always returns an odd number of points so the optimisers were insulated, but radial_spectrum is public and mis-centred the wavenumber grid for an even-sided array. It also assumed a square input while _taper_spectrum only warned about non-square. Both fixed. Bin edges were doubly inclusive, so cells landing exactly on a boundary were counted in two annuli; bins are now half-open. _FFT_spectrum additionally returns the per-bin cell count, surfaced through radial_spectrum(return_counts=True) for callers that need the standard error of the binned mean rather than the scatter of the cells. The default return is unchanged at three values. test_optimise.py::test_optimisation tolerances are widened. That test was already marginal -- zt sat at 0.091 of a 0.1 tolerance -- and the corrected wavenumbers pushed it over. The underlying cause is that the four-parameter Bouligand fit is poorly conditioned on this synthetic: sweeping window size and centroid gives dz between 6.4 and 11.9 km against a truth of 10.0, a spread far larger than the 0.33% wavenumber change. Recorded for the Bouligand workstream rather than papered over. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
Collected while finishing the Tanaka class: the four-parameter fit is poorly conditioned, min_func ignores sigma_Phi and so fits unweighted, sensitivity perturbs by the per-cell scatter rather than the standard error of the binned mean, calculate_CPD is a stub whose signature is incompatible with the Tanaka sibling, and Bouligand/Ex3 imports a scipy function that no longer exists. Also notes what v2-core changed underneath it: corrected wavenumbers, and the new parallel semantics. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
tests/synthetic.py generates a magnetic anomaly by filtering white noise with the square root of bouligand2009, so its expected radial power spectrum is the analytic model for any beta, zt and dz. Unlike Bouligand_forward.py this is fast enough to run in the test suite, and unlike test_mag_data.txt the true parameters can be varied. tests/test_recovery.py uses it to assert that CurieOptimiseBouligand recovers the parameters it was given, and that the generator really does produce the spectrum it claims -- without which the recovery tests would be circular. This also corrects the diagnosis recorded in docs/bouligand-findings.md. The four-parameter fit is not poorly conditioned: given 1024 km of data it recovers beta to +-0.1 and zt to +-0.07 in 0.1s. The instability seen earlier is a property of the legacy 305 km fixture, which is only ~30x the source thickness and has too few low-wavenumber bins to constrain the long wavelengths. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
Both spectral bands are fitted with curve_fit weighted by the measured scatter, so the covariance propagates into an uncertainty on the Curie depth. Depths are now returned positive downwards; previously optimise handed back raw gradients, so callers printed a "top of magnetic source" of -0.95 km. Two defects in the weighting: The centroid branch computed sigma as log(exp(sigma_Phi)/k), subtracting ln(k) from a standard deviation. Since ln(k) is deterministic the uncertainty is unchanged, and the erroneous term inflated sigma from ~0.6 to ~4.0 at the lowest wavenumbers -- precisely the bins the centroid fit depends on. Phi is the mean over each annulus, so the fit needs the standard error, not the scatter of the cells. That is not sigma/sqrt(N) either: the cells are not independent. Hermitian symmetry makes about half redundant, and tapering correlates neighbours. Monte Carlo over 150-250 realisations at n = 128, 256 and 512 puts N/N_eff at 2.0 with no taper -- exactly the Hermitian factor -- and 3.3 for hanning. Reported uncertainties were otherwise about 2x too small. Bands are now required arguments in rad/km. There are no defaults because the historically used ones are only accurate by coincidence: on an exact noiseless layer the old centroid band reaches |k|d = 3.1, biasing z0 low by 2.2 km, which happened to be cancelled by an opposing fractal bias. check_bands() reports point counts, wavelengths and |k|d so a choice can be checked rather than assumed. Since the previous units were cycles/km, and such a call still returns a plausible number rather than raising, a guard warns when a band sits implausibly low in the spectrum. Added: an optional beta argument removing the fractal contribution, which otherwise biases zt high by (beta-1)/2*kbar; and sensitivity(), which samples band placement as well as spectral noise, since band choice usually dominates. The /(2*pi) rescaling is removed. It cancelled in the gradient and in the gradient variance, but not in the returned intercepts, which came back 6.28x too small -- and the notebooks plot with those. tanaka1999 and ComputeTanaka are deprecated with FutureWarning rather than DeprecationWarning, which the default filter hides from exactly the Jupyter audience concerned. ComputeTanaka regains the abs() it lost earlier on this branch, which was making test_Tanaka return -10.34. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
It printed its diagnostics and also raised them as warnings, so each problem appeared twice in a notebook. It now prints when verbose and warns otherwise, the latter being what a script wants. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
Ex2 is largely rewritten. Its derivation of equation (4) had both
exponents as e^{-|k|(-d)}, making the bracket identically zero where it
should be 2*sinh(|k|d); the sinh step is now shown, along with the fact
that approximating it as 2|k|d is what requires |k|d << 1.
More importantly, the notebook no longer presents its own result as a
success. Working through it: the historical bands give a Curie depth of
11.4 km against a true 10.3 km, which looks good until check_bands warns
that the centroid band reaches |k|d = 3. Correcting for the fractal
magnetisation then recovers Z_t almost exactly (0.29 km against 0.305)
while the Curie depth collapses to 5.2 km -- because the two biases had
been cancelling. Narrowing the centroid band to where the approximation
holds leaves only two spectral estimates, so this 200 km window cannot
constrain the centroid depth at all. That is now the notebook's
conclusion rather than a footnote.
Ex3's sweeps are converted to rad/km and record what they demonstrate:
Curie depth varies from -0.8 to 11.5 km across band positions while the
reported uncertainty never exceeds 0.2 km. Failed fits are collected as
NaN rather than aborting the sweep.
Ex4 and Ex5 pass bands explicitly, since there are no defaults, and Ex4
gains a sensitivity section. Ex1 gains a note on why Tanaka needs power=1
where Bouligand uses power=2. Plots use rad/km throughout, dropping the
/(2*pi) scaling whose axis labels did not describe what was plotted, and
draw fitted lines from the returned intercepts rather than re-deriving
them.
README and 0-StartHere linked to Ex1-Plot-power-spectrum, which is the
Bouligand filename, and listed only three Tanaka notebooks when five
exist. The README paths also still pointed inside the package.
Outputs are regenerated by running each notebook headlessly under
nbconvert. Ex5 is excluded: it downloads EMAG2 and needs network access.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
Examples/data is a download cache: Ex5 fetches EMAG2 v3 (497 MB) and the Li et al. (2017) reference model (124 MB) through pycurious.download, both reproducible from their URL and checksum. Everything there is now ignored except the small synthetic anomaly the tests depend on. Also ignore the GeoTIFFs and npz grids the notebooks write as output. Personal working files -- collaborator data, scratch notebooks -- are excluded via .git/info/exclude instead, so they stay local rather than being imposed on other contributors. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
Ex5 referenced curie_depth_sphere when the variable is CPD_sphere, so its GeoTIFF export had never run; the notebook cannot have been executed to completion since the variable was renamed. A Curie depth colourbar was also labelled nT. mapping.export_geotiff and import_geotiff import osgeo, but GDAL was not declared anywhere in the project metadata. It is now a `geotiff` extra rather than part of `mapping`, since the bindings need a matching libgdal already present and would otherwise make pip install pycurious[mapping] fail for most users. Both import sites now explain what is missing rather than raising a bare ModuleNotFoundError. Verified by running Ex5 against the real EMAG2 v3 grid: 12 of its 14 code cells execute, covering the whole analysis -- cached data load, gridding to Irish Transverse Mercator, 780 centroids through optimise_routine with the rad/km bands, the cartopy maps and the Li et al. (2017) comparison. The two that remain are the GeoTIFF export and re-import, which need the GDAL bindings this machine cannot build. Outputs are therefore not stored for Ex5, unlike Ex1 to Ex4. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
optimise was solving a weighted nonlinear least-squares problem with a
general-purpose scalar minimiser and returning res.x, discarding the
residual structure that yields the covariance. It now reports a sigma per
parameter, and profile() gives an interval that does not assume symmetry.
Four defects had to be fixed first, since a covariance is only meaningful
once the misfit is a proper chi-squared.
min_func ignored sigma_Phi and passed a literal 1.0, so every bin carried
equal weight. Note the fix is *not* the one docs/bouligand-findings.md
prescribes: weighting by the raw within-annulus scatter measures worse than
not weighting at all (beta error 0.119 against 0.073 unweighted, and it
breaks test_optimisation), because sigma_Phi is nearly flat across the
spectrum at the theoretical 1.28 for a complex-Gaussian FFT and so carries
no information. The weight is the uncertainty of the binned mean. That
lands chi2_red at 1.13.
The degrees-of-freedom deflation was a single number per taper, which
overstated the information in the innermost bins -- the ones a centroid
depth leans on hardest. Monte Carlo over 400 realisations at n = 128, 256
and 512 shows a fixed number of cells is lost to correlation however few
the annulus holds:
dof(N) = dof_inf * N / max(N - lost, 1)
The asymptotes come back at 2.01-2.03, 3.22-3.30 and 3.03-3.11, confirming
the published table; what was missing is lost ~ 3.4, 4.9, 4.8. On an unseen
grid size the reported sigma of the lowest bin goes from 0.65x the true
scatter to 1.05x. This moves Tanaka's centroid depth, which is the point:
on the legacy fixture z0 6.5484 +/- 0.1262 -> 6.4112 +/- 0.1378 and CPD
12.150 -> 11.876. zt shifts by 0.0001, its band having no sparse bins.
Neighbouring radial bins are not independent either -- a taper spreads each
wavenumber over a main lobe several bins wide -- and nothing corrected for
that. Measured over 200 independent realisations, treating them as
independent understates every uncertainty by about 30%. The covariance is
now generalised least squares over a banded correlation estimated from the
residuals, which brings the ratio of true spread to reported sigma to
0.997, 0.989 and 0.991 for beta, zt and C.
dz keeps a ratio near 1.44 after that correction, and no symmetric interval
can fix it: its likelihood has a long upper tail (Mather & Fullea, 2019),
and a single realisation can land 75% high at a *lower* misfit than the
truth. Hence profile(), which traces the deviance with one parameter held
and the rest re-optimised, and reports where it crosses chi2_1. It also
profiles the Curie depth directly, by substituting dz = CPD - zt, rather
than propagating a symmetric sigma_dz.
Both the covariance and the profile describe the scatter of the spectrum at
a fixed window and model. They do not cover the choice of either, which on
a small grid is the larger term; the docstrings say so.
Also here:
- 1e99 as the overflow sentinel was a trap: its ULP is ~2e82, so adding a
prior term to it was exactly a no-op. Non-finite residuals are now
replaced elementwise, so the misfit still falls as the model returns to
where it can be evaluated.
- sensitivity redrew prior centres by mutating self.prior in place and
restoring it afterwards, leaving the instance corrupted if anything
raised in between. min_func takes a prior argument instead, and prior
values are tuples so it cannot happen again.
- seed= reached np.hanning and raised, so parallelise_routine's
per-centroid seeding was unreachable. It is an explicit parameter now,
and _taper_spectrum rejects unknown keywords rather than swallowing them
when taper=None, which is how it failed silently on the notebooks that
pass no taper.
optimise returns eight values rather than four, matching the shape of the
Tanaka sibling. The covariance matrix is available via return_cov, but
deliberately not from optimise_routine: parallel._collect dispatches on
result dimensionality and a 4x4 per centroid does not fit.
test_optimise.py::test_optimisation becomes a smoke test. Its fixture is
305 km across for a 10 km layer, so dz moves over 6.4-11.9 km across
overlapping centroids; a tolerance tight enough to mean anything would
break on any legitimate change. Recovery is asserted against generated
synthetics instead, and there the per-seed dz tolerance is replaced by two
better-founded claims: that the profile interval contains the truth, and
that the mean over eight seeds is within 10% (measured +6.6% and +6.9%,
both high, because the tail pulls it).
The calibration test is the one that matters and is marked slow. sensitivity
resamples each bin independently and metropolis_hastings shares the same
diagonal weighting, so both agree with the covariance to within 10% whether
or not it is right. Only an ensemble over independent realisations of the
field is an outside check.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PzxnSJBbkFRdHEF32gHvbo
Second half of the uncertainty work. The first commit gave
CurieOptimiseBouligand a covariance and a profile; this one makes the
synthetic generator available to callers, makes the MCMC sample, and
carries the correlation correction across to the Tanaka side.
fractal_anomaly moves from tests/ into the package, so a notebook or a
user can generate a field with a Curie depth they chose rather than
loading the one 305 km fixture. Its noise is now normalised by n. C
absorbs every constant between the noise and the spectrum, and without
that normalisation it also absorbed 2*ln(n), so the recovered C moved
with the size of the grid for no physical reason. What remains is fixed
and explicable -- the log-averaging offset of gamma, and the power a
taper removes, together about -2.5 under np.hanning and -0.6 with none --
and the docstring says so. C is a nuisance parameter, not something with
a true value to check against.
metropolis_hastings did not sample. Comparing exp(-F) directly underflows
to zero once F reaches a few hundred, which it does for any real
spectrum, so every proposal was rejected: 9 distinct states in 2000
draws at an acceptance rate of 0.004. Acceptance is decided in log space
now, the chain respects self.bounds, and it is seeded.
Two things mattered more than the acceptance arithmetic. The chain starts
at the mode rather than the caller's guess, because a random walk spends
its whole burn-in travelling there -- from a default start the posterior
mean sat at a misfit of 121 against the mode's 50. And the proposal is
drawn along the fit covariance rather than a diagonal: beta and zt
correlate at -0.92 and the marginal widths differ by thirty times, so no
per-parameter width can move along the ridge they lie on. Acceptance goes
from 0.004 to 0.23, and the posterior agrees with the sensitivity
ensemble and the profile interval.
Adapting the covariance from the burn-in instead, as Haario prescribes,
does not work here: a chain started at the mode with too small a step
learns a covariance narrower than the truth, proposes from it and
confirms itself. Measured that way it reported a sigma on dz of 0.8
against a true 8.7. The analytic covariance is already the right shape,
so the burn-in only tunes a scalar multiplier on it.
Tempering is implemented properly rather than left inert -- the previous
burn-in accepted only strict improvements, which is hill-climbing, so the
documented tempering did nothing -- but it is off by default, because
turning it on made every case worse. It existed to work around large
regions of zero probability, which were the underflow rather than a
property of the problem, and annealing fights the scale tuning: high
temperature makes everything acceptable, driving the scale up, and the
scale then collapses as the temperature falls and freezes the chain.
On the Tanaka side, _fit_band gets the same correlation correction, since
that is a property of the taper and not of the method. The gradients are
untouched. Against 80 independent synthetics the ratio of true spread to
reported sigma_zt improves from 1.48 to 1.12. sigma_z0 improves from 1.54
to 1.34 and stays understated, which no covariance can fix: the centroid
gradient is fitted over a handful of the longest wavelengths a window
resolves and its distribution is heavy-tailed, spanning 4.2 to 25.3 km
for a true 11.0. The docstring says to read it as a lower bound.
check_bands now quotes a bias in km instead of leaving the user to judge
whether |k|d is small enough. Measured against exact layer spectra,
Z_b bias = -0.17 * thickness * |k|d
to within 6% over thicknesses of 10-40 km. Both thresholds were wrong in
opposite directions: |k|d > 1 only fired once the Curie depth was already
3.3 km out, having said nothing at |k|d = 0.5 where it is 1.7 km out, and
the zt rule warned at a wavelength where the measured error is 0.05% of
the source thickness -- complaining about ten metres while silent about
kilometres. Now 0.3 and four times the thickness respectively.
The beta correction is documented as verified for zt only. Removing the
ln|k| term does not make the two models agree, because bouligand2009
carries beta inside its cosh/Bessel factor too, and the remainder lands
on the centroid. Fitting an exact spectrum with zt=1 and dz=20, the error
in Zb is +16.6 km at beta=1, +4.2 at 2, +0.6 at 3 and -0.5 at 4, while zt
comes back to three decimal places throughout. test_tanaka_recovers_curie_depth
passes because beta=3 is where that residual cancels the opposing |k|d
bias -- a coincidence of one value, now pinned by a test parametrised
over beta so it reads as one.
No profile() for Tanaka. Each band is a straight-line fit, so its misfit
is exactly quadratic in the gradient and the profile interval is provably
the covariance interval -- a deviance of 3.841459 against a threshold of
3.841459. It would return the number _fit_band already returns. Recorded
in the docstring so the omission reads as a conclusion.
docs/bouligand-findings.md is annotated rather than rewritten. Two of its
items were wrong on measurement and that is worth keeping: item 2's
prescribed weighting measures worse than the bug it replaces, and item
4's table misreads the Tanaka signature it compares against. None of the
seven anticipated the between-bin correlation, which was the largest of
them.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PzxnSJBbkFRdHEF32gHvbo
Ex1-Ex4 now generate their own anomaly with pycurious.fractal_anomaly instead of loading the 305 km fixture, so every number printed can be checked against a Curie depth we chose. Executed headlessly; Ex5 is left unexecuted as its Tanaka counterpart is, since it needs half a gigabyte of EMAG2. Ex1 keeps the analytic sweeps and gains a reason for each: zt sets the short-wavelength slope, dz acts only at long wavelengths and saturates, beta tilts the whole curve. That is why dz is the hard one. The measured-against-analytic comparison uses the parameters the field was generated from rather than a hard-coded guess, and the taper section now says why to taper: over 120 realisations an untapered fit returns beta = 3.057 and zt = 0.945 against truths of 3.0 and 1.0, where hanning gives 3.003 and 1.001. Ex2 is the centrepiece. It reports every parameter with its uncertainty against the truth, and states no truth for C, which is a level rather than a depth and absorbs constants the fit cannot recover. It shows the fitted spectrum lying almost on top of the true one -- a 17% error in dz is nearly invisible over most of the range, which is the difficulty of the method in one figure. The window sweep is the practical lesson: at 300 km, dz comes back at 46 km against a true 20, and reports +/- 67. The fit is not quietly wrong, it says so. Ex3 is the notebook this work existed to make possible. Four ways to characterise the uncertainty, on one dataset, in increasing cost: covariance, profile, spectral resampling, MCMC. Three agree closely and all four contain the truth. The point of the notebook is that the agreement of the middle two proves nothing: sensitivity resamples each bin independently and metropolis_hastings sums over bins as though they were independent, so they share an assumption and will agree whether or not it holds. It does not hold -- adjacent bins correlate at 0.36 under hanning -- which is why the corrected covariance reports a sigma_beta about 40% larger. Only an ensemble over independent realisations of the field can say which is right, and it says the covariance is. Ex4 maps a synthetic whose Curie depth is 21 km everywhere. The map runs from 5 to 35 km, with structure that would read as geology. That is the noise floor a real map has to beat before any of it means anything. It also shows the scatter across windows exceeding the average reported sigma, and one window whose 95% interval excludes the truth, which is what a 95% interval does one time in twenty and worth seeing once. Ex5 gets the eight-value return, calculate_CPD, an uncertainty map beside the depth map, and a colourbar that no longer labels a Curie depth in nT. The Tanaka notebooks are re-executed rather than merely re-read. Their numbers moved: the per-bin degrees-of-freedom correction and the correlation correction shift the centroid depth and widen its error bar, taking the Curie depth from 11.39 +/- 0.39 to 11.00 +/- 0.53 km. Ex2's prose is updated to match, and check_bands now quantifies the |k|d bias it warns about -- about 5 km on that band -- rather than saying "several". Also the duplicated "Varying window size" heading in Ex3, and the sentence in Ex5 crediting the Bouligand method in a Tanaka notebook. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PzxnSJBbkFRdHEF32gHvbo
All eleven notebooks execute end to end against the current code, with no errors. Their outputs are then stripped, so what lands here is the source only: outputs, execution counts, and the wall-clock timestamps nbconvert writes into cell metadata. That keeps the diffs readable and stops re-running a notebook from showing up as a change to it. Both Ex5 notebooks now run in full, which needed GDAL for the GeoTIFF round-trip at the end. The Tanaka one has never completed before: it referenced a variable that had been renamed, so the export cell had always raised, and the fix in 2bf1869 was made without being able to run past it. Fourteen of fourteen cells now. Stray empty code cells removed from five notebooks. The one on the landing page is deliberate and stays. Note the notebooks discuss their own results in prose -- "falling from 11.0 km to about 5 km" and the like -- so with the outputs stripped those numbers cannot be checked by reading the file on its own. They were checked against the executed outputs before stripping, and re-running any notebook reproduces them; every synthetic is seeded. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PzxnSJBbkFRdHEF32gHvbo
Cleanup pass over the uncertainty work. No behaviour is intended to change; the parameters and intervals it produces are identical to 1e-14, and the notebooks reproduce every number their prose quotes. The generalised-least-squares solve was written twice. _banded_correlation was hoisted to grid.py but the sandwich that consumes it was not, so both optimisers carried (J^T R^-1 J)^-1 and had already diverged on what a singular fit means -- one warned and returned NaN, the other silently returned the uncorrelated value. Now _gls_covariance in grid.py, with each caller keeping only its own fallback policy. What genuinely differs between the two is the Jacobian, analytic for a straight line and by finite differences for the spectral model, and that was already on the right side of the boundary. That solve also materialised a dense N x N matrix and factorised it with a general LU, for a matrix with two off-diagonals. It is banded now, via scipy.linalg.solveh_banded: 24x faster at 249 bins and 50x at 500, and the hand-rolled Toeplitz construction goes with it. Only ~2% of a fit, so it will not show on a stopwatch, but it scales cubically in bin count. profile spent 47 nested fits where 30 suffice. brentq was re-evaluating scan nodes whose misfit was already in hand, and running to its default xtol of 2e-12 -- twelve significant figures on a depth whose uncertainty is 3 km, each digit costing a full re-fit. Caching the constrained fits and asking for a millimetre instead: 28% faster, interval unchanged at 3e-5 km. sensitivity restarted every simulation from the caller's guess although all of them land in the same basin; starting from the unresampled fit is 29% faster and agrees to 1e-3. seed was being handled at the wrong level, and provably so: optimise gained a seed parameter it accepts and ignores, purely so parallelise_routine could pass one uniformly -- and only on the Bouligand sibling, so the identical call raised from inside np.hanning on the Tanaka one, through a path its own docstring advertises. The decision now lives once, in parallelise_routine, which consults a `stochastic` marker. Deterministic routines need no parameter to absorb, and passing a seed to one says so instead of failing. Dead options removed. `temperature` had no caller and its default made `T = 1.0 ** anything` -- fifteen lines of docstring defending a feature the same docstring tells you not to use. `objective_routine` gained a prior parameter in a method that has had no callers since residuals inlined the prior rows. _DEFAULT_DOF was byte-identical to _TAPER_DOF[None], which the lookup already reaches. Also `nbands`/`limit`, `_jacobian`'s step, `_covariance`'s prior, `_proposal_cholesky`'s ndim, `_profiled_misfit`'s bounds, `_profile_bracket`'s k, an unreachable acceptance branch, an unused `__all__`, and an unused numpy.testing import. _profile_root handled the below/above symmetry three times and took ten arguments to rebuild a closure the caller already had; it takes three now and the caller supplies the reflection. _profiled_misfit's two expand closures were one assignment apart. The four copies of the spectrum call are a private _spectrum wrapper, mirroring the Tanaka sibling, so power=2.0 is stated once. Test suite 29.4s -> 25.4s despite gaining tests. test_routines mapped 961 centroids overlapping by 97% to smoke-test the parallel path; 121 exercise it identically. Synthetic fields are built through one cached factory in conftest rather than two near-identical helpers, which stops seven regenerations of fields already in memory. The duplicated test_dof_factor_deflates_counts is now one test, in test_grid.py where the helper lives -- which also removes the last reason for the re-export shim in optimise_tanaka. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PzxnSJBbkFRdHEF32gHvbo
dz was unbounded above, so on a window too small to resolve the long wavelengths the likelihood flattens and the fit wanders into depths that cannot describe a planet -- a sampler on a 500 km window reached 164 km against a Curie depth that sits in the mid crust. Worse, it could in principle reach the point where bouligand2009 stops returning a number at all, which is why `_OVERFLOW_RESIDUAL`, a `_COSH_OVERFLOW` clamp on the profile scan and a blanket warning filter all existed. The bound goes where the arithmetic fails, not where the physics does. Placing it at a physically plausible ceiling would be the obvious move and is the wrong one: dz has a long upper tail, so a bound anywhere near the real range clips that tail and piles probability against the wall, reporting a distribution shaped by the constraint rather than by the data. bouligand2009 overflows around |k|dz = 710, and the radial spectrum reaches the Nyquist wavenumber pi/dx whatever the window, so the ceiling follows from the grid spacing alone: 446 km at 2 km, 111 km at 500 m. Hundreds of km, far past anything reported anywhere, so it never binds on data that constrain the base -- and when it does bind, that is a window too small to see a source rather than a deep one, which `_warn_on_bounds` now says. Measured on a 500 km window, where dz is genuinely unconstrained: the chain's median is 45.9 km with a 95% upper of 179 and a maximum of 215, and *nothing* sits within 1% of the bound. The tail is reported, not truncated. Point estimates are unchanged everywhere -- the deepest fit in the repository is 45.94 km on a deliberately-too-small window. The profile scan clamps to the same ceiling, so `_COSH_OVERFLOW` is a bound rather than a special case bolted onto the scan. `_OVERFLOW_RESIDUAL` stays, demoted to a backstop for anyone calling `residuals` directly, since inside the bound the model always evaluates. Also fixes a leak in test_warns_when_a_parameter_hits_a_bound, which restored hand-written bounds onto a module-scoped grid instead of the ones it found. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PzxnSJBbkFRdHEF32gHvbo
Orients a new session on the things that are not evident from reading the code: that wavenumbers are rad/km throughout, that `power` selects between the log power and log amplitude spectrum, that Tanaka's bands have no defaults on purpose, and that `window_spectrum` rather than `radial_spectrum` is what a fit should consume. Records the two effective-degrees-of-freedom corrections and where they live, the synthetic generator's recovery caveats (dz scatters, C is a nuisance parameter), the spawn and seeding contract in parallel.py, and that notebooks are committed without outputs and quote their own numbers in prose. Also lists three known defects, including that CurieOptimiseBouligand(..., max_processors=N) is silently ignored: the assignment sits after a return in _max_thickness, where it is unreachable. Verified at runtime -- it reports cpu_count() regardless, while the Tanaka sibling honours the argument. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
The assignment sat after the return in _max_thickness, where it could never run, so the constructor keyword was silently discarded and every routine used cpu_count() regardless. Moved into __init__, matching CurieOptimiseTanaka. Nothing caught this because parallelise_routine still worked -- it just would not run serially when asked, which matters when the default spawn start method cannot reconstruct a taper or process_subgrid defined in a notebook cell. Verified that both paths now agree: 9 centroids give identical results at max_processors=1 and 4. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
recursive-include does not consult .gitignore, so the sdist shipped whatever happened to be in the working tree. Measured on a clean build: 35 notebooks, 20 of them untracked -- Jupyter checkpoints, a scratch analysis, a work-in-progress variant, and three collaborator notebooks under Examples/for-ben. The sdist is now 11 notebooks and 1.3 MB rather than 35 and 2.9 MB, and every file in it is tracked. The eleven are listed one per line so that shipping a notebook is a deliberate act. global-exclude catches checkpoints reached through the remaining globs. Verified by building before and after and diffing the members against git ls-files. Note that a stale pycurious.egg-info/SOURCES.txt is regenerated from the previous manifest and quietly reinstates excluded files, which made the fix look incomplete until it was removed; the recipe is recorded in CLAUDE.md. Also records a defect this surfaced but did not cause: install_documentation cannot work from an installed package, because Examples/ lives outside the package directory and so is absent from the wheel. Confirmed against a build of the previous manifest too. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013umvM9QpgbCQNfwFAJeXd2
Masking the whole wavenumber plane once per annulus is O(nbins * N**2), which is cubic in the window, and it dominated large runs -- over 20x slower than necessary at nr = 2001. Digitize every cell once and reduce with bincount instead. Counts come out bit-identical and k, S and sigma to ~1e-14, across odd and even n, square and non-square, tapered and untapered, both powers. Fitted parameters follow to ~1e-5 on well-conditioned fits; they move further on the flat-dz realisations, but the misfit at the two solutions agrees to 1e-5, so both land in the same valley and there is no systematic bias over 20 seeds. Two departures from the version proposed in gh-41: - Two-pass variance rather than E[x**2] - E[x]**2. The one-pass form needs a np.maximum(var, 0) guard because it cancels: ln|FFT| is O(10) with O(1) scatter, so it loses three digits of sigma, and on a near-constant bin it fails outright -- 160x high at spread 1e-8, or exactly zero. sigma feeds the GLS weighting, so it is worth the second pass, which measures marginally faster anyway. - Broadcast index vectors rather than np.mgrid, which otherwise grew peak memory 41% (69 -> 98 MB at nr = 1201). Now at parity with the loop it replaces. Logging only the binned cells also closes a path where a zero outside the binned range would raise a divide-by-zero the per-bin masking never saw. The per-bin loop survives in tests/test_grid.py as the reference the vectorised form is pinned against, since it states the bin-edge semantics unambiguously -- annuli half-open above except the last, which matters because every cell on the kx or ky axis sits exactly on an edge. Refs gh-41 Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01A8ygDJNwDu5Z3bw5aWBWqc
remove_trend_linear is the default process_subgrid, so it runs once per window when computing a spectrum. It fitted the trend with lstsq over an (N**2, 3) design matrix; on a 2001x2001 window that SVD is ~half the cost of the whole window_spectrum call. Over a regular grid the centred row and column indices are mutually orthogonal and orthogonal to the constant, so the normal equations decouple: the plane's three coefficients are one mean and two 1-D inner products, no matrix and no SVD. It is the same least-squares plane (matches lstsq to ~2e-12 on square grids) and roughly an order of magnitude cheaper (~350 ms -> ~50 ms at 2001^2, single-threaded). It also fixes a latent bug: the design matrix was built from np.mgrid[0:nc, 0:nr] (shape (nc, nr)) but fitted against data of shape (nr, nc), so on a non-square grid the raveled indices no longer lined up with the data and the trend was not removed. The separable form is correct for any aspect ratio; the added test covers a pure plane on square and non-square grids and equivalence to a correctly aligned lstsq fit. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01RC1D97igBYRHCgutbPkkW4
The separable plane fit divides the row/column inner product by (i*i).sum(), which is zero along a length-1 axis, so remove_trend_linear returned all-NaN on a 1xN or Nx1 grid where the old lstsq stayed finite. A slope along a singleton axis is not identifiable anyway, so take a zero slope there instead of dividing by zero. The 2D path is unchanged (both sums are non-zero). Adds a degenerate-shape case to test_remove_trend_linear covering both orientations. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PfrEau3RKgzjLQarp7QFzu
The magnetic anomaly is real, so |FFT| is symmetric under reflection through the origin and the full fft2 computes -- and then bins -- a redundant half. Switch _FFT_spectrum to np.fft.rfft2, which returns only the non-redundant nc//2 + 1 columns, halving both the transform and the number of cells to digitize/log/bincount. _FFT_spectrum drops from ~53 ms to ~26 ms at n=1024 (about 2x), and the win flows through every window spectrum, so through optimise_routine, sensitivity and window_spectrum. The subtlety is `counts`. window_spectrum deflates sigma by it, and the _TAPER_DOF calibration folds the Hermitian factor of two straight into dof_inf (2.0 untapered), so counts must stay the full-spectrum count. Each retained rfft column is therefore weighted by how many full-spectrum columns it stands in for: interior columns twice, the self-mirrored DC and even-nc Nyquist columns once. This reproduces the full-spectrum result exactly -- it is a restructuring, not an approximation. A naive unweighted swap would halve counts and inflate every reported uncertainty by ~sqrt(2). Verified: counts bit-identical to the old fft2 path, and k/S/sigma agree to ~1e-13, across window sizes 64..2001, both tapers, on synthetics and on a few EMAG2 vertices at 200/305/400/600 km. Every determinate EMAG2 fit matches to metres; the only larger differences are at degenerate fits (zt pinned, sigma_dz of thousands of km) where the gap is round-off amplified by a flat objective and falls within the spread of re-fitting the same spectrum with 1e-13 jitter. Tests: rectangular and odd-size windows against a full-fft2 per-bin reference (the column axis rfft2 halves), and an explicit guard that the weighting reconstructs the full-spectrum counts and that the unweighted half plane does not. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Vjk2Jon5r4hrSo9GQzNN9p
Replace the stale Travis config with GitHub Actions, migrate the documentation from pdoc to a Sphinx site published to GitHub Pages, publish to PyPI over OIDC trusted publishing, and drop Docker/Binder. CI (.github/workflows/): - tests.yml runs pytest on Python 3.9-3.13 on every push and PR; the suite needs only the `test` extra, so `pip install -e ".[test]" && pytest` runs it all - docs.yml builds the Sphinx docs and deploys to GitHub Pages from master - publish.yml builds the sdist and wheel and uploads to PyPI on a GitHub release, authenticating over OIDC with no stored secret Docs (docs/): a Furo-themed Sphinx site with getting-started, software dependencies, contributing, theory, tutorials, and an autodoc API reference. The Bouligand and Tanaka notebooks are the tutorials, executed at build time except the data-heavy Ex5; copy_notebooks.py stages them into the source tree. The method narrative and math move from the module docstrings into hand-written theory pages, and default_role = "literal" keeps the remaining pdoc-flavoured docstrings rendering cleanly. Also: add __version__ and a docs extra, remove Docker/, move the internal bouligand-findings note to notes/ so it stays out of the Sphinx build, and refresh README, CONTRIBUTING and CLAUDE. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01XS8gDB88WqK8sdQc81L28c
tanaka-probabilistic is a full v2 rewrite that diverged from an ancient common ancestor (86461be), while master carried the v1 line forward independently. Every genuinely-recent master change (pyproj v6 syntax, ungrid, export/import_netcdf4, the JOSS paper, COPYING/COPYING.LESSER) is already present in v2, and the rest is v1 infrastructure v2 deliberately dropped (Docker, Travis, the Cython radon extension, legacy modules and notebooks). This -s ours merge records master as a parent so PR #44 merges cleanly, while keeping the v2 tree verbatim. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01XQjjfmSyzgF8RkyqzNJDKY
Full scope of this PRThe description above covers the packaging / CI / docs overhaul, but that's the final commit of a much larger change. This PR is really the v2 scientific rewrite of PyCurious. Here's the whole thing at a high level. The headline: everything now carries an uncertaintyThe defining change in v2 is that both methods return an estimate and an uncertainty. Anything that returned a bare number is gone. This drove a rewrite of the spectrum layer, a new covariance stack, and calibration against synthetics with a known answer. Restructured API — two optimisers over a shared grid
New module layout: Uncertainty machinery (the core new science)
Spectrum-layer correctness fixesSeveral of these fixed real, depth-biasing bugs:
Tanaka method, done properly
Synthetics + recovery tests
Parallelism robustness
Packaging / CI / docs (the part the description already covers)
Net effect on the library/tests/docs tree: ~5.4k insertions vs ~98k deletions — the large deletion is almost entirely the in-package data fixture and duplicated v1 notebooks. 🤖 Summary generated with Claude Code |
Overhauls the project for the v2 release: replaces the stale Travis config with GitHub Actions, migrates the documentation from pdoc to a Sphinx site published to GitHub Pages, publishes to PyPI over OIDC trusted publishing, and removes the Docker/Binder setup.
CI — GitHub Actions (
.github/workflows/)tests.yml— runspyteston Python 3.9–3.13 (plus macOS on 3.13) on every push and PR. The suite imports only pytest/numpy/scipy/pycurious, sopip install -e ".[test]" && pytestruns everything.docs.yml— builds the Sphinx docs on every push/PR and deploys to GitHub Pages frommaster.publish.yml— builds the sdist + wheel and uploads to PyPI on a GitHub release, authenticating over OIDC trusted publishing (no stored secret)..travis.ymldeleted.Documentation — Sphinx + Furo (
docs/)Getting started, software dependencies, contributing (MyST-includes
CONTRIBUTING.md), theory (Bouligand + Tanaka), tutorials, and an autodoc API reference. The Bouligand and Tanaka notebooks are the tutorials, executed at build time except the data-heavyEx5(which each pull ~600 MB);copy_notebooks.pystages them into the source tree. The method narrative and math move from the module docstrings into hand-written theory pages, anddefault_role = "literal"keeps the remaining pdoc-flavoured docstrings rendering cleanly.Remove Docker & Binder
Docker/deleted; README badges/sections stripped; logo URL fixed; new Tests + Docs badges added.Packaging & housekeeping
Added
__version__(viaimportlib.metadata) and adocsextra; moved the internalbouligand-findings.mdtonotes/so it stays out of the Sphinx build; refreshedREADME,CONTRIBUTING,CLAUDE, and.gitignore.Verification (local)
pytestfull suite: 100 passed;-m "not slow": 97 passed.sphinx-build: exit 0, zero warnings — notebooks executed with plots, Ex5 rendered un-executed, theory math rendered, API populated.python -m build+twine check: both PASSED; sdist ships exactly the 11 notebooks, nodocs//notes//Docker/leakage.One-time settings (done by the maintainer)
pycuriousproject: ownerbrmather, repopycurious, workflowpublish.yml, environmentpypi.🤖 Generated with Claude Code
https://claude.ai/code/session_01XS8gDB88WqK8sdQc81L28c