From 475083b9e3def404429eb95b1a682a4e941810d7 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 6 Nov 2026 15:51:51 +0000 Subject: [PATCH] docs: README range 0.13.2, drop the 0.9.0 EELS chain design note, reword the tools READMEs docs/design/eels-chain.md was the pre-port design, implemented and archived in the suite. The four tools READMEs point at the archived design notes instead of the registration wording. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01VkKZEV3frN2QPgMm19185h --- README.md | 2 +- docs/design/eels-chain.md | 279 --------------------------------- tools/edx-planted/README.md | 4 +- tools/eels-chain/README.md | 2 +- tools/eels-edge/README.md | 11 +- tools/eels-thickness/README.md | 12 +- 6 files changed, 16 insertions(+), 294 deletions(-) delete mode 100644 docs/design/eels-chain.md diff --git a/README.md b/README.md index 7219bdf..1d8cbaf 100644 --- a/README.md +++ b/README.md @@ -41,7 +41,7 @@ let package = Package( name: "MyTool", platforms: [.macOS(.v15)], dependencies: [ - .package(url: "https://github.com/mac4DSTEM/SpectroKit", from: "0.12.0"), + .package(url: "https://github.com/mac4DSTEM/SpectroKit", from: "0.13.2"), ], targets: [ .executableTarget(name: "MyTool", dependencies: [ diff --git a/docs/design/eels-chain.md b/docs/design/eels-chain.md deleted file mode 100644 index 901bcec..0000000 --- a/docs/design/eels-chain.md +++ /dev/null @@ -1,279 +0,0 @@ -# The EELS chain: zero-loss alignment, background, deconvolution, edge fit (design note) - -How SpectroEELS ports exSpy 0.3.2 (`exspy/signals/eels.py`, `exspy/models/eelsmodel.py`, `exspy/components/eels_cl_edge.py`) and -HyperSpy 2.4.0 (`hyperspy/_signals/signal1d.py`, `hyperspy/misc/math_tools.py`, `hyperspy/_components/power_law.py`), as installed by -`tools/reference/setup.sh` (NumPy 2.5.4, SciPy 1.18.1). Line numbers are from those files. Every result stays `validation: "none"`. -Written for 0.9.0 (2026-10-21); the code is the reference where they differ, and `docs/parity.md` holds each function's state and tolerance. - -## 1. Shared types (SpectroEELS) - -Inputs stay `LowLossSpectrumImage` and `CoreLossSpectrumImage` (same layout, `counts[(y*nx+x)*channels+c]`, uniform axis -`energyOffsetEV + energyScaleEV * c`, counts per channel). New, public: - -```swift -/// One EELS spectrum image regardless of region; LowLoss/CoreLoss images conform. Lets alignment, taper and -/// extrapolation return the caller's type. -public protocol EELSSpectrumImage: Sendable { nx, ny, channels, counts: [Double], energyOffsetEV, energyScaleEV; - init(nx:ny:channels:counts:energyOffsetEV:energyScaleEV:) throws } - -/// A per-pixel answer the method could not give, with the reason (the kit's rule: nil with a reason, never NaN). -public enum EELSUnavailable: String, Error, Sendable { case masked, nonFiniteResult, noPeak, divisionByZero, noTwoRoots } - -public enum EELSAlignmentError: Error, Equatable, Sendable { - case maskShapeMismatch(expected: Int, got: Int), alsoAlignShapeMismatch(index: Int), signalRangeOutsideAxis(fromEV: Double, toEV: Double), - emptyAfterCrop, nonUniformInput } -public enum EELSDeconvolutionError: Error, Equatable, Sendable { - case shapeMismatch(expected: Int, got: Int), kernelSizeMismatch(expected: Int, got: Int), invalidFWHM(Double), - invalidIterations(Int), extrapolationWindowTooWide(window: Int, channels: Int), tooFewChannels(Int) } -public enum EELSModelError: Error, Equatable, Sendable { - case noEdgeInRange(fromEV: Double, toEV: Double), edgeNotBundled(String), duplicateEdge(String), lowLossShapeMismatch(expected: Int, got: Int), - backgroundWindowEmpty, fitDidNotConverge(pixel: Int) } -``` - -`PixelMask = [Bool]` (SpectroCore) is the navigation mask everywhere: `true` = skip the pixel (HyperSpy's `mask`). -Progress and cancellation follow `EELSThickness.logRatio`: `cancelled: @Sendable () -> Bool = SpectroCancellation.task`, -`progress: (@Sendable (Double) -> Void)? = nil`, one check per block of pixels. Every method name string carries the -reference and its version, as `EELSThickness.methodName` does. - -## 2. Zero-loss alignment and mask - -Reference: `EELSSpectrum.align_zero_loss_peak` (exSpy `eels.py:318-483`), `estimate_zero_loss_peak_centre` (`:275-316`, -`valuemax`: first argmax), `get_zero_loss_peak_mask` (`:485-511`); HyperSpy `Signal1D.shift1D` (`signal1d.py:423-551`, -`_shift1D :252-274`: `make_interp_spline(k=1)` = linear interpolation at `axis - shift`, outside the axis NaN), -`estimate_shift1D` (`:626-736`, `_estimate_shift1D :236-249`, `interpolate1D :220-233`), `align1D` (`:738-853`). - -```swift -public struct ZeroLossAlignmentOptions: Sendable, Equatable { - public var calibrate: Bool = true // offset -= nanmean(zlpc) after each pass - public var subpixel: Bool = true // align1D over signalRangeEV ?? -3...3 eV, reference = pixel 0, linear x6 interpolation - public var signalRangeEV: ClosedRange? = nil // search window for the maximum (value slicing, HyperSpy isig end-exclusive) - public var crop: Bool = true // drop channels any shift left undefined (shift1D ilow/ihigh); false keeps size, NaN -> nil pixels - public var interpolationPoints: Int = 5 // HyperSpy number_of_interpolation_points; ip = points + 1 -} -public struct ZeroLossAlignment: Sendable { - public let aligned: LowLossSpectrumImage // new offset and channel count - public let alsoAligned: [CoreLossSpectrumImage] - public let initialCentreEV: [Double?] // zlpc before anything moved (the app's diagnostic); nil where masked - public let shiftEV: [Double?] // total shift applied per pixel (coarse + subpixel) - public let keptChannels: Range // of the input axis, after crop - public let methodName: String // "align_zero_loss_peak, exSpy 0.3.2 / HyperSpy 2.4.0" -} -public enum EELSAlignment { - public static func alignZeroLossPeak(_ image: LowLossSpectrumImage, alsoAlign: [CoreLossSpectrumImage] = [], mask: PixelMask? = nil, - options: ZeroLossAlignmentOptions = .init(), cancelled:..., progress:...) throws -> ZeroLossAlignment - public static func zeroLossPeakCentre(_ image: LowLossSpectrumImage, mask: PixelMask? = nil, signalRangeEV: ClosedRange? = nil) -> [Double?] - public static func zeroLossPeakMask(_ image: LowLossSpectrumImage, widthEV: Double = 5.0, signalMask: [Bool]? = nil) -> [Bool] // per channel - public static func estimateShift(_ image: LowLossSpectrumImage, fromEV: Double?, toEV: Double?, referencePixel: Int = 0, interpolationPoints: Int = 5, - mask: PixelMask? = nil) -> [Double?] // estimate_shift1D in eV; HyperSpy's max_shift is a no-op there (clip without assignment) and is not offered - public static func shift(_ image: I, byEV shifts: [Double?], crop: Bool = true) throws -> (image: I, keptChannels: Range) -} -``` - -Port exactly: coarse pass shifts every pixel by `mean(zlpc) - zlpc` (not to zero: to the mean), calibrate subtracts the -new nanmean centre from the offset, subpixel pass only when `pixelCount > 1`, window clamped to the axis as `:453-462`, -`np.correlate(ref, data, "full")` on mean-subtracted, linearly upsampled (`ch*ip-(ip-1)` points) signals, argmax ties to -the first index, shift = (argmax - len(ref) + 1) / ip × scale; then calibrate again. Cross-correlation on the CPU via -`vDSP_convD` in float64 (direct, so the argmax is reproducible). `zlpc` of an all-equal or all-NaN pixel: first index, as -NumPy. `interpolationPoints < 0` and an empty `signalRangeEV` are refused. - -Tolerances (upstream): `test_eels.py::TestAlignZLP` mean and std of the re-estimated centre at `assert_allclose` -defaults (rtol 1e-7, atol 0; with mask atol `eps`), `calibrate=False` std atol 1e-2, crop=False keeps size exactly, -`test_align_zero_loss_peak_with_spike_signal_range` max 8 at defaults; HyperSpy `test_1D_tools.py::TestAlignTools` -`test_estimate_shift` atol 1e-3 (eV), `test_shift1D` / `test_align` data at defaults, new offset exact. - -## 3. Power-law extrapolation and Hanning taper - -Reference: `EELSSpectrum.power_law_extrapolation(window_size=20, extrapolation_size=1024, add_noise=False, fix_neg_r=False)` -(exSpy `eels.py:1236-1331`; `PowerLaw.estimate_parameters(s, E1=index2value(size-window), E2=index2value(size-1), out=True)` -on the zero-extended signal, binned axis multiplies by the scale, `A = where(r <= 0, 0, A)` for `fix_neg_r`); -HyperSpy `Signal1D.hanning_taper(side="both", channels=None, offset=0)` (`signal1d.py:1453-1531`: `channels = max(20, -round(0.02 n))`, `np.hanning(2 channels)` halves, in place, returns `channels`). - -```swift -public enum EELSSpectrumOps { - /// exSpy power_law_extrapolation; add_noise is not offered (random). The window is the existing PowerLawFit (two-area) over the last `windowChannels`. - public static func powerLawExtrapolation(_ image: I, windowChannels: Int = 20, extrapolationChannels: Int = 1024, - fixNegativeR: Bool = false) throws -> (image: I, fits: [PowerLawParameters?]) // nil where the window sum is not positive (kit rule) - public static func hanningTaper(_ image: I, side: TaperSide = .both, channels: Int? = nil, offset: Int = 0) throws -> (image: I, taperedChannels: Int) -} -public struct PowerLawParameters: Sendable, Equatable { public let a: Double, r: Double } // HyperSpy A, r -public enum TaperSide: Sendable { case left, right, both } -``` - -`PowerLawFit` (package) becomes the engine of both; its two-area formula is already pinned (`EdgePinsTests`). The -extrapolation reuses `PowerLawWindow(axis:from:to:)` with HyperSpy's own `value2index` rounding on the extended axis. -Tolerances: `test_eels.py::TestPowerLawExtrapolation` unbinned and binned, extrapolated data atol 1e-2 against the exact -power law (`atol=10e-3`); HyperSpy `tests/component/test_power_law.py` A and r at defaults. Hanning: exact (a product with -`np.hanning`, which is `0.5 - 0.5 cos(2πk/(2c-1))`; compute it the same way, no `vDSP_hann_window` whose normalisation -differs). - -## 4. Fourier-log and Fourier-ratio deconvolution - -Reference: `fourier_log_deconvolution(zlp, add_zlp=False, crop=False)` (`eels.py:872-977`), `fourier_ratio_deconvolution(low_loss, -fwhm=None, threshold=None, extrapolate_lowloss=True, extrapolate_coreloss=True)` (`:979-1097`), `optimal_fft_size` -(HyperSpy `math_tools.py:148`, `scipy.fft.next_fast_len(size, real)`), `estimate_peak_width(factor=0.5)` (`signal1d.py:1612-1725`, -cubic `UnivariateSpline(s=0)` roots of `y - 0.5 max`, exactly two roots else NaN), HyperSpy `Gaussian` (`A/(σ√2π) exp(-(x-c)²/2σ²)`). - -```swift -public struct FourierLogOptions: Sendable, Equatable { public var addZeroLoss = false; public var crop = false } -public struct FourierRatioOptions: Sendable, Equatable { - public var fwhmEV: Double? = nil // nil: estimate_peak_width of pixel 0 of the low loss (exSpy uses the current = first spectrum) - public var threshold: ThresholdRule = .estimated() // exSpy: estimate_elastic_scattering_threshold() default; a fixed eV is .fixed - public var extrapolateLowLoss = true, extrapolateCoreLoss = true // (window 100, size 100) and (20, 100) channels, as exSpy -} -public struct DeconvolvedSpectrumImage: Sendable { - public let image: I; public let unavailable: [EELSUnavailable?] // per pixel - public let fftSize: Int; public let fwhmEV: Double?; public let taperedChannels: Int; public let methodName: String -} -public enum EELSDeconvolution { - public static func fourierLog(_ spectrum: LowLossSpectrumImage, zeroLoss: LowLossSpectrumImage, options: FourierLogOptions = .init(), - cancelled:..., progress:...) throws -> DeconvolvedSpectrumImage - public static func fourierRatio(coreLoss: CoreLossSpectrumImage, lowLoss: LowLossSpectrumImage, options: FourierRatioOptions = .init(), - cancelled:..., progress:...) throws -> DeconvolvedSpectrumImage - /// scipy.fft.next_fast_len: smallest 2,3,5-smooth (real) or 2,3,5,7,11-smooth (complex) integer >= n. Public so the pins can assert it. - public static func optimalFFTSize(_ n: Int, real: Bool) -> Int -} -``` - -Faithful details: Fourier-log takes `size = next_fast_len(zlp + self - 1, real=True)`, `j1 = z · nan_to_num(log(j/z))`, -keeps the first `self` channels, `addZeroLoss` adds the ZLP over its channels, `crop` removes the last `taperedChannels`. -Fourier-ratio takes `size = next_fast_len(ll + cl - 1)` (complex rule: `optimal_fft_size(size)` without `real`), Gaussian on -`linspace(ll.offset, ll.offset + ll.scale (size-1), size)`, result `irfft(z · jk / jl) · I0` cropped to the core-loss size. -The FFT runs on `ImagoCore.FFT2DDouble(nx: size, ny: 1)` (vDSP for 2^n × {1,3,5,15}, Bluestein otherwise) so every -`next_fast_len` size is exact in length (if the `ny: 1` plan proves awkward, an any-length `FFT1DDouble` is a one-task -ImagoKit change, not a SpectroKit copy); `irfft` of a half spectrum follows NumPy's Hermitian rule (imaginary parts of bins 0 and N/2 dropped). - -Tolerances: upstream has none (`TestFourierRatioDeconvolution::test_running` only runs). Fallback: `assert_allclose` -defaults are rtol 1e-7 with atol 0, which a deconvolved channel near zero can never meet across two FFT libraries. Kit -decision, recorded in `parity.md`: rtol 1e-7 **and** atol 1e-7 × max|pinned result| per pixel; `fftSize`, `taperedChannels` -and the FWHM root positions exact to 1e-9 eV (`UnivariateSpline` roots are a cubic solve; port `scipy`'s `sproot`-free -route: cubic not-a-knot spline through the points, roots per interval by Cardano, as `InterpolatingSpline` already builds -the spline). The planted truth test checks that deconvolving a planted plural-scattering spectrum returns -the planted single-scattering distribution within its counting noise. - -## 5. Richardson-Lucy deconvolution (ordered last) - -Reference: `richardson_lucy_deconvolution(psf, iterations=15)` (`eels.py:1099-1170`): per pixel `imax = argmax(kernel)`, -`first = convolve(kernel, result)[imax : imax+n]`, `result *= convolve(kernel[::-1], signal/first)[n-1-imax : 2n-1-imax]`, -`n = psf_size`, which the slices make equal to the spectrum length (exSpy never checks; the kit refuses a mismatch). - -```swift -public enum EELSDeconvolution { // continued - public static func richardsonLucy(_ spectrum: I, kernel: LowLossSpectrumImage, iterations: Int = 15, - cancelled:..., progress:...) throws -> DeconvolvedSpectrumImage // kernel: one spectrum (pixelCount 1) or one per pixel -} -``` - -Full linear convolution per iteration via `vDSP_convD` (float64, direct; 2048 channels × 15 iterations × 4096 pixels is -about 1 s single-threaded, parallel over pixels as `quantifyMaps` does). A division by a zero `first` gives inf/NaN in -exSpy; the kit marks the pixel `.divisionByZero` and leaves its channels unchanged. Tolerance: none upstream; same -rule as section 4 (rtol 1e-7, atol 1e-7 × max|result|); `iterations < 1` refused. - -## 6. Edge model: GOS tables, background fit, edge intensities - -### 6.1 Tables and edges - -`CoreLossEdge` carries any subshell: `tools/eels-edge/extract_gos.py` writes every K, L and M table exSpy builds with the -dft GOS into one float32 pack (`Resources/GOS/dft-gos-tables.bin` with `dft-gos-index.json`; the GOSH file's md5 checked): -383 edges of 95 elements since 0.12.0 (0.9.0 bundled 20). Onsets from exSpy's `eels_binding_energies.json` (exact, `tests/misc/test_gos.py::test_binding_energy_database`). -`EdgeCrossSection` already integrates any subshell's table (`_integrate_GOS`, `function(E)` at rtol 5e-7, `test_eels_cl_edge.py`); -the only new rule is exSpy's twin: the subshells of one element share one intensity, the master is the lexicographically -last name (`L3` over `L2` over `L1`, `eelsmodel.py:283-327`), the others' onsets ride on the master's with the fixed -table difference. - -```swift -public struct CoreLossEdge { public let element: String, subshell: String, onsetEV: Double; public static let bundled: [CoreLossEdge] - public static func bundled(_ element: String, _ subshell: String) -> CoreLossEdge?; public var family: String /* "K", "L" */ } -``` - -### 6.2 Background fit, model fit, quantify - -Reference: `EELSSpectrum.create_model(auto_background=True, auto_add_edges=True, GOS="dft")` (`eels.py:1628`), `EELSModel` -(`eelsmodel.py`): `_preedge_safe_window_width = 2` eV, `fit_background` (`:509-546`: range `[start, firstOnset - 2)`, edges -at or above it disabled), `two_area_background_estimation` (`:548-592`), `smart_fit` (`:444-480`: background, then edges one -by one over `[start, nextOnset - 2)`), `fit` = HyperSpy `Model1D.fit` (unbounded Levenberg-Marquardt, "ls" loss), -convolution with the low loss when given, `quantify` (`:673-697`: prints `intensity.value` per element, twins skipped). - -```swift -public struct EELSModelOptions: Sendable, Equatable { - public var startEnergyEV: Double? = nil // exSpy start_energy; nil = axis start - public var preEdgeSafeWindowEV: Double = 2 // _preedge_safe_window_width - public var background: BackgroundModel = .powerLaw // .powerLaw | .none (auto_background False) - public var fitKind: EELSFitKind = .smart // .smart (cascade) | .standard (all at once, from the two-area start) - public var freeOnsetEnergy: Bool = false // enable_free_onset_energy: not offered in 0.9.0; the field is reserved and refused if true -} -public struct EdgeIntensityPixel: Sendable, Equatable { - public let intensity: Double? // exSpy's intensity.value (counts × barn⁻¹ scale); nil with `unavailable` - public let arealDensityPerNM2: Double? // n from intensity by the ratioMethod unit rule (1 barn = 1e-10 nm²); see the pin - public let background: PowerLawParameters?; public let unavailable: EELSUnavailable? -} -public struct EELSModelResult: Sendable { - public let edges: [CoreLossEdge] // masters only, in fit order - public let pixels: [[EdgeIntensityPixel]] // [pixel][edge] - public let fitRangeEV: [ClosedRange] // per edge, the smart-fit window used - public let methodName: String; public let validation: String // "none" -} -public enum EELSModelFit { - public static func fit(coreLoss: CoreLossSpectrumImage, edges: [CoreLossEdge], microscope: EELSMicroscope, lowLoss: LowLossSpectrumImage? = nil, - mask: PixelMask? = nil, options: EELSModelOptions = .init(), cancelled:..., progress:...) throws -> EELSModelResult - public static func fitBackground(coreLoss:..., edges:..., options:...) throws -> [PowerLawParameters?] // fit_background alone - public static func twoAreaBackground(coreLoss:..., edges:..., startEnergyEV: Double? = nil) throws -> [PowerLawParameters?] -} -``` - -Solver: the model is linear in the edge intensities and `A` for a fixed `r`, so the kit solves `r` by a bounded 1-D -minimiser (Brent, `r ∈ [0, 10]`, start from the two-area estimate) with intensities and `A` by the existing -`LeastSquaresFit` (NNLS, variable projection) at each `r`. This is the EDX precedent (`LeastSquaresFit.swift` header: -HyperSpy's LM stops short of the minimum by 1e-4 relative) and is the `DEVIATION` to write inline. Low-loss convolution: -exSpy convolves each component with the raw low-loss counts over the model's convolution axis; port that convolution -(linear, `vDSP_convD`, axis as `set_convolution_axis`) and record the factor between exSpy's `intensity` and `n` on the -planted set before fixing `arealDensityPerNM2` (the pin script writes both numbers). - -Tolerances: `test_eelsmodel.py::test_two_area_powerlaw_estimation_{BC,C,no_edge}` A and r at defaults (pinned numbers -2.1451237089380295 / 3.0118980767392736 and the two others); `test_fit_background_{B_C,C}` offset at defaults; -`test_quantify` exact intensity text; `test_eels_cl_edge.py` rtol 5e-7 for the component. For `fit` upstream states no -tolerance; kit decision: fitted intensities and `r` within rtol 1e-3 of exSpy's converged values on the planted set (the -LM shortfall is 1e-4), and within the counting-statistics sigma of the planted truth (`EdgeTruthTests` rule). - -## 7. Reference functions and tolerances, one table - -| kit function | reference (file:lines) | upstream test and tolerance | -|---|---|---| -| `alignZeroLossPeak` | exSpy `eels.py:318-483`; HyperSpy `shift1D :423`, `estimate_shift1D :626`, `align1D :738` | `test_eels.py::TestAlignZLP` defaults / eps / atol 1e-2; `test_1D_tools.py::TestAlignTools` atol 1e-3 (shift), defaults (data) | -| `zeroLossPeakCentre` | `eels.py:275-316` | `test_estimate_zero_loss_peak_centre` defaults | -| `zeroLossPeakMask` | `eels.py:485-511` | `test_get_zero_loss_peak_mask*` exact | -| `powerLawExtrapolation` | `eels.py:1236-1331`; `power_law.py` `estimate_parameters` | `TestPowerLawExtrapolation` atol 1e-2; `test_power_law.py` defaults | -| `hanningTaper` | `signal1d.py:1453-1531` | no numeric test; exact by construction | -| `fourierLog` | `eels.py:872-977`; `math_tools.py:148` | none upstream; kit: rtol 1e-7 + atol 1e-7·max, sizes exact | -| `fourierRatio` | `eels.py:979-1097`; `signal1d.py:1612` | `test_running` only; kit rule as above, FWHM roots 1e-9 eV | -| `richardsonLucy` | `eels.py:1099-1170` | none upstream; kit rule as above | -| `twoAreaBackground`, `fitBackground` | `eelsmodel.py:548-592`, `:509-546` | `test_two_area_powerlaw_estimation_*`, `test_fit_background_*` defaults | -| `EELSModelFit.fit` | `eelsmodel.py:444-480`, `:594-671`; `eels_cl_edge.py:441-497` | component rtol 5e-7; fit: none upstream, kit rtol 1e-3 + planted sigma | -| `CoreLossEdge` L tables | `gosh_gos.py` `read_gos_data`; `eels_binding_energies.json` | `test_gos.py::test_binding_energy_database` exact; tables bit for bit | - -Pin environment: exSpy 0.3.2 and HyperSpy 2.4.0 from `References/python/bin/python -I -B` (`tools/reference/setup.sh`; Python 3.14, -NumPy 2.5.4, SciPy 1.18.1), the dft GOSH 1.5.0 file by md5 `7fee8891c147a4f769668403b54c529b`, as the existing scripts. - -## 8. Performance: Accelerate, no Metal - -Every path is per pixel on a 1-D axis of 1k to 4k channels: linear interpolation, one correlation, two FFTs and a -division, an NNLS of three to six columns. At a typical map size (64×64×2048, 8 M values) that is tens of milliseconds -per step on the CPU in float64 with the pixel loop parallel over chunks (the pattern of `EDXQuantification.quantifyMaps`, -at most 8 threads). A Metal float32 path would spend its dispatch and copy on work smaller than the copy; the "where -size pays" rule does not hold for anything in this task, so no kernel, no `.metal` option, and the result records -`backend: "cpu"` only. Richardson-Lucy at 15 iterations is the heaviest (about 1 s at that size); if a spectrum -image ten times larger arrives, the FFT-based convolution on the same plan is the first step, Metal the second. The FFT comes -from ImagoKit (`FFT2DDouble`, vDSP for 2^n × {1,3,5,15}, Bluestein for the other `next_fast_len` sizes); writing a second -Bluestein in SpectroKit was the alternative and is not worth a copy. - -## 9. Decisions a reviewer should see - -1. ImagoKit becomes a dependency of SpectroKit (FFT only). Reverse by vendoring the Bluestein plan; nothing else changes. -2. Deconvolution tolerance: rtol 1e-7 plus atol 1e-7 × max|result| where exSpy tests nothing numerically (the `assert_allclose` default - with an atol, because allclose's atol 0 is unattainable across FFT libraries for near-zero channels). -3. The model fit is an exact minimiser (variable projection + NNLS), not HyperSpy's LM; pinned at rtol 1e-3, truth at - the planted sigma. Same deviation as the EDX fit, same inline `DEVIATION` note. -4. Non-finite per-pixel outcomes are nil with a reason (`EELSUnavailable`), never NaN, as every EELS result today. -5. `max_shift` and `add_noise` are not offered (a no-op and a random draw in the reference). -6. Richardson-Lucy came last (0.10.0). -7. Since 0.12.0 every edge exSpy builds with the dft GOS is bundled, so there is no list to pick. diff --git a/tools/edx-planted/README.md b/tools/edx-planted/README.md index f4c83a9..33d3180 100644 --- a/tools/edx-planted/README.md +++ b/tools/edx-planted/README.md @@ -1,6 +1,6 @@ # tools/edx-planted: generators and pins for the the EDX checks -Written before the Swift tests, from the EDX design note. Reference venv `References/python` (tools/reference/setup.sh; HyperSpy 2.4.0, NumPy 2.5.4, +Fixtures from the EDX design note (archived in the suite's private folder). Reference venv `References/python` (tools/reference/setup.sh; HyperSpy 2.4.0, NumPy 2.5.4, SciPy 1.18.1, h5py 3.16.0), run `-I -B`. Fixtures: `Tests/SpectroEDXTests/Fixtures/` (15 MB). ## Run order (all deterministic; two full runs wrote identical bytes, sha256 compared) @@ -32,4 +32,4 @@ Docstrings of each script define every parameter; this file lists readings made ## Findings **C (counts).** NumPy decode of `fei_SI_EDS-HAADF-4detectors_2frames.emd`: 916174 / 867647 / 913682 / 865236 events per stream in name order, 3562739 total, 512 markers per stream: equals the published constants (asserted, the script stops otherwise). Streams decoded for 8 public files (4detectors, three SuperX 10x50, fei_emd_si, EELS_EDS, SI_empty selection, Test SI 16x16; the ReducedData file has an empty stream group: skipped). `fei_emd_spectrum.emd` has a stream group whose grid is unresolved (skipped, reported). `example_velox_EELS_EDS.emd`: ScanSize x ScanArea 16 x 20 vs rsciio's last-image 128 x 128 (known deviation). Planted files: decode equals the draw for every frame range (self-check, passed). exSpy `sum()` of a uint16 image: dtype uint64, equals NumPy uint64 (max 131077 from two 65535), `live_time` 2.0 -> 40.0 over 20 pixels. Predicted outcome met. **E (axis).** rsciio's `_get_dispersion_offset` (ast-cut, sha256 asserted) agrees with the hand-written expectations on all six planted files: 4 streams offsets -1.0/-1.0/-0.996/-1.0 keV (calibration must disagree), decoy picks (0.01, -0.5), axis-order picks Detector-10 (0.02, -1.9) by FILE order. Possible Swift deviation (read, not run): `VeloxEMDReader` sorts `detectors` by index before the "first whose name contains" rule, which would pick Detector-2 (10, -0.5) on `axis-order`. E3, exSpy `calibrate_energy_axis("offset")` (offset change vs -delta, eV): N 1e6 (288 fits) mean +0.03, max |err| 0.71; N 1e7 (96) mean +0.02, max 0.62. By gain: -3e-4 mean +0.46, 0 mean +0.03 (max 0.29 at 1e6, 0.15 at 1e7), +3e-4 mean -0.41. **Registered prediction "within +-0.2 eV at N_Al >= 1e6" is FALSE when the gain is off (max 0.71 eV) and marginally false at gain 0, N 1e6 (0.29 eV; measured 0.1 in the probe at one cell).** `"scale"` mode returns scale 0.0199577 for every planted gain (true 0.0200060 / 0.0200000 / 0.0199940): it does not recover gain here (and alters the cropped offset). `"resolution"` raised "FWHM_MnKa ... smaller than physically possible" in 2 of 384 fits. The Swift refinement bound (1 eV, 1e-4) is untested yet. -**P (peaks).** P1: **exSpy's `get_xray_lines_near_energy` order for exactly tied |dE| lines is not reproducible** (it iterates a hash-ordered element dict; the pins list differed between two processes). The registered "same order including ties" cannot be a contract; pins store canonical (|dE|, name) order (45 tied adjacent pairs in 360 cases). Window-method bias (exSpy windows, Mg net at rho_Mg = 0, noise-free mean of the generator): U130 family s -134 / -1335 / -13353 at N 1e5 / 1e6 / 1e7; family L -129 / -1292 / -12922; M/s -117 / -1174 / -11744. Registered 7.2 (U, 200 realisations): -127 / -1300 / -12893: reproduced. Realisations agree (e.g. U130/L 1e7 -12951 +- 534, n = 20). P3: `find_peaks1D_ohaver` with defaults finds one peak per spectrum and Mg in 0/5 at every cell; with amp 0.2 % of max, medfilt 3 it finds Mg (within 20 eV) 5/5 at rho >= 0.01 for N >= 1e5 and at rho 0.002 for N >= 1e6 (2/5 at 1e5), with 8.6 to 22 peaks per spectrum. Surrogate (probe3 design, SciPy NNLS) on the generator's noise-free means at N 1e7, Mg net/truth for rho 0.002 / 0.01: matched 0.99 / 0.99; shift +1 0.45 / 0.88, +2 0 / 0.72, +3 0 / 0.55, +5 0 / 0.16, +-10 0 / 0; Al width x1.02 0 / 0.77, x1.06 0 / 0.28; Fano (rho_c .0008) 0 / 0.10, (rho_c .008) 0 / 0.39; tail 2 % 0 / 0.77. This reproduces the 7.3 table (matched 1.00 / 0.98, +1 0.43 / 0.87, +2 0 / 0.71, +3 0 / 0.54, +5 0 / 0.15, x1.02 0 / 0.77, x1.06 / 0.28, Fano 0 / 0.39) to within 0.02, except the tail row (0.77 here against 0.81 registered) and with the matched L_D of 2009 against 2.2e3 registered. No Swift proposer was run. +**P (peaks).** P1: **exSpy's `get_xray_lines_near_energy` order for exactly tied |dE| lines is not reproducible** (it iterates a hash-ordered element dict; the pins list differed between two processes). The design note's "same order including ties" cannot be a contract; pins store canonical (|dE|, name) order (45 tied adjacent pairs in 360 cases). Window-method bias (exSpy windows, Mg net at rho_Mg = 0, noise-free mean of the generator): U130 family s -134 / -1335 / -13353 at N 1e5 / 1e6 / 1e7; family L -129 / -1292 / -12922; M/s -117 / -1174 / -11744. Design note 7.2 (U, 200 realisations): -127 / -1300 / -12893: reproduced. Realisations agree (e.g. U130/L 1e7 -12951 +- 534, n = 20). P3: `find_peaks1D_ohaver` with defaults finds one peak per spectrum and Mg in 0/5 at every cell; with amp 0.2 % of max, medfilt 3 it finds Mg (within 20 eV) 5/5 at rho >= 0.01 for N >= 1e5 and at rho 0.002 for N >= 1e6 (2/5 at 1e5), with 8.6 to 22 peaks per spectrum. Surrogate (probe3 design, SciPy NNLS) on the generator's noise-free means at N 1e7, Mg net/truth for rho 0.002 / 0.01: matched 0.99 / 0.99; shift +1 0.45 / 0.88, +2 0 / 0.72, +3 0 / 0.55, +5 0 / 0.16, +-10 0 / 0; Al width x1.02 0 / 0.77, x1.06 0 / 0.28; Fano (rho_c .0008) 0 / 0.10, (rho_c .008) 0 / 0.39; tail 2 % 0 / 0.77. This reproduces the 7.3 table (matched 1.00 / 0.98, +1 0.43 / 0.87, +2 0 / 0.71, +3 0 / 0.54, +5 0 / 0.15, x1.02 0 / 0.77, x1.06 / 0.28, Fano 0 / 0.39) to within 0.02, except the tail row (0.77 here against the design note's 0.81) and with the matched L_D of 2009 against the design note's 2.2e3. No Swift proposer was run. diff --git a/tools/eels-chain/README.md b/tools/eels-chain/README.md index 20ec31f..1d8c0ba 100644 --- a/tools/eels-chain/README.md +++ b/tools/eels-chain/README.md @@ -1,6 +1,6 @@ # tools/eels-chain: exSpy / HyperSpy pins for the EELS chain, sections 2, 3 and 4 -Written before the Swift code, so the port of `docs/design/eels-chain.md` started with numbers. +Fixtures for the EELS chain (the 0.9.0 design note is archived in the suite's private folder, `process/archive/SpectroKit/`). Covers zero-loss alignment and mask (design section 2), power-law extrapolation and Hanning taper (section 3), Fourier-log, Fourier-ratio and `optimal_fft_size` (section 4). Not covered: GOS tables, the edge model fit (section 6), Richardson-Lucy. diff --git a/tools/eels-edge/README.md b/tools/eels-edge/README.md index aef9050..3497076 100644 --- a/tools/eels-edge/README.md +++ b/tools/eels-edge/README.md @@ -1,6 +1,7 @@ # tools/eels-edge: truth and pins for EELS edge quantification (power-law background, C K areal density) -Written before the Swift code, as the EELS edge quantification design note sections 3.3, 4 and 5 fix them. Four scripts: +Four scripts. The "design note" below is the EELS edge quantification design note (eels-slice-2), archived in the suite's +private folder; its section and gate numbers are kept so the findings can be read against it. | script | needs | writes | |---|---|---| @@ -59,7 +60,7 @@ two-area estimate's I1, I2, r, A computed from the components. ## Findings -**1. Registered prediction false: the single-scattering bias is not below 1 %.** (Gate 3 predicted "the bias is below 1 %".) +**1. The design note's prediction is false: the single-scattering bias is not below 1 %.** (Gate 3 predicted "the bias is below 1 %".) The bias is absolute (atoms/nm^2) and, for single scattering, independent of n_true: `(I_b,true - I_b,est) * 1e10 / (I_0 sigma_k)`. Default window, Delta 50, sigma_k = 3850.67 barn, I_0 = 1e7: bias -0.428 (A=1e9, r=3.0), -3.548 (A=5e8, r=2.5), -0.494 (A=2e10, r=3.5) atoms/nm^2. In units of n_true the 15 pixels read (row by row, n = 2, 5, 10, 20, 50): @@ -99,7 +100,7 @@ and `energy_axis` as float32 (NumPy 2 keeps float32 + Python float float32), and factor `4 pi a0^2 R^2 / E / T * subshell_factor * 1e28` is a float32 array; only `qint` itself (and everything after the `searchsorted`/interpolation, which mix float32 data with float64 `qmin`/`qmax`) is float64. Replacing the arrays by their float64 widenings gives `qint` and `function(E)` that differ from exSpy's by up to 1.9e-7 and 1.0e-7 relative. A Double port will -miss the registered 1e-12 / 1e-10 bounds by five orders of magnitude unless it reproduces exSpy's float32 steps in the same places +miss the design note's 1e-12 / 1e-10 bounds by five orders of magnitude unless it reproduces exSpy's float32 steps in the same places (or the bounds are amended). The pins hold `tabulatedEnergies` (the float32 axis, widened), `qint`, and for C K 200/10/20 the `qmin`/`qmax` passed to `get_qaxis_and_gos` per energy (`splineC.integrateq`), to help find the first divergent step. @@ -108,7 +109,7 @@ reused across E0 at alpha = 0 (effective angle = beta) keeps a stale `qint`. The The effective angle is computed with `GOS.onset_energy` (the table's onset), not the editable onset. **8. Counting statistics, first set.** The noisy images (seed 20261015) have the edge (7.7 to 192 counts over 50 eV) on 1717 to 15038 background counts: the -ss default-window Delta-50 n of the 15 noisy pixels runs from -141.7 to +138.3 atoms/nm^2 for n_true 2-50 (sigma_n ~ 36-130 atoms/nm^2). The first set stays as registered; the +ss default-window Delta-50 n of the 15 noisy pixels runs from -141.7 to +138.3 atoms/nm^2 for n_true 2-50 (sigma_n ~ 36-130 atoms/nm^2). The first set stays as the design note defines it; the realistic set below is the one the bias and noise gates can meaningfully use. ## The second, realistic set (R): n_true 200, 500, 1000, 2000, 5000 atoms/nm^2 (a ~2-50 nm carbon film) @@ -125,7 +126,7 @@ relative to n_true, default window, Delta 50, by row (n = 200, 500, 1000, 2000, -0.214, -0.086, -0.043, -0.021, -0.009 %; -1.774, -0.71, -0.355, -0.177, -0.071 %; -0.247, -0.099, -0.049, -0.025, -0.010 %. Min/max -1.774 % / -0.0086 %: **"<= 2 % for Delta 50 with the default window" holds at every pixel** (worst: A = 5e8, r = 2.5, n = 200). It holds for all three windows and every Delta 30-150 too (worst -1.951 %, narrow window, Delta 30; default window: -1.827 % at Delta 30 ... -1.547 % at Delta 150). -The registered "below 1 %" fails at that one pixel only, at Delta 50 default (by 0.8 percentage points). +The design note's "below 1 %" fails at that one pixel only, at Delta 50 default (by 0.8 percentage points). **R2. Plural scattering (gate 4), I_low(Delta) route, default window, Delta 50**, relative bias (min .. max over the 15 pixels): tau 0.3: -9.40 % .. -7.64 %; tau 0.7: -19.35 % .. -17.57 %; tau 1.3: -32.57 % .. -30.70 %. diff --git a/tools/eels-thickness/README.md b/tools/eels-thickness/README.md index 43e6716..867cbb5 100644 --- a/tools/eels-thickness/README.md +++ b/tools/eels-thickness/README.md @@ -1,11 +1,11 @@ # tools/eels-thickness: truth and pins for EELS thickness (the t/λ thickness map) -Written before any Swift code, as the EELS thickness design note sections 4 and 5 fix them. Two generators, one -fixture directory, `Tests/SpectroEELSTests/Fixtures/` (written before the `SpectroEELS` target existed; `SpectroEELSTests` reads it). +Two generators, one fixture directory, `Tests/SpectroEELSTests/Fixtures/`. The "design note" below is the EELS thickness +design note (eels-slice-1), archived in the suite's private folder; its section numbers are kept so the findings can be read against it. **Agreement with exSpy is not truth.** `pins.json` is what exSpy 7185a4d1 answers; a port that equals it is a port, not a measurement. The truth is `planted-truth.json`, made by `make_planted.py` from NumPy alone, with no exSpy, HyperSpy or SciPy -in it. The two answer different gates (design note section 7: gate 1 against the pins, gates 2 and 3 against the truth). +in it. The two answer different tests (design note section 7: gate 1 against the pins, gates 2 and 3 against the truth). Nothing here makes a thickness validated: `validation: "none"` stays until a sample of known thickness is measured. ## Run @@ -126,7 +126,7 @@ is at or above the last channel's energy the crop ends at `size − 1`, so the l stop = `value2index(θ)` (above), end exclusive, so I_0 is the sum of channels `[0, value2index(θ))`; the channel nearest θ is not in it. If `value2index` raises because θ is above the last energy the slice runs to the end (I_0 = I_t, `axes.py:462-466`); θ below the first energy raises `IndexError` (`:457-461`); θ exactly at the last energy gives stop = `size − 1`. For the -registered 2.5, 4.5, 6.0 eV the channels are 250, 290, 320 (planted) and 205, 209, 212 (line scan). +design note's 2.5, 4.5, 6.0 eV the channels are 250, 290, 320 (planted) and 205, 209, 212 (line scan). **estimated threshold, per pixel** (`_eels.py:608-621`): NaN gives NaN; else `ind = axis.value2index(threshold)` and `I_0 = data[:ind].sum()` (numpy pairwise sum of the float64 channels `[0, ind)`; `pins.py` asserts exSpy's I_0 equals exactly @@ -180,9 +180,9 @@ noiseless image (≤ 5.6e-16: the generator and exSpy agree to rounding); and `g per pixel and route, with the pixels above 4. For the estimated threshold it is **7.4 at the worst pixel (pixels 7 and 9 exceed 4)**, for every fixed threshold at most 2.01: see the findings below. -## Findings while writing these (the design note is not amended here; these are for its amendment) +## Findings while writing these (read against the design note) -1. **Gate 3 as registered cannot hold for the estimated threshold.** With the estimated threshold the noisy image's thresholds +1. **Gate 3 as the design note states it cannot hold for the estimated threshold.** With the estimated threshold the noisy image's thresholds differ from the noiseless ones (noiseless 1.2 to 1.35 eV; noisy 1.3 to 8.65 eV, because the global `tol` is set by the noisiest pixel), so |t/λ_noisy − t/λ_noiseless| is a threshold effect, not counting noise: 6.2 σ_counting at the τ = 1.3 pixel and 7.4 σ_counting at τ = 2.0 (gate: ≤ 4). Over 200 noise realisations (NumPy seeds 1 to 200, run once with exSpy, not stored) the standard