fix: Fourier-ratio leaves a pixel with a non-finite power-law continuation unanswered

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<EELSUnavailable?>(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)

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<I: EELSSpectrumImage>(
+        _ 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<I: EELSSpectrumImage>(
+        _ 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<PowerLawParameters?>(count: pixels, initial: nil)
+        let nonFinite = EELSSlots<Bool>(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..<extended {
                         let x = axis.energy(ofChannel: c)
                         // factor * A * x ** (-r), NumPy's order of operations
                         out[c] = (scale * a) * Libm.pow(x, -fit.r)
+                        if !out[c].isFinite { finite = false }
+                    }
+                    if zeroNonFinite && !finite {
+                        for c in n..<extended { out[c] = 0 }
+                        nonFinite.set(p, true)
                     }
                 }
             }
         }
         let result = try I(nx: image.nx, ny: image.ny, channels: extended, counts: counts, energyOffsetEV: image.energyOffsetEV, energyScaleEV: scale)
-        return (result, fits.values)
+        return (result, fits.values, nonFinite.values)
     }
 
     // MARK: Hanning taper

Tests/SpectroEELSTests/EELSDeconvolutionTests.swift

@@ -93,6 +93,101 @@ struct EELSDeconvolutionTests {
         }
     }
 
+    /// A 4 x 3 low loss (256 channels from -5 eV at 0.5 eV) and core loss (256 channels from 250 eV) built in closed form. `bad` gets the
+    /// window that makes the two-area power law overflow on the continuation (A 8e-296, r -143.09 in HyperSpy): its last 100 low-loss channels
+    /// are 5.4e-15 then 50.
+    private func overflowPair(bad: Int?) throws -> (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()

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 |