diff --git a/Sources/SpectroEELS/EELSDeconvolution.swift b/Sources/SpectroEELS/EELSDeconvolution.swift index 8190218..c5abc71 100644 --- a/Sources/SpectroEELS/EELSDeconvolution.swift +++ b/Sources/SpectroEELS/EELSDeconvolution.swift @@ -161,13 +161,22 @@ public enum EELSDeconvolution { let stops = try zeroLossStops(options.threshold, lowLoss: lowLoss, axis: originalAxis, cancelled: cancelled, progress: progress.map { p in { @Sendable f in p(0.2 * f) } }) + // A pixel whose continuation is not finite (A x^-r overflows) is NaN in exSpy, which deconvolves the others; here it is flagged, + // its continuation zeroed so the image can be formed, and it is left unanswered below. + var overflowed = [Bool](repeating: false, count: pixels) let cl: CoreLossSpectrumImage if options.extrapolateCoreLoss { - cl = try EELSSpectrumOps.powerLawExtrapolation(coreLoss, windowChannels: 20, extrapolationChannels: 100, cancelled: cancelled).image + let e = try EELSSpectrumOps.powerLawExtrapolationFlaggingNonFinite(coreLoss, windowChannels: 20, extrapolationChannels: 100, cancelled: cancelled) + cl = e.image + for p in 0 ..< pixels where e.nonFinite[p] { overflowed[p] = true } } else { cl = coreLoss } let ll: LowLossSpectrumImage + var firstLowLossOverflowed = false if options.extrapolateLowLoss { - ll = try EELSSpectrumOps.powerLawExtrapolation(lowLoss, windowChannels: 100, extrapolationChannels: 100, cancelled: cancelled).image + let e = try EELSSpectrumOps.powerLawExtrapolationFlaggingNonFinite(lowLoss, windowChannels: 100, extrapolationChannels: 100, cancelled: cancelled) + ll = e.image + for p in 0 ..< pixels where e.nonFinite[p] { overflowed[p] = true } + firstLowLossOverflowed = e.nonFinite[0] } else { ll = lowLoss } let (llT, llTaper) = try EELSSpectrumOps.hanningTaper(ll) let (clT, _) = try EELSSpectrumOps.hanningTaper(cl) @@ -179,6 +188,8 @@ public enum EELSDeconvolution { let fwhm: Double if let f = options.fwhmEV { fwhm = f } else { + // exSpy measures the first pixel; with its continuation non-finite the width is NaN and so is every pixel's result. + if firstLowLossOverflowed { throw EELSDeconvolutionError.noPeakWidth } let first = Array(llT.counts[0 ..< llChannels]) guard let w = try peakWidth(first, offsetEV: llT.energyOffsetEV, scaleEV: llT.energyScaleEV) else { throw EELSDeconvolutionError.noPeakWidth @@ -202,7 +213,7 @@ public enum EELSDeconvolution { var out = [Double](repeating: 0, count: pixels * origin) let why = EELSSlots(count: pixels, initial: nil) - let clSource = clT.counts, llSource = llT.counts + let clSource = clT.counts, llSource = llT.counts, unanswered = overflowed try out.withUnsafeMutableBufferPointer { dst in nonisolated(unsafe) let target = dst try clSource.withUnsafeBufferPointer { csrc in @@ -211,6 +222,7 @@ public enum EELSDeconvolution { try EELSParallel.forEachPixel(count: pixels, cancelled: cancelled, progress: progress.map { p in { @Sendable f in p(0.2 + 0.8 * f) } }) { p in guard let stopChannel = stops[p], stopChannel <= llChannels else { why.set(p, .noThreshold); return } + if unanswered[p] { why.set(p, .nonFiniteResult); return } let i0 = NumpySum.pairwise(linput.baseAddress! + p * llChannels, count: stopChannel) let jk = realFFT(cinput.baseAddress! + p * clChannels, count: clChannels, size: size, plan: forward) let jl = realFFT(linput.baseAddress! + p * llChannels, count: llChannels, size: size, plan: forward) diff --git a/Sources/SpectroEELS/EELSSpectrumOps.swift b/Sources/SpectroEELS/EELSSpectrumOps.swift index df32689..dd49d70 100644 --- a/Sources/SpectroEELS/EELSSpectrumOps.swift +++ b/Sources/SpectroEELS/EELSSpectrumOps.swift @@ -34,6 +34,27 @@ public nonisolated enum EELSSpectrumOps { _ image: I, windowChannels: Int = 20, extrapolationChannels: Int = 1024, fixNegativeR: Bool = false, cancelled: @Sendable () -> Bool = SpectroCancellation.task, progress: (@Sendable (Double) -> Void)? = nil ) throws -> (image: I, fits: [PowerLawParameters?]) { + let r = try extrapolate(image, windowChannels: windowChannels, extrapolationChannels: extrapolationChannels, fixNegativeR: fixNegativeR, + zeroNonFinite: false, cancelled: cancelled, progress: progress) + return (r.image, r.fits) + } + + /// The same, for the Fourier-ratio deconvolution: where a pixel's continuation is not finite (`A x^-r` overflows although A and r are + /// finite), its new channels are zeros and `nonFinite[p]` is true, so the caller can leave that pixel unanswered, as exSpy's NaN does, + /// instead of the whole image failing. The other pixels are exactly those of `powerLawExtrapolation`. + package static func powerLawExtrapolationFlaggingNonFinite( + _ image: I, windowChannels: Int, extrapolationChannels: Int, + cancelled: @Sendable () -> Bool = SpectroCancellation.task + ) throws -> (image: I, nonFinite: [Bool]) { + let r = try extrapolate(image, windowChannels: windowChannels, extrapolationChannels: extrapolationChannels, fixNegativeR: false, + zeroNonFinite: true, cancelled: cancelled, progress: nil) + return (r.image, r.nonFinite) + } + + private static func extrapolate( + _ image: I, windowChannels: Int, extrapolationChannels: Int, fixNegativeR: Bool, zeroNonFinite: Bool, + cancelled: @Sendable () -> Bool, progress: (@Sendable (Double) -> Void)? + ) throws -> (image: I, fits: [PowerLawParameters?], nonFinite: [Bool]) { let n = image.channels guard extrapolationChannels >= 1 else { throw EELSDeconvolutionError.extrapolationWindowTooWide(window: windowChannels, channels: n) } guard windowChannels >= 3, windowChannels <= n else { throw EELSDeconvolutionError.extrapolationWindowTooWide(window: windowChannels, channels: n) } @@ -45,6 +66,7 @@ public nonisolated enum EELSSpectrumOps { let pixels = image.nx * image.ny var counts = [Double](repeating: 0, count: pixels * extended) let fits = EELSSlots(count: pixels, initial: nil) + let nonFinite = EELSSlots(count: pixels, initial: false) let scale = image.energyScaleEV let source = image.counts try counts.withUnsafeMutableBufferPointer { dst in @@ -58,16 +80,22 @@ public nonisolated enum EELSSpectrumOps { guard case .success(let fit) = PowerLawFit.estimate(UnsafeBufferPointer(start: row, count: n), window: window) else { return } fits.set(p, PowerLawParameters(a: fit.a, r: fit.r)) let a = fixNegativeR && fit.r <= 0 ? 0.0 : fit.a + var finite = true for c in n.. (ll: LowLossSpectrumImage, core: CoreLossSpectrumImage) { + let n = 256, pixels = 12 + var ll = [Double](repeating: 0, count: pixels * n), cl = ll + for p in 0 ..< pixels { + for c in 0 ..< n { + let e = -5 + 0.5 * Double(c) + var s = 1000 * exp(-0.5 * (e / 0.6) * (e / 0.6)) + s += (80 + 5 * Double(p)) * exp(-0.5 * ((e - 20) / 4) * ((e - 20) / 4)) + if e > 1 { s += 300 * pow(max(abs(e), 1), -1.5) } + ll[p * n + c] = s + let ec = 250 + 0.5 * Double(c) + cl[p * n + c] = (2 + 0.1 * Double(p)) * (ec > 284 ? pow(ec / 284, -3) : 0) + 1e-3 * pow(ec / 250, -2) + } + } + if let bad { + for c in (n - 100) ..< (n - 50) { ll[bad * n + c] = 5.4e-15 } + for c in (n - 50) ..< n { ll[bad * n + c] = 50 } + } + return (try LowLossSpectrumImage(nx: 4, ny: 3, channels: n, counts: ll, energyOffsetEV: -5, energyScaleEV: 0.5), + try CoreLossSpectrumImage(nx: 4, ny: 3, channels: n, counts: cl, energyOffsetEV: 250, energyScaleEV: 0.5)) + } + + /// exSpy 0.3.2 gives NaN at a pixel whose low-loss continuation overflows (with a RuntimeWarning) and deconvolves the rest; the kit leaves + /// that pixel unanswered (zeros, `nonFiniteResult`) and answers the other eleven, which do not change by one bit against the same image + /// without the bad pixel. + @Test func testFourierRatioPixelWithNonFiniteContinuationIsUnavailable() throws { + let clean = try overflowPair(bad: nil), broken = try overflowPair(bad: 5) + // the bad window really overflows: the public extrapolation of that image cannot be formed + #expect(throws: (any Error).self) { try EELSSpectrumOps.powerLawExtrapolation(broken.ll, windowChannels: 100, extrapolationChannels: 100) } + let options = FourierRatioOptions(threshold: .fixed(2.5)) + let before = try EELSDeconvolution.fourierRatio(coreLoss: clean.core, lowLoss: clean.ll, options: options) + let after = try EELSDeconvolution.fourierRatio(coreLoss: broken.core, lowLoss: broken.ll, options: options) + let n = clean.core.channels + #expect(after.fwhmEV == before.fwhmEV && after.fftSize == before.fftSize) + for px in 0 ..< 12 { + let row = after.image.counts[px * n ..< (px + 1) * n] + if px == 5 { + #expect(after.unavailable[px] == .nonFiniteResult && row.allSatisfy { $0 == 0 }) + } else { + #expect(after.unavailable[px] == nil, Comment(rawValue: "pixel \(px)")) + #expect(Array(row) == Array(before.image.counts[px * n ..< (px + 1) * n]), Comment(rawValue: "pixel \(px) moved")) + } + } + // exSpy 0.3.2 on the same inputs (estimated FWHM, threshold 2.5 eV): pixel 1, channels 60...63, within the kit's Fourier-ratio tolerance + let exspy = [-0.017163108412498216, -0.00234262241292717, -0.03682816896247529, 0.038049721723035654] + for (k, want) in exspy.enumerated() { + let got = after.image.counts[1 * n + 60 + k] + #expect(abs(got - want) <= 1e-7 * abs(want) + 1e-7 * 5.330745866192701, Comment(rawValue: "channel \(60 + k): \(got)")) + } + } + + /// The FWHM is estimated from the first pixel (exSpy's current spectrum). When that pixel's continuation overflows, exSpy's width is NaN + /// and so is every pixel; the kit refuses with `noPeakWidth` (pass `fwhmEV`). With a FWHM given, only that pixel is unanswered. + @Test func testFourierRatioFirstPixelNonFiniteContinuation() throws { + let broken = try overflowPair(bad: 0) + #expect(throws: EELSDeconvolutionError.noPeakWidth) { + try EELSDeconvolution.fourierRatio(coreLoss: broken.core, lowLoss: broken.ll, options: FourierRatioOptions(threshold: .fixed(2.5))) + } + let given = FourierRatioOptions(fwhmEV: 1.2, threshold: .fixed(2.5)) + let out = try EELSDeconvolution.fourierRatio(coreLoss: broken.core, lowLoss: broken.ll, options: given) + let ref = try EELSDeconvolution.fourierRatio(coreLoss: broken.core, lowLoss: try overflowPair(bad: nil).ll, options: given) + let n = broken.core.channels + #expect(out.unavailable[0] == .nonFiniteResult && out.image.counts[0 ..< n].allSatisfy { $0 == 0 }) + #expect(out.unavailable[1 ..< 12].allSatisfy { $0 == nil }) + #expect(Array(out.image.counts[n...]) == Array(ref.image.counts[n...])) + } + + /// The core loss's own continuation (window 20) overflowing at pixel 0: that pixel is unanswered, the FWHM (from the low loss) is still + /// estimated, and the other pixels are bit-identical. The window (ten channels at 0.31, ten at 50 on 5 + 0.5 c eV) gives A 2.4e-296, + /// r -140.65 in HyperSpy, whose continuation overflows. + @Test func testFourierRatioCoreLossNonFiniteContinuation() throws { + let n = 256 + func core(bad: Bool) throws -> CoreLossSpectrumImage { + var cl = [Double](repeating: 0, count: 12 * n) + for p in 0 ..< 12 { for c in 0 ..< n { cl[p * n + c] = (100 + Double(p)) * pow((5 + 0.5 * Double(c)) / 5, -1) } } + if bad { + for c in (n - 20) ..< (n - 10) { cl[c] = 0.31 } + for c in (n - 10) ..< n { cl[c] = 50 } + } + return try CoreLossSpectrumImage(nx: 4, ny: 3, channels: n, counts: cl, energyOffsetEV: 5, energyScaleEV: 0.5) + } + let low = try overflowPair(bad: nil).ll, broken = try core(bad: true) + #expect(throws: (any Error).self) { try EELSSpectrumOps.powerLawExtrapolation(broken, windowChannels: 20, extrapolationChannels: 100) } + let options = FourierRatioOptions(threshold: .fixed(2.5)) + let out = try EELSDeconvolution.fourierRatio(coreLoss: broken, lowLoss: low, options: options) + let ref = try EELSDeconvolution.fourierRatio(coreLoss: try core(bad: false), lowLoss: low, options: options) + #expect(out.fwhmEV == ref.fwhmEV) + #expect(out.unavailable[0] == .nonFiniteResult && out.image.counts[0 ..< n].allSatisfy { $0 == 0 }) + #expect(out.unavailable[1 ..< 12].allSatisfy { $0 == nil }) + #expect(Array(out.image.counts[n...]) == Array(ref.image.counts[n...])) + } + /// A zero in the zero-loss spectrum's transform: the logarithm leaves the real line, the pixel says so. @Test func testFourierLogOfZeroSpectrumIsUnavailable() throws { let t = try plural() diff --git a/docs/parity.md b/docs/parity.md index c8d5d3c..2191892 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -201,7 +201,7 @@ differ from raw SciPy where the CPU-dependent tie rule says so. | exSpy / HyperSpy | SpectroKit | state | code | priority | upstream tolerance (test file) | evidence | |---|---|---|---|---|---|---| | `fourier_log_deconvolution` (add_zlp, crop) | `EELSDeconvolution.fourierLog` | present | TOL | study | none upstream (`test_running` only runs); kit decision: per pixel rtol 1e-7 + atol 1e-7 x max\|pinned\| (`assert_allclose` defaults cannot hold near zero across two FFT libraries); fftSize, tapered channels, channel count exact | DeconvolutionPinsTests `testFourierLogEqualsExSpy` (4 variants, worst 6e-7 of the allowance); EELSDeconvolutionTests `testFourierLogReturnsPlantedSingleScattering` (planted plural scattering returns Z x S / I0 within 0.1 %). Quirk kept: `irfft` returns 2 (m - 1) channels | -| `fourier_ratio_deconvolution` (fwhm, threshold, both extrapolations) | `EELSDeconvolution.fourierRatio` | present | TOL | study | none upstream; kit tolerance as Fourier-log; fftSize, FWHM (1e-9 eV), tapered channels exact | DeconvolutionPinsTests `testFourierRatioEqualsExSpy` (6 variants, worst 0.11 of the allowance); EELSDeconvolutionTests `testFourierRatioReturnsPlantedEdge` (planted edge within 5 %). Quirks kept: the FFT size comes from the ORIGINAL core-loss length while the FFT input is the extrapolated one; an odd complex size gives an `irfft` one channel short. A pixel with no threshold gets zeros and `noThreshold` | +| `fourier_ratio_deconvolution` (fwhm, threshold, both extrapolations) | `EELSDeconvolution.fourierRatio` | present | TOL | study | none upstream; kit tolerance as Fourier-log; fftSize, FWHM (1e-9 eV), tapered channels exact | DeconvolutionPinsTests `testFourierRatioEqualsExSpy` (6 variants, worst 0.11 of the allowance); EELSDeconvolutionTests `testFourierRatioReturnsPlantedEdge` (planted edge within 5 %). Quirks kept: the FFT size comes from the ORIGINAL core-loss length while the FFT input is the extrapolated one; an odd complex size gives an `irfft` one channel short. A pixel with no threshold gets zeros and `noThreshold`. A pixel whose power-law continuation is not finite (exSpy: NaN there, the rest deconvolved) gets zeros and `nonFiniteResult` (`testFourierRatioPixelWithNonFiniteContinuationIsUnavailable`, the other pixels bit-identical); when that pixel is the first and the FWHM is estimated, exSpy is NaN everywhere and the kit throws `noPeakWidth` | | `optimal_fft_size` (scipy `next_fast_len`, real and complex) | `EELSDeconvolution.optimalFFTSize` | present | EXACT | study | `tests/misc/test_math_tools.py` | DeconvolutionPinsTests `testOptimalFFTSizeEqualsScipy` (217 sizes, both rules) | | `estimate_peak_width` (factor, FWHM roots of the interpolating cubic spline) | `EELSDeconvolution.peakWidth` | present | TOL | study | `tests/signals/test_signal1d.py` (width allclose) | DeconvolutionPinsTests `testPeakWidthEqualsExSpy` (12 pixels, 1e-9 eV; noise gives nil as exSpy's NaN); EELSDeconvolutionTests `testPeakWidthOfPlantedPeak`. `window` argument not offered | | `richardson_lucy_deconvolution` (psf, iterations) | `EELSDeconvolution.richardsonLucy` | present | TOL | study | none upstream (`exspy/tests` only checks that it is refused on a non-uniform axis); `assert_allclose` defaults for float64: rtol 1e-7, atol 0, per channel | RichardsonLucyPinsTests `testRichardsonLucyEqualsExSpy` (10 cases: shared and per-pixel PSF, noise-free and Poisson, 0 / 5 / 15 / 40 iterations, flat-topped PSF (first argmax), PSF centre at channel 3, zeros giving NaN; worst 7.7e-8 of the allowance); RichardsonLucyTests (blurred peaks recovered, one iteration by hand, delta PSF and 0 iterations are the identity, refusals). Direct convolutions as numpy (no FFT, no wrap). Needs the PSF to have the spectrum's channel count (exSpy: broadcast error otherwise). A pixel exSpy returns NaN for is zeros with a reason |