diff --git a/Sources/SpectroEDX/Fit/EDSFit.swift b/Sources/SpectroEDX/Fit/EDSFit.swift index ff9a9f8..84d7c01 100644 --- a/Sources/SpectroEDX/Fit/EDSFit.swift +++ b/Sources/SpectroEDX/Fit/EDSFit.swift @@ -27,15 +27,18 @@ public nonisolated struct FitSettings: Sendable, Equatable { package var method: FitMethod = .leastSquares /// Opt-in low-energy tail on every line, escape and sum peak; nil (default) = Gaussians, every number unchanged. public var lineTail: LineTail? = nil + /// The table tying each family's sub-lines to its alpha line. Default `.exspy` (the parity path); `.epq` ties the K families at NIST + /// EPQ's Kβ weights. + public var lineWeights: EDSLineWeights = .exspy package init(elements: [String], resolutionMnKaEV: Double, beamEnergy: Double, fitFrom: Double? = nil, fitTo: Double? = nil, background: FitBackground, escapePeaks: Bool = true, referenceShapes: ReferenceShapes = .none, - method: FitMethod = .leastSquares, lineTail: LineTail? = nil) { + method: FitMethod = .leastSquares, lineTail: LineTail? = nil, lineWeights: EDSLineWeights = .exspy) { self.elements = elements; self.resolutionMnKaEV = resolutionMnKaEV; self.beamEnergy = beamEnergy self.fitFrom = fitFrom; self.fitTo = fitTo; self.background = background self.escapePeaks = escapePeaks; self.referenceShapes = referenceShapes; self.method = method - self.lineTail = lineTail + self.lineTail = lineTail; self.lineWeights = lineWeights } /// The upper end of the default fitted range, keV. @@ -70,7 +73,7 @@ public nonisolated struct FitSettings: Sendable, Equatable { background: .continuum(ContinuumForm(beamEnergy: beamEnergy, edges: policy.continuumEdges, orders: policy.continuumOrders, nonNegative: policy.continuumNonNegative, efficiency: efficiency, edgeStep: policy.continuumEdgeStep)), - escapePeaks: policy.escapePeaks, lineTail: policy.lineTail) + escapePeaks: policy.escapePeaks, lineTail: policy.lineTail, lineWeights: policy.lineWeights) } } @@ -155,7 +158,7 @@ package nonisolated enum EDSFit { guard !channels.isEmpty else { throw EDSFitError.noChannels } let model = EDSLineModel.build(elements: settings.elements, axis: axis, beamEnergy: settings.beamEnergy, resolutionMnKaEV: settings.resolutionMnKaEV, escapePeaks: settings.escapePeaks, - lineTail: settings.lineTail) + lineTail: settings.lineTail, lineWeights: settings.lineWeights) try validate(settings.referenceShapes, axis: axis, model: model) let design = LinearDesign.build(model: model, axis: axis, channels: channels, background: settings.background, shapes: settings.referenceShapes, counts: counts) diff --git a/Sources/SpectroEDX/Fit/EDSModel.swift b/Sources/SpectroEDX/Fit/EDSModel.swift index d907132..d241294 100644 --- a/Sources/SpectroEDX/Fit/EDSModel.swift +++ b/Sources/SpectroEDX/Fit/EDSModel.swift @@ -73,14 +73,20 @@ public nonisolated struct FitLineGroup: Sendable, Equatable { package let escapes: [FitLine] } -/// Which table ties a family's sub-lines to its alpha line. -public nonisolated enum EDSLineWeights: Sendable, Equatable { - /// eXSpy's `xray_lines.json` weights: the quantification fit (parity, D17). +/// Which table ties a family's sub-lines to its alpha line (`FitSettings.lineWeights`, `FitPolicy.lineWeights`). +/// Encoded as "exspy" or "epq". +public nonisolated enum EDSLineWeights: String, Sendable, Equatable, Codable { + /// eXSpy's `xray_lines.json` weights: the parity path (D17). case exspy - /// The identification fit (docs/design/auto-id-tem.md): the K family's Kβ from NIST EPQ `LineWeights.csv` (public domain), every K - /// transition other than K-L2 and K-L3 over those two (Al 0.0315 against eXSpy's 0.0132; two measured Al-alloy map sums read 0.031 and - /// 0.034 with Kβ free). L and M families keep eXSpy's weights (EPQ's are per subshell, with no clean map onto eXSpy's lines). - case epqK + /// NIST EPQ relative line weights for the K family (`LineWeights.csv`, public domain; N. W. M. Ritchie, NIST DTSA-II / EPQ library): + /// Kβ over Kα is every K transition other than K-L2 and K-L3, over those two (Al 0.0315, Mg 0.0174, Si 0.0571, Cu 0.146 against eXSpy's + /// 0.0132, 0.01, 0.0278, 0.132). On the public Leeds pure-Al wedge (DOI 10.5518/1826) a free Al Kβ fits to 0.028-0.030 of Kα. L and M families keep + /// eXSpy's weights (EPQ's are per subshell, with no clean map onto eXSpy's lines). The identification fit always uses it + /// (docs/design/auto-id-tem.md); `epq` is the same table. + case epqK = "epq" + + /// The NIST EPQ table (`epqK`). + public static let epq: EDSLineWeights = .epqK } public nonisolated struct EDSLineModel: Sendable, Equatable { diff --git a/Sources/SpectroEDX/Fit/EnergyAxisRefinement.swift b/Sources/SpectroEDX/Fit/EnergyAxisRefinement.swift index ece7801..ff3bc20 100644 --- a/Sources/SpectroEDX/Fit/EnergyAxisRefinement.swift +++ b/Sources/SpectroEDX/Fit/EnergyAxisRefinement.swift @@ -83,7 +83,7 @@ package nonisolated enum EnergyAxisRefinement { try EDSFit.validate(settings.referenceShapes, axis: fileAxis, model: EDSLineModel.build(elements: settings.elements, axis: fileAxis, beamEnergy: settings.beamEnergy, resolutionMnKaEV: fileResolutionMnKaEV, escapePeaks: settings.escapePeaks, - lineTail: settings.lineTail)) + lineTail: settings.lineTail, lineWeights: settings.lineWeights)) // Units: u0 = offset shift / 1 eV, u1 = gain shift / 1e-4, u2 = (R - R_file) / 1 eV. func rss(_ u: [Double]) -> Double { @@ -93,7 +93,8 @@ package nonisolated enum EnergyAxisRefinement { guard r > options.minimumResolutionEV, r < options.maximumResolutionEV else { return .infinity } let axis = EnergyAxis(offset: fileAxis.offset + d, scale: fileAxis.scale * (1 + g), size: fileAxis.size) let model = EDSLineModel.build(elements: settings.elements, axis: axis, beamEnergy: settings.beamEnergy, - resolutionMnKaEV: r, escapePeaks: settings.escapePeaks, lineTail: settings.lineTail) + resolutionMnKaEV: r, escapePeaks: settings.escapePeaks, lineTail: settings.lineTail, + lineWeights: settings.lineWeights) let design = LinearDesign.build(model: model, axis: axis, channels: channels, background: settings.background, shapes: settings.referenceShapes, counts: counts) diff --git a/Sources/SpectroEDX/Fit/SumEscapeFit.swift b/Sources/SpectroEDX/Fit/SumEscapeFit.swift index 87adc58..75b80b6 100644 --- a/Sources/SpectroEDX/Fit/SumEscapeFit.swift +++ b/Sources/SpectroEDX/Fit/SumEscapeFit.swift @@ -151,7 +151,8 @@ public nonisolated enum EDSSumEscapeFit { var s = settings if model.commonEscapeFraction { s.escapePeaks = true } // the line model needs the escape lists; the free columns are removed below let lineModel = EDSLineModel.build(elements: s.elements, axis: axis, beamEnergy: s.beamEnergy, - resolutionMnKaEV: s.resolutionMnKaEV, escapePeaks: s.escapePeaks, lineTail: s.lineTail) + resolutionMnKaEV: s.resolutionMnKaEV, escapePeaks: s.escapePeaks, lineTail: s.lineTail, + lineWeights: s.lineWeights) var design = LinearDesign.build(model: lineModel, axis: axis, channels: channels, background: s.background, shapes: s.referenceShapes, counts: counts) if model.commonEscapeFraction { diff --git a/Sources/SpectroEDX/Proposer/ElementProposer.swift b/Sources/SpectroEDX/Proposer/ElementProposer.swift index aa1087e..fe93f9b 100644 --- a/Sources/SpectroEDX/Proposer/ElementProposer.swift +++ b/Sources/SpectroEDX/Proposer/ElementProposer.swift @@ -5,7 +5,7 @@ // method (Newbury and Ritchie's automatic qualitative analysis, as AutoEMX documents its defaults; no code copied): // // 1. fit the listed elements (the landed `EDSFit`: Kramers x Bernstein continuum with the Al K step, Si-escape columns, unweighted -// non-negative least squares) with the NIST EPQ Kβ weights (`EDSLineWeights.epqK`; the quantification fit keeps eXSpy's); +// non-negative least squares) with the NIST EPQ Kβ weights (`EDSLineWeights.epqK`; the quantification fit keeps eXSpy's by default, `FitPolicy.lineWeights`); // 2. top-hat filter the residual and take its strongest peak at or above 5 σ (1.5 times that below 1 keV), the variance of a channel // being y + (0.01 mu)^2, so a strong line's small shape error is not a peak; // 3. every pool element with a line of weight >= 0.1 within 2.5 σ + 25 eV of the peak is a candidate, and so is a pile-up pair of the diff --git a/Sources/SpectroEDX/Quant/PooledQuantification.swift b/Sources/SpectroEDX/Quant/PooledQuantification.swift index e99cf35..97bebbf 100644 --- a/Sources/SpectroEDX/Quant/PooledQuantification.swift +++ b/Sources/SpectroEDX/Quant/PooledQuantification.swift @@ -262,6 +262,7 @@ public nonisolated enum PooledQuantifier { s.method = method.estimator == .leastSquares ? .leastSquares : .poissonML(policy.fit.irls) s.fitTo = range(axis).toKeV s.lineTail = policy.fit.lineTail + s.lineWeights = policy.fit.lineWeights return s } @@ -411,7 +412,8 @@ public nonisolated enum PooledQuantifier { progress?(0.95) // Plot curves: the model, and the background as the fitted continuum columns. let lineModel = EDSLineModel.build(elements: active, axis: axis, beamEnergy: beam, resolutionMnKaEV: resolution, - escapePeaks: fitSettings.escapePeaks, lineTail: fitSettings.lineTail) + escapePeaks: fitSettings.escapePeaks, lineTail: fitSettings.lineTail, + lineWeights: fitSettings.lineWeights) let design = LinearDesign.build(model: lineModel, axis: axis, channels: fit.channels, background: fitSettings.background, shapes: fitSettings.referenceShapes, counts: counts) var plotModel = [Double](repeating: 0, count: axis.size), plotBackground = plotModel diff --git a/Sources/SpectroEDX/Quant/QuantificationPolicy.swift b/Sources/SpectroEDX/Quant/QuantificationPolicy.swift index ed59b63..630e150 100644 --- a/Sources/SpectroEDX/Quant/QuantificationPolicy.swift +++ b/Sources/SpectroEDX/Quant/QuantificationPolicy.swift @@ -109,6 +109,9 @@ public nonisolated struct FitPolicy: Sendable, Equatable, Codable { /// reduced chi2 (or deviance) when that is above 1. Default 3: a stated choice, not a measured false-keep rate (c is clipped at 0, /// and the rate holds only for a model that fits; a real line under the self pair is read as pile-up whatever the level). public var sumPeakSignificance: Double + /// The table tying each family's sub-lines to its alpha line (`EDSLineWeights`). Default `.exspy`, the parity path; `.epq` ties the K + /// families at NIST EPQ's Kβ weights (Al Kβ 0.0315 of Kα against eXSpy's 0.0132). A stored policy without the field reads `.exspy`. + public var lineWeights: EDSLineWeights public init(fitFromKeV: Double = 0.2, defaultUpperLimitKeV: Double = 20, continuumOrders: [Int] = [9, 5], continuumEdges: [Double] = [FitPolicy.alKEdgeKeV], continuumNonNegative: Bool = true, @@ -116,8 +119,10 @@ public nonisolated struct FitPolicy: Sendable, Equatable, Codable { irls: IRLSSettings = IRLSSettings(), escapeRatioWarning: Double = 0.05, axisRefinement: AxisRefinementSettings = AxisRefinementSettings(), lineTail: LineTail? = nil, sumPeaks: SumPeakRule = .whenSignificant, sumPeakSignificance: Double = 3, - continuumEdgeStep: ContinuumEdgeStep = .broadened, continuumDeadLayer: ContinuumDeadLayer? = nil) { + continuumEdgeStep: ContinuumEdgeStep = .broadened, continuumDeadLayer: ContinuumDeadLayer? = nil, + lineWeights: EDSLineWeights = .exspy) { self.continuumDeadLayer = continuumDeadLayer + self.lineWeights = lineWeights self.fitFromKeV = fitFromKeV; self.defaultUpperLimitKeV = defaultUpperLimitKeV self.continuumOrders = continuumOrders; self.continuumEdges = continuumEdges self.continuumNonNegative = continuumNonNegative; self.continuumEfficiency = continuumEfficiency @@ -134,7 +139,7 @@ public nonisolated struct FitPolicy: Sendable, Equatable, Codable { private enum CodingKeys: String, CodingKey { case fitFromKeV, defaultUpperLimitKeV, continuumOrders, continuumEdges, continuumNonNegative, continuumEfficiency, escapePeaks, lineTail, maximumPolynomialOrder, irls, escapeRatioWarning, axisRefinement, sumPeaks, sumPeakSignificance, continuumEdgeStep, - continuumDeadLayer + continuumDeadLayer, lineWeights } /// The synthesized encoding in a fixed key order. Every field with a value is written, `continuumEdgeStep` and `sumPeaks` @@ -157,10 +162,11 @@ public nonisolated struct FitPolicy: Sendable, Equatable, Codable { try c.encode(sumPeakSignificance, forKey: .sumPeakSignificance) try c.encode(continuumEdgeStep, forKey: .continuumEdgeStep) try c.encodeIfPresent(continuumDeadLayer, forKey: .continuumDeadLayer) + try c.encode(lineWeights, forKey: .lineWeights) } /// A policy stored before sum peaks were the default has no `sumPeaks`: it reads `.off`, the fit it was stored with (and gets the - /// non-default footer line). One stored before `continuumEdgeStep` existed reads `.sharp`, the form it was run with. Every other + /// non-default footer line). One stored before `continuumEdgeStep` existed reads `.sharp`, the form it was run with; one without `lineWeights` reads `.exspy`. Every other /// field is required, as before. public init(from decoder: any Decoder) throws { let c = try decoder.container(keyedBy: CodingKeys.self) @@ -179,7 +185,8 @@ public nonisolated struct FitPolicy: Sendable, Equatable, Codable { sumPeaks: try c.decodeIfPresent(SumPeakRule.self, forKey: .sumPeaks) ?? .off, sumPeakSignificance: try c.decodeIfPresent(Double.self, forKey: .sumPeakSignificance) ?? 3, continuumEdgeStep: try c.decodeIfPresent(ContinuumEdgeStep.self, forKey: .continuumEdgeStep) ?? .sharp, - continuumDeadLayer: try c.decodeIfPresent(ContinuumDeadLayer.self, forKey: .continuumDeadLayer)) + continuumDeadLayer: try c.decodeIfPresent(ContinuumDeadLayer.self, forKey: .continuumDeadLayer), + lineWeights: try c.decodeIfPresent(EDSLineWeights.self, forKey: .lineWeights) ?? .exspy) } /// The Al K absorption edge, keV (X-ray Data Booklet, Table 1-1): the default continuum split. diff --git a/Sources/SpectroEDX/Quant/UnlistedLineCheck.swift b/Sources/SpectroEDX/Quant/UnlistedLineCheck.swift index 1c5b413..ab8b4e0 100644 --- a/Sources/SpectroEDX/Quant/UnlistedLineCheck.swift +++ b/Sources/SpectroEDX/Quant/UnlistedLineCheck.swift @@ -255,7 +255,8 @@ public nonisolated enum UnlistedLineChecker { // axis and width; nets from the fit's group values. let fs = q.fitSettings let listed = EDSLineModel.build(elements: fs.elements, axis: q.usedAxis, beamEnergy: fs.beamEnergy, - resolutionMnKaEV: fs.resolutionMnKaEV, escapePeaks: fs.escapePeaks, lineTail: fs.lineTail).groups + resolutionMnKaEV: fs.resolutionMnKaEV, escapePeaks: fs.escapePeaks, lineTail: fs.lineTail, + lineWeights: fs.lineWeights).groups var nets: [String: Double] = [:] for (id, v) in zip(q.fit.groupIDs, q.fit.values) { nets[id] = v } var edges: [Double] = [] diff --git a/Tests/SpectroCoreTests/Fixtures/spectro-golden/demo-edx-tiny-dm4.txt b/Tests/SpectroCoreTests/Fixtures/spectro-golden/demo-edx-tiny-dm4.txt index 7a2c21f..87743fa 100644 --- a/Tests/SpectroCoreTests/Fixtures/spectro-golden/demo-edx-tiny-dm4.txt +++ b/Tests/SpectroCoreTests/Fixtures/spectro-golden/demo-edx-tiny-dm4.txt @@ -664,6 +664,7 @@ whole.quantification.fitSettings.escapePeaks = true whole.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325 whole.quantification.fitSettings.method = .leastSquares whole.quantification.fitSettings.lineTail = nil +whole.quantification.fitSettings.lineWeights = .exspy whole.quantification.sumEscapeCounts = nil whole.quantification.fitRange.fromKeV = 3fc999999999999a 0.20000000000000001 whole.quantification.fitRange.toKeV = 40239d945918854d 9.8077724306519034 @@ -906,6 +907,7 @@ region-locked-axis.quantification.fitSettings.escapePeaks = true region-locked-axis.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325 region-locked-axis.quantification.fitSettings.method = .leastSquares region-locked-axis.quantification.fitSettings.lineTail = nil +region-locked-axis.quantification.fitSettings.lineWeights = .exspy region-locked-axis.quantification.sumEscapeCounts = nil region-locked-axis.quantification.fitRange.fromKeV = 3fc999999999999a 0.20000000000000001 region-locked-axis.quantification.fitRange.toKeV = 402380c28065dbf0 9.7514839290942916 diff --git a/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-realistic-l-sim.txt b/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-realistic-l-sim.txt index 55d1077..6de2e09 100644 --- a/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-realistic-l-sim.txt +++ b/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-realistic-l-sim.txt @@ -319,6 +319,7 @@ seed0.seven.quantification.fitSettings.escapePeaks = true seed0.seven.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325 seed0.seven.quantification.fitSettings.method = .leastSquares seed0.seven.quantification.fitSettings.lineTail = nil +seed0.seven.quantification.fitSettings.lineWeights = .exspy seed0.seven.quantification.sumEscapeCounts = nil seed0.seven.quantification.fitRange.fromKeV = 3fc999999999999a 0.20000000000000001 seed0.seven.quantification.fitRange.toKeV = 4034000000000000 20 @@ -691,6 +692,7 @@ seed0.seven-60keV.quantification.fitSettings.escapePeaks = true seed0.seven-60keV.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325 seed0.seven-60keV.quantification.fitSettings.method = .leastSquares seed0.seven-60keV.quantification.fitSettings.lineTail = nil +seed0.seven-60keV.quantification.fitSettings.lineWeights = .exspy seed0.seven-60keV.quantification.sumEscapeCounts = nil seed0.seven-60keV.quantification.fitRange.fromKeV = 3fc999999999999a 0.20000000000000001 seed0.seven-60keV.quantification.fitRange.toKeV = 404e000000000000 60 @@ -990,6 +992,7 @@ seed0.four.quantification.fitSettings.escapePeaks = true seed0.four.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325 seed0.four.quantification.fitSettings.method = .leastSquares seed0.four.quantification.fitSettings.lineTail = nil +seed0.four.quantification.fitSettings.lineWeights = .exspy seed0.four.quantification.sumEscapeCounts = nil seed0.four.quantification.fitRange.fromKeV = 3fc999999999999a 0.20000000000000001 seed0.four.quantification.fitRange.toKeV = 4034000000000000 20 diff --git a/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-unlisted-sim.txt b/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-unlisted-sim.txt index bbbb8b8..cb4b7f2 100644 --- a/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-unlisted-sim.txt +++ b/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-unlisted-sim.txt @@ -293,6 +293,7 @@ seed0.four.quantification.fitSettings.escapePeaks = true seed0.four.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325 seed0.four.quantification.fitSettings.method = .leastSquares seed0.four.quantification.fitSettings.lineTail = nil +seed0.four.quantification.fitSettings.lineWeights = .exspy seed0.four.quantification.sumEscapeCounts = nil seed0.four.quantification.fitRange.fromKeV = 3fc999999999999a 0.20000000000000001 seed0.four.quantification.fitRange.toKeV = 4034000000000000 20 @@ -1507,6 +1508,7 @@ seed0.seven.quantification.fitSettings.escapePeaks = true seed0.seven.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325 seed0.seven.quantification.fitSettings.method = .leastSquares seed0.seven.quantification.fitSettings.lineTail = nil +seed0.seven.quantification.fitSettings.lineWeights = .exspy seed0.seven.quantification.sumEscapeCounts = nil seed0.seven.quantification.fitRange.fromKeV = 3fc999999999999a 0.20000000000000001 seed0.seven.quantification.fitRange.toKeV = 4034000000000000 20 diff --git a/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-weakline-sim.txt b/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-weakline-sim.txt index 5adc0ad..71071c9 100644 --- a/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-weakline-sim.txt +++ b/Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-weakline-sim.txt @@ -213,6 +213,7 @@ level0.quantification.fitSettings.escapePeaks = true level0.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325 level0.quantification.fitSettings.method = .leastSquares level0.quantification.fitSettings.lineTail = nil +level0.quantification.fitSettings.lineWeights = .exspy level0.quantification.sumEscapeCounts = nil level0.quantification.fitRange.fromKeV = 3fc999999999999a 0.20000000000000001 level0.quantification.fitRange.toKeV = 4034050539b865ad 20.019610984343377 @@ -443,6 +444,7 @@ level4.quantification.fitSettings.escapePeaks = true level4.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325 level4.quantification.fitSettings.method = .leastSquares level4.quantification.fitSettings.lineTail = nil +level4.quantification.fitSettings.lineWeights = .exspy level4.quantification.sumEscapeCounts = nil level4.quantification.fitRange.fromKeV = 3fc999999999999a 0.20000000000000001 level4.quantification.fitRange.toKeV = 4034037c2c18d3a3 20.013613468228538 diff --git a/Tests/SpectroEDXTests/LineWeightTableTests.swift b/Tests/SpectroEDXTests/LineWeightTableTests.swift new file mode 100644 index 0000000..c417c77 --- /dev/null +++ b/Tests/SpectroEDXTests/LineWeightTableTests.swift @@ -0,0 +1,157 @@ +// +// LineWeightTableTests.swift +// Role: The opt-in NIST EPQ line-weight table of the spectrum fit (`EDSLineWeights.epq`, `FitPolicy.lineWeights`) next to +// eXSpy's (the default, the parity path), on PLANTED TRUTH ON SYNTHETIC DATA: a noise-free Al film generated here with its +// K-beta lines planted at either table (EPQ by default), and the demo absorption-truth pool re-planted at EPQ's K-beta. +// +// Why: eXSpy ties Al K-beta at 0.0132 of K-alpha, NIST EPQ at 0.0315. On the public Leeds pure-Al wedge (Snelson, Collins, +// MacArthur; DOI 10.5518/1826; CC BY 4.0) the EPQ tie lowers chi2/n over 0.8-2.6 keV from 3.15 to 2.51 and the data excess at +// 1.53-1.65 keV from +893 to -228 counts. A fit tied at one table on a spectrum planted at the other misreads a weak line beside +// Al K-alpha by tens of percent, so the table is a choice the footer names. +// + +import Foundation +import Testing +import SpectroCore +@testable import SpectroEDX + +struct LineWeightTableTests { + + // MARK: the table and the policy + + @Test func epqIsTheEPQTableAndEncodesAsEpq() throws { + #expect(EDSLineWeights.epq == .epqK) + #expect(EDSLineWeights.epq.rawValue == "epq" && EDSLineWeights.exspy.rawValue == "exspy") + #expect(try JSONDecoder().decode(EDSLineWeights.self, from: Data("\"epq\"".utf8)) == .epq) + } + + @Test func theDefaultStaysExspyAndTheModelTiesEachTable() { + #expect(FitPolicy.default.lineWeights == .exspy) + #expect(FitPolicy.exspyEquivalent.lineWeights == .exspy) + let axis = EnergyAxis(offset: 0, scale: 0.005, size: 4000) + #expect(FitSettings.standard(elements: ["Al"], axis: axis, resolutionMnKaEV: 128, beamEnergy: 200).lineWeights == .exspy) + #expect(FitSettings.standard(elements: ["Al"], axis: axis, resolutionMnKaEV: 128, beamEnergy: 200, + policy: FitPolicy(lineWeights: .epq)).lineWeights == .epq) + func kb(_ el: String, _ t: EDSLineWeights) -> Double? { + EDSLineModel.build(elements: [el], axis: axis, beamEnergy: 200, resolutionMnKaEV: 128, lineWeights: t) + .groups.first { $0.id == el + "_Ka" }?.lines.first { $0.id == el + "_Kb" }?.weight + } + #expect(kb("Al", .exspy) == 0.0132) // eXSpy xray_lines.json, unchanged (parity) + #expect(abs((kb("Al", .epq) ?? 0) - 0.0474 / 1.505) < 1e-12) // EPQ K-M3 + K-M2 + ... over K-L3 + K-L2 + #expect(abs((kb("Mg", .epq) ?? 0) - 0.017386) < 1e-5) + #expect(abs((kb("Si", .epq) ?? 0) - 0.057143) < 1e-5) + } + + @Test func aNonDefaultTableIsNamedInTheFooterAndOldPoliciesReadExspy() throws { + let epq = QuantificationPolicy(fit: FitPolicy(lineWeights: .epq)) + #expect(QuantificationPolicy.default.footerLine == nil) + #expect(epq.footerLine == "policy (non-default): fit.lineWeights = epq", "\(epq.footerLine ?? "nil")") + #expect(epq.canonicalJSON.contains("\"lineWeights\":\"epq\"")) + #expect(try JSONDecoder().decode(QuantificationPolicy.self, from: Data(epq.canonicalJSON.utf8)) == epq) + // A policy stored before the field existed reads eXSpy's table, the one it was run with. + var json = try JSONSerialization.jsonObject(with: Data(QuantificationPolicy.default.canonicalJSON.utf8)) as! [String: Any] + var fit = json["fit"] as! [String: Any] + fit["lineWeights"] = nil + json["fit"] = fit + let old = try JSONDecoder().decode(QuantificationPolicy.self, from: try JSONSerialization.data(withJSONObject: json)) + #expect(old == .default) + } + + // MARK: a synthetic Al film, K-beta planted at either table + + private static let axis = EnergyAxis(offset: 0, scale: 0.005, size: 4000) + private static let energies = (0.. Double in + guard e >= 0.1 else { return 0 } + let mu = (try? MassAbsorptionTable.bundled!.mac(element: "Al", energyKeV: e)) ?? 0 + return 9.0 * (200 - e) / e * exp(-pow(0.35 / e, 2.2)) * exp(-9e-5 * mu) + } + return energies.indices.map { i in + let e = energies[i] + let sigma = (XRayLines.fwhm(resolutionMnKaEV: resolution, atEnergy: max(e, 0.3)) ?? resolution / 1000) / (8 * log(2.0)).squareRoot() + let reach = Int((5 * sigma / axis.scale).rounded(.up)) + var sum = 0.0, ws = 0.0 + for j in max(0, i - reach)...min(axis.size - 1, i + reach) { + let w = exp(-0.5 * pow((energies[j] - e) / sigma, 2)); sum += w * raw[j]; ws += w + } + return sum / ws + } + }() + + /// Al K 3e5 counts, Si K and Mg K 2e3 each (K-alpha areas), every K-beta planted at `planted` (default the physical EPQ weights). + private static func film(planted: EDSLineWeights = .epq) -> [Double] { + func column(_ el: String) -> [Double] { + let g = EDSLineModel.build(elements: [el], axis: axis, beamEnergy: 200, resolutionMnKaEV: resolution, lineWeights: planted) + .groups.first { $0.id == el + "_Ka" }! + return EDSLineModel.column(of: g.lines, energies: energies, scale: axis.scale) + } + let al = column("Al"), si = column("Si"), mg = column("Mg") + return energies.indices.map { continuum[$0] + 3e5 * al[$0] + 2e3 * si[$0] + 2e3 * mg[$0] } + } + + private static func siBias(_ counts: [Double], fit table: EDSLineWeights) throws -> (si: Double, al: Double) { + let s = FitSettings.standard(elements: ["Al", "Si", "Mg"], axis: axis, resolutionMnKaEV: resolution, beamEnergy: 200, + policy: FitPolicy(lineWeights: table)) + let f = try EDSFit.run(counts: counts, axis: axis, settings: s, computeCovariance: false) + return (f.values[f.groupIDs.firstIndex(of: "Si_Ka")!] / 2e3 - 1, f.values[f.groupIDs.firstIndex(of: "Al_Ka")!] / 3e5 - 1) + } + + /// Red before `FitPolicy.lineWeights` reached the fit: the fit tied eXSpy's 0.0132 whatever the policy, and read Si K-alpha + /// -43.6 % on the film planted at EPQ's weights. Tied at the planted table the fit reads the same -25.4 to -25.8 % on either film + /// (the broadened-edge bias of this film, `ContinuumEdgeBroadeningTests`), so what is left is not the K-beta. + @Test func theFitTiedAtThePlantedTableReadsAWeakSiLineAsWellOnEitherFilm() throws { + let onEPQ = try Self.siBias(Self.film(), fit: .epq) + let mismatched = try Self.siBias(Self.film(), fit: .exspy) + let onExspy = try Self.siBias(Self.film(planted: .exspy), fit: .exspy) + FileHandle.standardError.write(Data(String(format: "KBETA film: EPQ-planted fit epq Si %+.4f (Al %+.5f), fit exspy %+.4f (Al %+.5f); exspy-planted fit exspy %+.4f\n", + onEPQ.si, onEPQ.al, mismatched.si, mismatched.al, onExspy.si).utf8)) + #expect(abs(onEPQ.si - onExspy.si) < 0.01, "matched ties: \(onEPQ.si) vs \(onExspy.si)") + #expect(mismatched.si < -0.40, "the eXSpy tie on an EPQ film: \(mismatched.si)") + // Al K-alpha: -0.09 % tied at the planted table, -0.91 % with the eXSpy tie (the missing K-beta area moves into the continuum). + #expect(abs(onEPQ.al) < 0.003 && mismatched.al < -0.006, "Al K-alpha: \(onEPQ.al) vs \(mismatched.al)") + } + + // MARK: the demo absorption-truth pool, K-beta re-planted at EPQ's weights + + /// The pool's noise-free expected spectrum plus (w_EPQ - w_eXSpy) x the K-alpha area of Al, Mg and Si in a Gaussian at each + /// K-beta, Cu included (the generator planted eXSpy's weights; the areas are the noise-free default fit's, within 1 % of the + /// planted ones; no K-beta here is above the Si K edge but Cu's, whose escape copy is not re-planted: 1 % of a weak line). + private static func replantedAtEPQ(_ c: DemoEDXPoolTests.Case) throws -> DemoEDXPoolTests.Case { + let base = try PooledQuantifier.run(DemoEDXPoolTests.input(c, noiseFree: true), tables: DemoEDXPoolTests.tables) + let e = (0...py ` (`tools/reference/setup.sh` first). Every -output is byte-identical to the committed fixture (checked 2026-10-28). +Run each as `EDX_SIM_KBETA=exspy References/python/bin/python -I -B tools/edx-sim/