From 6792a47f61ad54f0bc9cb178c3c031237ec13e50 Mon Sep 17 00:00:00 2001 From: Paullo9 Date: Sun, 1 Nov 2026 18:15:07 +0100 Subject: [PATCH] fix: MLLS composition counts a negative amount as not detected and keeps the other elements Found in an app EELS room: with C K, N K and O K picked and N absent, the unconstrained MLLS amount of N came out below zero, and composition(_: EELSMLLSResult, symbols:) refused the whole pixel, so C and O lost their at% and wt% too. An amount at or below zero now counts as 0 at% and 0 wt%, is flagged in the new EELSCompositionPixel.notDetected, and the others are normalised among themselves. This is a kit decision consistent with the non-negative model fit, which gives an absent element 0. exSpy has no EELS composition; its EDS Cliff-Lorimer passes a negative line through as a negative fraction. noComposition now means no amount is above zero. Pixels with every amount positive are bit-identical. The model-fit composition fills notDetected; its numbers do not change. Done: per-element not-detected rule for MLLS composition, red-first test. Next: the app shows not-detected elements in the footer and the table. Co-Authored-By: Claude Opus 5.5 --- Sources/SpectroEELS/EELSComposition.swift | 10 +++- Sources/SpectroEELS/EELSMLLS.swift | 19 +++++-- Tests/SpectroEELSTests/EELSMLLSTests.swift | 55 ++++++++++++++++++- .../SpectroEELSTests/EELSModelFitTests.swift | 4 +- docs/parity.md | 2 +- 5 files changed, 80 insertions(+), 10 deletions(-) diff --git a/Sources/SpectroEELS/EELSComposition.swift b/Sources/SpectroEELS/EELSComposition.swift index 89ddf8c..7599e26 100644 --- a/Sources/SpectroEELS/EELSComposition.swift +++ b/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 }) } } } diff --git a/Sources/SpectroEELS/EELSMLLS.swift b/Sources/SpectroEELS/EELSMLLS.swift index cd9c868..86b1a21 100644 --- a/Sources/SpectroEELS/EELSMLLS.swift +++ b/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[.. 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) } } } diff --git a/Tests/SpectroEELSTests/EELSMLLSTests.swift b/Tests/SpectroEELSTests/EELSMLLSTests.swift index 88f8027..999f785 100644 --- a/Tests/SpectroEELSTests/EELSMLLSTests.swift +++ b/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]) } } diff --git a/Tests/SpectroEELSTests/EELSModelFitTests.swift b/Tests/SpectroEELSTests/EELSModelFitTests.swift index da4334e..3310dd0 100644 --- a/Tests/SpectroEELSTests/EELSModelFitTests.swift +++ b/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) } diff --git a/docs/parity.md b/docs/parity.md index c8d5d3c..25196b9 100644 --- a/docs/parity.md +++ b/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) |