Neutron Spectrum Unfolding

Neutron spectrum unfolding from activation-foil rates: SAND-II, least-squares, GRAVEL, and MAXED adjustment.

Nucleide implements neutron spectrum unfolding from activation-type measurements: a caller-supplied guess spectrum is adjusted until folding it through a caller-supplied response matrix reproduces the measured detector/reaction rates. Four methods share one iterator and diagnostics shape, landed one per cycle: SAND-II (Spectrum Analysis by Neutron Detectors II; McElroy et al., AFWL-TR-67-41, 1967)^mcelroy-1967, written clean-room from the public-domain report — the adjustment form reproduced across the open unfolding literature; STAYSL-class damped least-squares (Perey, ORNL/TM-6062, 1977); GRAVEL chi-square-weighted adjustment (Matzke, PTB-N-19, 1994); and MAXED maximum entropy (Reginatto & Goldhagen, Health Phys. 77 (1999) 579)^{reginatto-1999}, from the journal equations only — the closed UMG package is never touched. Every response value, measured rate, guess entry, and energy-group bound is caller-supplied data: evaluated libraries such as IRDFF are IAEA-copyright and are never vendored — the IRDFF-II response pack below ships as a runtime download, parsed into caller-ready rows. The IRDFF-II analytical benchmark-field shapes used in the correctness checks are re-derived from their published closed forms^trkov-2020 — the library’s tabulated group spectra are never read. The S1–S3 labels are the SAND-II pinned equations, T1–T4 the STAYSL-class set, G1–G3 the GRAVEL set, M1–M6 the MAXED set; U1–U23 label the correctness checks.

Overview

Foil-activation spectrometry — the JET/ITER diagnostics workhorse — measures a handful of detector/reaction rates and asks for the neutron spectrum behind them, one value per energy group. The forward problem is a matrix multiplication (the fold); the inverse problem has no unique answer whenever there are fewer detectors than energy groups, which is the normal case. Unfolding therefore does not invent the unseen part: it iteratively adjusts a caller-supplied guess until the fold reproduces the measurements, moving only the components the detectors can see. It is inverse-problem iteration, not a transport solve.

The Python surface nucleide.unfold (a thin facade over the crate) exposes five entry points: sandii(...), staysl(...), gravel(...), and maxed(...), which run their adjustment and return the converged spectrum with diagnostics, and forward_fold(...), the shared forward operator.

Code shape

The crate mirrors the other solver crates: sandii.rs, staysl.rs, gravel.rs, and maxed.rs hold one adjustment iterator each behind the unfold driver, lib.rs holds the shared forward operator (forward_fold) and the input validation, and error.rs is a named error set (a #[non_exhaustive] enum with a Result alias).

Mathematical details

Definitions: a response matrix R with one row per detector/reaction (entry R_ij ≥ 0 per energy group j), measured rates N_i ≥ 0, and a guess spectrum φ_j > 0 (strictly positive — the update is multiplicative).

S1 — Forward fold

c_i = \sum_j R_{ij}\,\phi_j \tag{S1}

The calculated rate c_i of detector i. This same map is how the correctness checks synthesize measured rates from a known spectrum.

W_{ji} = R_{ij}\,\phi_j / c_i \tag{S2}

The share of detector i’s calculated rate contributed by group j; the shares sum to one over j for each detector with a nonzero fold.

S3 — Multiplicative adjustment

\phi_j \leftarrow \phi_j \cdot \exp\!\left( \frac{\sum_i W_{ji}\,\ln(N_i / c_i)}{\sum_i W_{ji}} \right) \tag{S3}

One adjustment per iteration: the new spectrum is the old one times a weighted geometric mean of the per-detector correction factors N_i / c_i, weighted by the S2 rate shares. Three properties follow directly:

Positivity: a strictly positive guess stays non-negative — and stays strictly positive when all measured rates are positive (a zero measurement pins its groups to exactly zero).
Fixed point: a guess that already reproduces the rates is unchanged (every correction factor is 1).
Unseen groups keep the guess: a group no detector responds to has Σ_i W_ji = 0; its adjustment factor is exactly 1.

W_ji is the canonical base SAND-II weight. The original implementation’s optional per-detector statistics factor (N_i²/σ_i) is not applied, so SAND-II adjustments are not weighted by detector uncertainty^mcelroy-1967 — that weighting is what the GRAVEL sibling below adds.

T1–T4 — STAYSL-class damped least-squares

Where SAND-II multiplies the guess by correction factors, the STAYSL-class method adjusts it toward the least-squares fit of the measured rates. The forward fold is linear, so the objective is the quadratic anchored-damped form:

c_i = \sum_j R_{ij}\,\phi_j \tag{T1}

\phi^\dagger = \arg\min_\phi \sum_i w_i\,(c_i(\phi) - N_i)^2 + \lambda \sum_j (\phi_j - \phi_{\mathrm{guess},j})^2, \qquad w_i = 1/\sigma_i^2 \tag{T2}

Each iterator cycle re-solves (T2) about the current iterate as one weighted least-squares system through the shared linalg least-squares kernel — there is no second least-squares implementation — solving for the increment Δ = φ† − φ_k:

\Delta_A = \arg\min_\Delta \sum_i w_i \left(\sum_{j \in A} R_{ij}\,\Delta_j - r_i\right)^2 + \lambda \sum_{j \in A} (\Delta_j - (\phi_{\mathrm{guess},j} - \phi_{k,j}))^2, \qquad r_i = N_i - c_i(\phi_k) \tag{T3}

with groups the solve drives non-positive clipped to zero and pinned there for the rest of the run (T4, the active-set convention). Sigmas are caller-supplied, one finite positive value per detector; anything else raises a clear error. The damping λ is absolute Tikhonov damping on the distance to the guess: directions the data constrains much stronger than λ fit the rates, directions constrained much weaker keep the guess. Once a cycle produces no clip and no movement above the tolerance, the next cycle is a bit-for-bit no-op and the run has converged. Unlike SAND-II the update is additive rather than multiplicative, and the fixed point reproduces the rates only up to the damping slack.

G1–G3 — GRAVEL chi-square weights

GRAVEL is the uncertainty-aware sibling of SAND-II: the same positivity-preserving multiplicative adjustment, but each detector’s correction is weighted by its measurement precision:

c_i = \sum_j R_{ij}\,\phi_j \tag{G1}

W_{ji} = \frac{R_{ij}\,\phi_j}{c_i} \cdot \frac{N_i^2}{\sigma_i^2} \tag{G2}

\phi_j \leftarrow \phi_j \cdot \exp\!\left( \frac{\sum_i W_{ji}\,\ln(N_i / c_i)}{\sum_i W_{ji}} \right) \tag{G3}

(G2) is the SAND-II base rate share times the measurement-weight factor N_i²/σ_i²: a precisely measured detector dominates the groups it responds to, a sloppy one barely pulls. At a fixed point every N_i/c_i is 1, so the weights are irrelevant there — a SAND-II fixed point is a GRAVEL fixed point and vice versa. Two pinned divergences from SAND-II: zero measurements carry zero weight and are simply not fitted (they pin nothing to zero — an all-zero measurement set returns the guess unchanged after one confirming no-op), and sigmas are required input (zero, non-finite, or missing sigmas raise a clear error, never silent uniform weighting).

M1–M6 — MAXED maximum entropy

Where SAND-II, STAYSL-class, and GRAVEL spread each detector’s correction across groups or fit the rates directly, MAXED maximizes the Shannon relative entropy against the guess-as-prior subject to the chi-square fit — among the spectra that reproduce the rates it keeps the one closest to the guess in the information sense:

c_i = \sum_j R_{ij}\,\phi_j \tag{M1}

S(\phi) = \sum_j \left(\phi_j - \phi^0_j - \phi_j \ln(\phi_j / \phi^0_j)\right) \tag{M2}

\chi^2(\phi) = \sum_i \left(\frac{c_i - N_i}{\sigma_i}\right)^2 \tag{M3}

(M2) is non-positive with equality only at the guess φ⁰. Its Lagrange stationary form under rate constraints is the exponential family in the per-detector multipliers λ:

\phi_j(\lambda) = \phi^0_j \cdot \exp\!\left(-\sum_i R_{ij}\,\lambda_i\right) \tag{M4}

Each cycle linearizes the fold about the current multipliers (∂c_i/∂λ_k = −J_ik) and takes the dual step — one weighted least-squares solve through the shared linalg kernel over the fitted detectors only:

J_{ik} = \sum_j R_{ij}\,\phi_j\,R_{kj} \tag{M5}

\min_\delta \sum_{i} \frac{\left(\sum_k H_{ik}\,\delta_k - v_i\right)^2}{{\sigma'}_i^2} + \mu\,\|\Delta\|^2, \qquad \Delta_k = \delta_k / q_k \tag{M6}

(M6) equilibrates the solve (max-abs row/column scales, with the row scales compensated into the weights so the caller 1/σ_i² weights and the fixed points survive bit-for-bit) plus an adaptive Levenberg-Marquardt damping term on the physical multiplier step, and accepts only a strictly better fit that moves no group by more than a factor of e. For consistent rates the multipliers converge to the maximum-entropy exact-fit spectrum; for inconsistent (redundant disagreeing) readings they converge to the minimum-chi-square member. The full derivation lives in the nucleide-unfold crate rustdoc (maxed module). Like GRAVEL, zero measurements carry zero weight and are simply not fitted; unlike the other three methods the run also takes a caller chi-square target (default: the detector count) — converging in relative change to a fit whose chi-square still exceeds the target is a hard error, never a partial spectrum.

The IRDFF-II runtime response pack

Response rows need not be hand-built: nucleide.data.fetch_irdff downloads the official IAEA IRDFF-II_g725.zip distribution at runtime and checks its SHA-256 hash (a mismatch or offline run raises a clear error naming the URL or both digests), and parse_irdff_g725 parses the IRDFF-II.g725 member into caller-ready rows over the SAND-II 725-group structure for the v1 foil-activation subset (Au-197 and In-115 capture, U-235 and U-238 fission, Fe-56 and Ni-58 (n,p), Al-27 (n,α), Na-23 (n,2n)) plus the v2 extension of 26 further dosimetry sections (threshold (n,p)/(n,a)/(n,2n) monitors, capture foils, fission chambers, high-threshold bismuth monitors). Nothing IAEA-copyright is vendored — the zip is fetched, never committed — and the rows feed sandii, staysl, gravel, and maxed unchanged. The pointwise isomer reactions stay out of the pack.

Underdetermined systems

With fewer detectors than energy groups the system has infinitely many solutions. Groups no detector responds to (zero-response columns) keep the guess exactly, but otherwise the adjustment is nonlinear — its weights depend on the current spectrum — so a general unseen direction is not preserved the way a linear least-squares null space would be. In the U4 check (6 detectors, 24 groups) the rates are reproduced to 1e-6 and the total guess error is cut by more than half (measured at about 2.5x): the recoverable bias converges while the unseen bias keeps the guess.

Convergence and diagnostics

After each adjustment the largest per-group relative change max_j |φ_new − φ_old| / φ_old is compared against the tolerance; the run converges when it drops strictly below it. The defaults are 1e-3 and 200 iterations — library defaults, not values from the report. MAXED additionally compares the chi-square against its target (default: the detector count, the chi-square expectation): the run converges only when the change is below tolerance and the fit sits at or below the target. Each adjustment also reports the refolded rates and the per-detector rate factors N_i / c_i: all ones at a fixed point, and 0.0 for a pinned detector (a zero measurement folding to zero — the 0/0 factor spelled as 0.0, never NaN).

One SAND-II cycle (S1–S3): fold, weigh rate ratios by detector shares, adjust, repeat to tolerance.

Algorithm

Construction validates the whole input set once, naming the offending input on every rejection: the response matrix must be non-empty, rectangular, finite, and non-negative; rates finite and non-negative, one per detector; the guess strictly positive, one value per group; the tolerance finite and positive; the iteration cap at least 1.

Each iterator step folds the working spectrum (S1), applies the S2/S3 adjustment group by group, and yields the post-adjustment spectrum, refolded rates, rate factors, and the largest relative change. The weighting has named, deliberate zero-handling:

A detector folding to zero contributes no weight (it is skipped).
A zero response entry R_ij = 0 skips detector i for group j — this also avoids the 0 · ln 0 = NaN hazard.
With no active weight a group keeps the guess (factor 1, exactly).

Zero measurements pin their groups to zero (their adjustment factor is 0). A positive rate against a zero response row is rejected up front; if pinning zeroes out the entire support of another detector with a nonzero rate, the run stops mid-iteration with the same clear error naming that detector. The unfold driver iterates to convergence and returns the solution with its diagnostics — a halted or non-converged run raises a clear error, never a silent partial spectrum.

Worked example: one detector, three groups

A single detector responds to only one of three energy groups:

import nucleide.unfold as unfold

response = [[1.0, 0.0, 0.0]]  # one detector; it responds only to group 0
rates = [4.0]  # its measured rate
guess = [2.0, 7.0, 3.0]  # strictly positive starting spectrum

out = unfold.sandii(response, rates, guess)
out["spectrum"]  # [4.0, 7.0, 3.0]
out["iterations"]  # 2
unfold.forward_fold(response, out["spectrum"])  # [4.0]

The first adjustment, by hand: the fold is c = 1·2 + 0·7 + 0·3 = 2. Group 0 carries weight W = R·φ/c = 1·2/2 = 1 from detector 0 only, so its correction factor is N/c = 4/2 = 2 and its new value is 2 · exp(1 · ln 2 / 1) = 4. Groups 1 and 2 receive no weight, so their factor is exactly 1 and they keep the guess values 7.0 and 3.0 bit-for-bit. The largest relative change is |4 − 2| / 2 = 1.0, not below the default 1e-3 tolerance, so a second adjustment runs: the fold is now 4, the factor is exp(ln(4/4)) = 1, the change is 0.0, and the run converges with spectrum = [4.0, 7.0, 3.0] and rate_factors = [1.0].

This is the underdetermined behavior in miniature: one equation, three unknowns. The seen group moves to its measured value; the unseen groups keep the guess — SAND-II cannot and does not invent information no detector recorded.

Assumptions and limitations

Caller-supplied data, plus one runtime pack. Response matrix, measured rates, guess, and energy-group bounds are all inputs. Evaluated libraries (IRDFF and other IAEA-copyright data) are never vendored, and the IRDFF-II tabulated group spectra are never used — the benchmark shapes below are re-derived from published closed forms with citation^trkov-2020. The one exception to hand-built responses is the IRDFF-II runtime pack above: fetched and hash-pinned at runtime, never committed.
Strictly positive guess required; response entries and rates must be non-negative and finite. Every malformed input raises a clear error naming the cause.
Detector-uncertainty weighting is the caller’s choice of method: SAND-II applies no per-detector statistics factor (see S3); STAYSL-class, GRAVEL, and MAXED consume caller-supplied sigmas as weights (see T2/G2/M6: 1/σ², N²/σ², and 1/σ² respectively).
Underdetermination is real. With fewer detectors than groups the result is not unique: the unseen component of the guess survives unchanged, so the unfolded spectrum depends on the guess there. Iteration reduces the error; it does not remove it.
Zero measurements pin their groups to zero under SAND-II, and a positive rate whose support is pinned away fails with a clear error, up front or mid-iteration. GRAVEL and MAXED instead skip zero measurements (zero weight, never fitted).
Non-convergence is a hard error: exhausting the iteration cap raises (a ValueError from Python), never returning a partial spectrum — and for MAXED, so does converging in relative change to a fit whose chi-square still exceeds the target.
No outlier-foil rejection (the original code’s σ-discard feature) and no cover-material corrections — recorded as out of scope.
Four methods, one at a time: SAND-II, STAYSL-class, GRAVEL, and MAXED share the iterator and diagnostics shape but run separately — there is no hybrid that mixes their updates within one run. MAXED closes the recorded one-method-per-cycle order; no fifth method is scoped.

Validation

The crate and validation/unfold_vs_analytic.py replay twenty-three pinned checks, U1–U23, all on synthetic forward-fold-then-recover round-trips — the response matrices are hand-built Gaussian bumps over log-spaced energy groups, never evaluated data:

U1: the forward fold matches a hand-written matrix-vector product to 1e-12.
U2: an exact guess is a fixed point — one confirming no-op adjustment with a change of exactly 0.0.
U3: a determined 6-detector, 6-group system is recovered from a uniformly 3x-biased guess to 1e-6.
U4: the underdetermined case — 6 detectors, 24 groups. Rates reproduced to 1e-6; spectrum error cut by more than half, with the recoverable bias converging and the unseen bias keeping the guess.
U5: the IRDFF-II analytical benchmark-field shapes named in Trkov et al.^trkov-2020 — a thermal Maxwellian at 293.6 K, φ(E) ∝ E·e^(−E/kT) with kT = 8.617333262·10⁻⁵ eV/K × 293.6 K ≈ 0.0253 eV; a pure 1/E field, φ(E) ∝ 1/E; and a Maxwellian at 25 keV, φ(E) ∝ √E·e^(−E/kT) with kT = 25 keV — re-derived from the published closed forms and recovered to 1e-5 in a determined probe. The IAEA-copyright tabulated spectra are deliberately not used.
U6: exhausting the iteration cap raises a clear error.
U7: unconstrained groups keep the guess bit-for-bit.
U8: malformed inputs raise clear errors naming the cause.
U9–U13 (STAYSL-class): caller sigmas weight the detectors (a redundant reading with 100x the sigma loses to the precise one); a determined 6x6 system is recovered from a 3x-biased guess to 1e-6; the underdetermined case reproduces rates; the IRDFF-II shapes recover to 1e-5; a SAND-II fixed point is a least-squares fixed point; and the named-contract errors (non-convergence, sigma shape, zero/weak responses) raise in four probes.
U14–U18 (GRAVEL): the same shape — sigma weighting with precise-wins, determined recovery to 1e-6, underdetermined rate reproduction to 1e-6 with spectrum error reduced under counting-statistics sigmas, IRDFF-II shape recovery to 1e-5, the cross-method fixed-point check in both directions (at a fixed point every rate ratio is 1, so the weights are irrelevant), and five named-contract error probes.
U19–U23 (MAXED): the same shape — 1/σ² chi-square weighting with precise-wins, determined recovery to 1e-6, underdetermined rate reproduction to 1e-6 with spectrum error reduced under counting-statistics sigmas, IRDFF-II shape recovery to 1e-5 with group bounds bracketing each shape’s support (the solve-based house rule), the cross-method fixed-point check in both directions (at a fixed point the dual residual is zero), and seven named-contract error probes — including the chi-square-target hard fail on disagreeing redundant readings.

No external reference implementation exists — PyNE and OpenMC ship no SAND-II/STAYSL/GRAVEL/MAXED counterpart — so these analytic checks are the whole correctness story: they always run, and none can be skipped.

No fifth unfolding method is scoped after MAXED, which closes the recorded one-method-per-cycle order. Because no external unfolding counterpart exists in PyNE or OpenMC, cross-code comparison is not available for this capability.