fix: MLLS composition counts a negative amount as not detected and keeps the other elements

Sources/SpectroEELS/EELSComposition.swift

@@ -12,6 +12,13 @@ public struct EELSCompositionPixel: Sendable, Equatable {
     public let atomicPercent: [Double]
     public let weightPercent: [Double]
     public let unavailable: EELSUnavailable?
+    /// Per symbol: true where the element's amount is at or below zero, so it is not detected and counts as 0 at% and 0 wt% (the others are
+    /// normalised among themselves). Empty where the pixel has no answer.
+    public let notDetected: [Bool]
+
+    init(atomicPercent: [Double], weightPercent: [Double], unavailable: EELSUnavailable?, notDetected: [Bool] = []) {
+        self.atomicPercent = atomicPercent; self.weightPercent = weightPercent; self.unavailable = unavailable; self.notDetected = notDetected
+    }
 }
 
 public enum EELSQuantification {
@@ -46,7 +53,8 @@ public enum EELSQuantification {
                 return EELSCompositionPixel(atomicPercent: [], weightPercent: [], unavailable: .noComposition)
             }
             return EELSCompositionPixel(atomicPercent: amounts.map { 100 * $0 / atoms },
-                                        weightPercent: zip(amounts, masses).map { 100 * $0 * $1 / weights }, unavailable: nil)
+                                        weightPercent: zip(amounts, masses).map { 100 * $0 * $1 / weights }, unavailable: nil,
+                                        notDetected: amounts.map { $0 <= 0 })
         }
     }
 }

Sources/SpectroEELS/EELSMLLS.swift

@@ -247,8 +247,13 @@ public enum EELSMLLS {
 
 extension EELSQuantification {
     /// at% and wt% of `symbols` at every pixel of an MLLS `result`, as `composition(_: EELSModelResult, symbols:)` does for the model fit: an
-    /// element's amount is its master edge's intensity (the first independent edge in fit order when an element has K and L). Negative
-    /// intensities (the fit is unconstrained) make the pixel unavailable (`noComposition`). No uncertainty on the percentages (open item).
+    /// element's amount is its master edge's intensity (the first independent edge in fit order when an element has K and L). The fit is
+    /// unconstrained, so an absent element's amount scatters around zero: an amount at or below zero is not detected (`notDetected`), counts
+    /// as 0 at% and 0 wt%, and the others are normalised among themselves (kit decision, consistent with the non-negative model fit, which
+    /// gives an absent element 0; exSpy has no EELS composition, and its EDS Cliff-Lorimer passes a negative line through as a negative
+    /// fraction). `notDetected == false` only means the amount is above zero, not above a detection limit. Clamping per pixel biases an
+    /// average of pixel compositions upward for an absent element: quantify a summed region instead. The pixel is unavailable
+    /// (`noComposition`) only when no amount is above zero. No uncertainty on the percentages (open item).
     public static func composition(_ result: EELSMLLSResult, symbols: [String]) throws -> [EELSCompositionPixel] {
         var columns: [Int] = []
         if let dup = symbols.enumerated().first(where: { i, s in symbols[..<i].contains(s) })?.element { throw EELSModelError.symbolNotFitted(dup + " (listed twice)") }
@@ -261,14 +266,18 @@ extension EELSQuantification {
             guard pixel.unavailable == nil else {
                 return EELSCompositionPixel(atomicPercent: [], weightPercent: [], unavailable: pixel.unavailable)
             }
-            let amounts = columns.map { pixel.intensities[$0].value }
+            let fitted = columns.map { pixel.intensities[$0].value }
+            // An amount at or below zero is not detected and counts as 0, as the non-negative model fit gives an absent element.
+            let notDetected = fitted.map { !($0 > 0) }
+            let amounts = fitted.map { $0 > 0 ? $0 : 0 }
             let atoms = amounts.reduce(0, +)
             let weights = zip(amounts, masses).map { $0 * $1 }.reduce(0, +)
-            guard atoms > 0, weights > 0, amounts.allSatisfy({ $0 >= 0 }) else {
+            guard atoms > 0, weights > 0 else {
                 return EELSCompositionPixel(atomicPercent: [], weightPercent: [], unavailable: .noComposition)
             }
             return EELSCompositionPixel(atomicPercent: amounts.map { 100 * $0 / atoms },
-                                        weightPercent: zip(amounts, masses).map { 100 * $0 * $1 / weights }, unavailable: nil)
+                                        weightPercent: zip(amounts, masses).map { 100 * $0 * $1 / weights }, unavailable: nil,
+                                        notDetected: notDetected)
         }
     }
 }

Tests/SpectroEELSTests/EELSMLLSTests.swift

@@ -207,7 +207,7 @@ struct EELSMLLSTests {
         #expect(bad.pixels[1].unavailable == nil)
     }
 
-    /// at% and wt% from the MLLS intensities equal the ratios of the intensities, and a negative intensity makes the pixel unavailable.
+    /// at% and wt% from the MLLS intensities equal the ratios of the intensities, with nothing flagged as not detected.
     @Test func testCompositionFromMLLS() throws {
         let t = try planted("fecu.full.given")
         let result = try EELSMLLS.fit(coreLoss: t.core, edges: t.edges, microscope: t.microscope, options: EELSMLLSOptions(exponent: .perPixel(t.truthR)))
@@ -221,5 +221,58 @@ struct EELSMLLSTests {
         }
         #expect(throws: EELSModelError.symbolNotFitted("Si")) { try EELSQuantification.composition(result, symbols: ["Si"]) }
         #expect(throws: EELSModelError.self) { try EELSQuantification.composition(result, symbols: ["Fe", "Fe"]) }
+        #expect(comp.allSatisfy { $0.notDetected == [false, false] })
+    }
+
+    /// An element whose unconstrained MLLS amount is at or below zero is not detected: 0 at% and 0 wt%, flagged in `notDetected`, and the
+    /// others are normalised among themselves (as the non-negative model fit gives such an element 0). The pixel keeps its composition.
+    /// Pixels with every amount positive are unchanged to the bit.
+    ///
+    /// The image: the planted C, N and O image with its N amplitude moved to +-1.5 % of C's by linearity. With H_CO the hat matrix of the
+    /// fit without N (power law, C, O), the CNO fit of Y - k (H_CNO Y - H_CO Y) has N amplitude (1 - k) N(Y), because H_CO Y lies in the
+    /// CNO design's span with no N part (C and O move towards the fit without N; only the sign of N matters here).
+    @Test func testCompositionWithNegativeAbsentElementKeepsTheOthers() throws {
+        let t = try planted("cno.full.given")
+        let options = EELSMLLSOptions(exponent: .perPixel(t.truthR))
+        let all = try EELSMLLS.fit(coreLoss: t.core, edges: t.edges, microscope: t.microscope, options: options)
+        let co = try EELSMLLS.fit(coreLoss: t.core, edges: [t.edges[0], t.edges[2]], microscope: t.microscope, options: options)
+        #expect(all.edges.map(\.element) == ["C", "N", "O"])
+        let n = t.core.channels
+        var counts = t.core.counts
+        for p in 0 ..< t.core.pixelCount {
+            let c = all.pixels[p].intensities[0].value, nk = all.pixels[p].intensities[1].value
+            let target = (p % 2 == 0 ? 0.015 : -0.015) * c
+            let k = 1 - target / nk
+            let full = try #require(all.modelCurve(pixel: p)), without = try #require(co.modelCurve(pixel: p))
+            for ch in 0 ..< n { counts[p * n + ch] -= k * (full[ch] - without[ch]) }
+        }
+        let image = try CoreLossSpectrumImage(nx: t.core.nx, ny: t.core.ny, channels: n, counts: counts,
+                                              energyOffsetEV: t.core.energyOffsetEV, energyScaleEV: t.core.energyScaleEV)
+        let result = try EELSMLLS.fit(coreLoss: image, edges: t.edges, microscope: t.microscope, options: options)
+        let comp = try EELSQuantification.composition(result, symbols: ["C", "N", "O"])
+        let masses = ["C", "N", "O"].map { EELSQuantification.atomicMasses[$0]! }
+        for p in 0 ..< t.core.pixelCount {
+            let a = result.pixels[p].intensities.map(\.value), c = comp[p]
+            #expect(c.unavailable == nil, Comment(rawValue: "pixel \(p)"))
+            if p % 2 == 0 {
+                #expect(a[1] > 0, Comment(rawValue: "pixel \(p): N \(a[1])"))
+                // the old rule, exactly
+                let atoms = a.reduce(0, +), weights = zip(a, masses).map { $0 * $1 }.reduce(0, +)
+                #expect(c.atomicPercent == a.map { 100 * $0 / atoms })
+                #expect(c.weightPercent == zip(a, masses).map { 100 * $0 * $1 / weights })
+                #expect(c.notDetected == [false, false, false])
+            } else {
+                #expect(a[1] < 0, Comment(rawValue: "pixel \(p): N \(a[1])"))
+                #expect(c.notDetected == [false, true, false])
+                #expect(c.atomicPercent[1] == 0 && c.weightPercent[1] == 0)
+                #expect(abs(c.atomicPercent[0] - 100 * a[0] / (a[0] + a[2])) < 1e-9)
+                #expect(abs(c.weightPercent[2] - 100 * a[2] * masses[2] / (a[0] * masses[0] + a[2] * masses[2])) < 1e-9)
+                #expect(abs(c.atomicPercent.reduce(0, +) - 100) < 1e-9 && abs(c.weightPercent.reduce(0, +) - 100) < 1e-9)
+            }
+        }
+        // every element at or below zero: no composition
+        let none = try EELSQuantification.composition(result, symbols: ["N"])
+        #expect(none[1].unavailable == .noComposition && none[1].atomicPercent.isEmpty)
+        #expect(none[0].unavailable == nil && none[0].atomicPercent == [100])
     }
 }

Tests/SpectroEELSTests/EELSModelFitTests.swift

@@ -130,7 +130,7 @@ struct EELSModelFitTests {
         let cu = ModelFixtures.nums((c.record["truth"] as! [String: Any])["intensity"].flatMap { ($0 as! [String: Any])["Cu_L3"] })
         for p in 0 ..< 6 {
             let px = comp[p]
-            #expect(px.unavailable == nil && px.atomicPercent.count == 2)
+            #expect(px.unavailable == nil && px.atomicPercent.count == 2 && px.notDetected == [false, false])
             let atFe = 100 * intensity[p] / (intensity[p] + cu[p])
             #expect(abs(px.atomicPercent[0] - atFe) < 0.05, Comment(rawValue: "pixel \(p): \(px.atomicPercent) vs \(atFe)"))
             #expect(abs(px.atomicPercent.reduce(0, +) - 100) < 1e-9 && abs(px.weightPercent.reduce(0, +) - 100) < 1e-9)
@@ -142,7 +142,7 @@ struct EELSModelFitTests {
         var mask = [Bool](repeating: true, count: 6); mask[1] = false
         let masked = try EELSModelFit.fit(coreLoss: c.core, edges: c.edges, microscope: c.microscope, mask: mask)
         let m = try EELSQuantification.composition(masked, symbols: ["Fe"])
-        #expect(m[1].unavailable == .masked && m[1].atomicPercent.isEmpty)
+        #expect(m[1].unavailable == .masked && m[1].atomicPercent.isEmpty && m[1].notDetected.isEmpty)
         #expect(abs(m[0].atomicPercent[0] - 100) < 1e-12 && abs(m[0].weightPercent[0] - 100) < 1e-12)
     }
 

docs/parity.md

@@ -218,7 +218,7 @@ differ from raw SciPy where the CPU-dependent tie rule says so.
 |---|---|---|---|---|---|---|
 | `create_model`, `EELSModel.fit(kind="smart")` (power-law background + dft edges, twins, low-loss convolution of the edges, `start_energy`) | `EELSModelFit.fit`, `EELSModelResult.modelCurve(pixel:)` | present | DEV | study | none upstream for `fit`; kit decision: intensities within rtol 1e-3 of exSpy's converged values, model curve within 1e-4 of its maximum | ModelPinsTests `testSmartFitEqualsExSpy` (5 images, noise-free and Poisson-noisy, K and L, two masters, convolved: intensities 3e-5, curve 5e-6 of the maximum measured); EELSModelFitTests `testFitRecoversPlantedIntensities` (planted intensities 5e-4, noise-free 1e-9). DEVIATION: variable projection (NNLS for A and the intensities, grid + Brent for r in [0, 10]) instead of exSpy's LM, which stops 1e-4 short; only the last cascade stage is solved (the whole range, last channel left out, as exSpy). Twins by element and family (exSpy: by element). An edge onset outside the fit range is refused (exSpy turns one below the range into background). Not offered: fine structure, free onset energy, `append`/`remove`/`enable`/`disable`/`fix`/`free`, other background components |
 | `EELSModel.two_area_background_estimation`, `fit_background`, `quantify` | `EELSModelFit.twoAreaBackground`, `.fitBackground`, `EELSModelResult.quantifyText(pixel:)` | present | TOL | study | `test_eelsmodel.py::test_two_area_powerlaw_estimation_*` allclose defaults; `test_quantify` text | ModelPinsTests `testTwoAreaBackgroundEqualsExSpy` (A, r 1e-7), `testFitBackgroundEqualsExSpy` (noise-free 1e-6; noisy: r 0.02, A 10 %: A and r trade off along a valley exSpy's LM leaves early), the quantify lines in `testSmartFitEqualsExSpy`. `resolve_fine_structure`, `remove_fine_structure_data` and the fine-structure switches: none, missing |
-| `Model1D.fit(optimizer="lstsq")` on the EELS model with the power-law exponent fixed (the multiple linear least squares; `p_std` from `_calculate_covariance`; `set_signal_range`) | `EELSMLLS.fit`, `EELSMLLSResult.modelCurve(pixel:)`, `EELSQuantification.composition(_: EELSMLLSResult, ...)` | present | TOL | study | `exspy/tests/models/test_linear_model.py::test_convolved_and_std_error` (lstsq against LM, `assert_allclose` defaults on values and `p_std`); the same defaults here, rtol 1e-7 | MLLSPinsTests `testMLLSEqualsHyperSpyLstsq` (14 runs: Fe L twin, C N O K, Fe + Cu two masters, Fe L with the low-loss convolution; whole axis and an explicit window; exponent planted or from the two-area estimate; A, intensities, `p_std`, RSS, dof, window indices, curve: worst 1.1e-4 of the allowance; the noise-free planted-exponent run, whose residual is rounding noise, compared at the noise floor: RSS < 1e-20 of the data's sum of squares, standard error < 1e-9 of the amplitude); EELSMLLSTests (planted recovery, errors against the scatter of 300 Gaussian-noise realisations within 12 %, badly scaled columns, references, window, refusals, masks, composition). DEVIATIONS: columns scaled to unit length (numpy's `lstsq(rcond=None)` drops the power law for a low loss of area 1e6); a singular design is `.fitFailed` (HyperSpy: `LinAlgError`); the unconstrained fit can return a negative intensity (as `lstsq`). The standard errors assume constant variance (HyperSpy's too): on Poisson data they overstate the background's amplitude error (7x on the pinned Fe L spectrum). Reference shapes (`references`) have no HyperSpy counterpart pinned (planted truth only). Not offered: weights, `ols` / `nnls` / `ridge`, a free exponent (that is `EELSModelFit.fit`) |
+| `Model1D.fit(optimizer="lstsq")` on the EELS model with the power-law exponent fixed (the multiple linear least squares; `p_std` from `_calculate_covariance`; `set_signal_range`) | `EELSMLLS.fit`, `EELSMLLSResult.modelCurve(pixel:)`, `EELSQuantification.composition(_: EELSMLLSResult, ...)` | present | TOL | study | `exspy/tests/models/test_linear_model.py::test_convolved_and_std_error` (lstsq against LM, `assert_allclose` defaults on values and `p_std`); the same defaults here, rtol 1e-7 | MLLSPinsTests `testMLLSEqualsHyperSpyLstsq` (14 runs: Fe L twin, C N O K, Fe + Cu two masters, Fe L with the low-loss convolution; whole axis and an explicit window; exponent planted or from the two-area estimate; A, intensities, `p_std`, RSS, dof, window indices, curve: worst 1.1e-4 of the allowance; the noise-free planted-exponent run, whose residual is rounding noise, compared at the noise floor: RSS < 1e-20 of the data's sum of squares, standard error < 1e-9 of the amplitude); EELSMLLSTests (planted recovery, errors against the scatter of 300 Gaussian-noise realisations within 12 %, badly scaled columns, references, window, refusals, masks, composition). DEVIATIONS: columns scaled to unit length (numpy's `lstsq(rcond=None)` drops the power law for a low loss of area 1e6); a singular design is `.fitFailed` (HyperSpy: `LinAlgError`); the unconstrained fit can return a negative intensity (as `lstsq`); its composition counts an amount at or below zero as not detected (0 at% and 0 wt%, `notDetected`) and normalises the others, a kit decision consistent with the non-negative model fit (exSpy has no EELS composition; its EDS Cliff-Lorimer passes a negative line through as a negative fraction once two lines exceed its 0.1 floor) (`testCompositionWithNegativeAbsentElementKeepsTheOthers`). The standard errors assume constant variance (HyperSpy's too): on Poisson data they overstate the background's amplitude error (7x on the pinned Fe L spectrum). Reference shapes (`references`) have no HyperSpy counterpart pinned (planted truth only). Not offered: weights, `ols` / `nnls` / `ridge`, a free exponent (that is `EELSModelFit.fit`) |
 | (kit) `EELSQuantification.composition`, `LowLossSpectrumImage.sum(mask:)`, `CoreLossSpectrumImage.sum(mask:)` | same | present | NP | study | not in exSpy (macTEM design 2.3): at% and wt% from the master intensities with the kit's atomic masses; the masked-in spectra added channel by channel | EELSModelFitTests `testCompositionFromFittedIntensities`, `testRegionSumAndItsFit` |
 | `EELSCLEdge.function` (energy-differential cross section, dft GOS, K shell) | `EELSEdgeQuantification.differentialCrossSection` (public), `EdgeCrossSection` (package) | present | TOL | study | `tests/components/test_eels_cl_edge.py` rtol 5e-7 (model restore); GOS values exact tables | EdgePinsTests `testCrossSectionEqualsExSpyOnTheGrid` (192 cases x 3 scales x 4 widths: cross section 1e-10 bound, 7e-16 measured), `testSplineAndQLimitsEqualExSpy`; second oracle `testCrossSectionWithinThreePercentOfPyEELSMODEL` (72 cases). The 17 task-16 edges (K of Mg, Al, Si; L1, L2, L3 of Mg, Al, Si, Fe, Cu; Al as the combined L2,3): EdgeTablePinsTests `testBundledTablesAreBitExact` (table sha256, float32 bytes) and `testNewEdgesEqualExSpy` (2 settings each; angle 0, qint 5e-16, function 6e-16 measured against the 1e-14 / 1e-12 / 1e-10 bounds). Every other K, L and M edge exSpy builds with the dft GOS (0.12.0: 383 edges from 231 tables, the pack `dft-gos-tables.bin`): EdgePackPinsTests `testBundledTablesAreBitExact` (all 231 tables, sha256) and `testOneEdgePerShellEqualsExSpy` (B K; Ti L1, L2, L3; Zr M2, M3; Ba M4, M5; 2 settings each: angle 0, qint 3e-16, function 6e-16 measured). Registration: design note eels-slice-2, archived in the suite's private folder |
 | `EELSCLEdge` with `GOS="Hartree-Slater"` | none | not planned | NP | not planned: Gatan's Hartree-Slater tables are never shipped | `tests/misc/test_gos.py::test_hartree_slater_gos` | Gatan's tables are never shipped (design note eels-slice-2, section 1, archived in the suite's private folder) |