N-spline unfolding (Islamgulov & Lartsev, 2008)#

The unfold_nspline method implements the neutron spectrum unfolding approach of

R. F. Islamgulov, V. D. Lartsev, Reconstruction of neutron spectra from activation measurements in the form of N-splines, Atomic Energy 104(5), 295–302 (May 2008) (RFNC — VNIITF named after E. I. Zababakhin).

Method outline#

The activation measurements are described by the integral system

\[Q_i = \int_0^\infty \sigma_i(E)\,\varphi(E)\,dE, \qquad i = 1,\dots,N, \qquad (1)\]

which is a classical ill-posed inverse problem. Instead of solving for \(n\) bin values of \(\varphi(E)\) directly, the paper parameterises the spectrum by a specialised neutron spline (N-spline) with basis functions

\[N_k(E) = \exp(a_k + q_k \ln E + r_k E), \qquad E_k \le E \le E_{k+1},\; k = 1,\dots,M, \qquad (2)\]

i.e. piecewise functions whose logarithm is linear both in \(\ln E\) and in \(E\). This family contains the classical model spectra — \(1/E\), Maxwellian evaporation \(\exp(-E/T)\), fission-like \(\sqrt{E}\exp(-bE)\), two-component Maxwell + slowed-down representations — as particular members, so the basis is close to complete for reactor and accelerator spectra while using only \(3M\) parameters.

Continuity of the spline value and of its derivative at the interior knots (Eqs. 3–4 of the paper),

\[ \begin{align}\begin{aligned}a_k - a_{k+1} + u_k (q_k - q_{k+1}) + E_k (r_k - r_{k+1}) = 0, \qquad u_k = \ln E_k,\\(q_k - q_{k+1}) + E_k (r_k - r_{k+1}) = 0,\end{aligned}\end{align} \]

is assembled into the block matrix \(D\) (Eq. 5), \(X = (a, q, r)^T\), giving the constraint system \(DX = 0\) (build_continuity_matrix()). The pointwise approximation of a tabulated spectrum (Eqs. 6–7) is a weighted least-squares problem in the log domain with weights \(w_{ki} = 1/\varepsilon_{ki}\) subject to \(DX = 0\):

\[GX = Y, \qquad DX = 0,\]

solved via the KKT (Lagrange multiplier) system (fit_nspline()).

The activation equations are then solved by the directed divergence minimisation loop (generalised MIRD algorithm of Lartsev, Preprint RFNC-VNIITF No. 216, 2005; Tarasko, Preprint FEI No. 1446, 1983). With normalised measured activations \(p_i = Q_i/\sum_j Q_j\) the functional

\[H = \sum_i \left[ pN_i \ln \frac{pN_i}{p_i} - pN_i + p_i \right] \ge 0\]

(measured \(p_i\) vs calculated \(pN_i\)) is decreased by the flux-conserving gradient iteration

\[ \begin{align}\begin{aligned}\varphi_{n+1}(E) = \varphi_n(E)\,\bigl[1 - \Delta\mu_n (R_n(E) - \bar R_n)\bigr],\\R_n(E) = \sum_i \frac{p_i}{Q_i}\,\sigma_i(E) \ln\frac{pN_i}{p_i},\end{aligned}\end{align} \]

where \(\bar R_n\) is the flux-weighted mean of \(R_n\) and the step \(\Delta\mu_n\) starts from the paper’s conservative value \(0.1/\sup|R_n - \bar R_n|\) and is halved while \(H\) increases (backtracking). After every iteration the current spectrum is re-fitted by the N-spline — the paper’s regularisation, which makes the loop act on \(3M\) spline parameters instead of \(n\) bin values and avoids the nonlinearity / local-minimum difficulties of a direct spline fit to the activation integrals.

Stopping criteria and quality control#

Iterations stop when \(H\) reaches the level corresponding to the measurement errors,

\[H \le H_{\text{target}} = \tfrac12 \sum_i p_i \left(\frac{\Delta Q_i}{Q_i}\right)^2,\]

or when the relative decrease of \(H\) per iteration falls below tol. The acceptability of the reconstruction is measured by the paper’s mean-squared residual

\[nev = \sqrt{\frac{1}{N-1} \sum_i \left(\frac{Qr_i - Q_i} {\Delta Q_i}\right)^2}, \qquad \text{acceptable when } nev \le 1 + \frac{2}{\sqrt{N}},\]

with \(Qr_i\) the activation integrals recalculated from the unfolded spectrum.

Knots#

Knot sets used in the paper for the BARS-5, IGRIK (channel and surface) and YAGUAR reactors are available as NSPLINE_KNOT_PRESETS ("BARS5_channel", "IGRIK_channel", "IGRIK_surface", "YAGUAR_channel"). knots=None (default) builds a log-uniform grid, and explicit knot sequences (MeV) are accepted as well. User and preset knots are clipped to the energy grid range and the outer knots are extended so that the spline domain always spans the whole grid, as in the paper.

API#

bssunfold.core.unfold_nspline.unfold_nspline(detector_names: list[str], n_energy_bins: int, E_MeV: ndarray, sensitivities: dict[str, ndarray], cc_icrp116: dict[str, ndarray], save_result_callback, readings: dict[str, float], initial_spectrum: ndarray | None = None, knots: str | Sequence[float] | None = None, continuity: str = 'C0C1', relative_uncertainty: float = 0.1, max_iterations: int = 200, tol: float = 0.001, step_theta: float = 0.1, smoothing: bool = True, n_segments: int | None = None, calculate_errors: bool = False, noise_level: float = 0.01, n_montecarlo: int = 100, save_result: bool = False, random_state: int | None = None) → dict[str, Any][source]#

Unfold neutron spectrum using the N-spline method (2008).

Detector-facing wrapper of the Islamgulov & Lartsev N-spline / directed-divergence unfolding algorithm (Atomic Energy 104(5), 2008).

Parameters:

  • detector_names (List[str]) – Names of available detectors.

  • n_energy_bins (int) – Number of energy bins.

  • E_MeV (np.ndarray) – Energy grid (MeV).

  • sensitivities (Dict[str, np.ndarray]) – Detector sensitivity arrays.

  • cc_icrp116 (Dict[str, np.ndarray]) – ICRP-116 conversion coefficients.

  • save_result_callback (callable) – Callback to save result to history.

  • readings (Dict[str, float]) – Detector readings.

  • initial_spectrum (Optional[np.ndarray], optional) – Initial spectrum guess (the paper recommends a Monte-Carlo calculated spectrum); flat when None.

  • knots (str / Sequence[float] / None, optional) – Knot preset name (NSPLINE_KNOT_PRESETS: “BARS5_channel”, “IGRIK_channel”, “IGRIK_surface”, “YAGUAR_channel”), explicit knot sequence or None for an automatic log-uniform grid.

  • continuity (str, optional) – Spline continuity: "C0C1" (default), "C0" or "none".

  • relative_uncertainty (float, optional) – Relative measurement uncertainty dQ/Q used in the stopping criteria and the nev statistic (default: 0.1).

  • max_iterations (int, optional) – Iteration budget (default: 200).

  • tol (float, optional) – Relative H-decrease stopping tolerance (default: 1e-3).

  • step_theta (float, optional) – Conservative step factor (default: 0.1).

  • smoothing (bool, optional) – Per-iteration N-spline smoothing (default: True).

  • n_segments (int, optional) – Number of spline segments when knots=None.

  • calculate_errors (bool, optional) – Calculate Monte-Carlo errors (default: False).

  • noise_level (float, optional) – Noise level for Monte-Carlo (default: 0.01).

  • n_montecarlo (int, optional) – Number of Monte-Carlo samples (default: 100).

  • save_result (bool, optional) – Save result to history (default: False).

  • random_state (int, optional) – Random seed for reproducibility.

Returns:

Standardized unfolding results dictionary enriched with the method diagnostics (H, H_history, H_target, nev, nev_limit, acceptable, stop_reason, fluence, mean_energy, knots, knots_source).

Return type:

Dict[str, Any]

bssunfold.core.unfold_nspline.solve_nspline(A: ndarray, b: ndarray, x0: ndarray | None = None, E_MeV: ndarray | None = None, knots: str | Sequence[float] | None = None, sigma_rel: ndarray | None = None, continuity: str = 'C0C1', max_iterations: int = 200, tol: float = 0.001, step_theta: float = 0.1, smoothing: bool = True, n_segments: int | None = None) → tuple[ndarray, int, bool][source]#

Solve the unfolding problem using the N-spline method.

Thin standard-API wrapper around solve_nspline_full() returning the usual (spectrum, iterations, converged) tuple used by the bssunfold solver protocol.

Parameters:

  • A (np.ndarray) – Response matrix (m, n).

  • b (np.ndarray) – Measurement vector (m,).

  • x0 (np.ndarray, optional) – Initial spectrum guess (n,) (flat when None).

  • E_MeV (np.ndarray, optional) – Energy grid (MeV), positive; required.

  • knots (str / Sequence[float] / None, optional) – Knot preset name, explicit knots or None (auto grid).

  • sigma_rel (np.ndarray, optional) – Relative measurement uncertainties (default: 0.1).

  • continuity (str, optional) – Spline continuity: "C0C1" (default), "C0" or "none".

  • max_iterations (int, optional) – Iteration budget (default: 200).

  • tol (float, optional) – Relative H-decrease stopping tolerance (default: 1e-3).

  • step_theta (float, optional) – Conservative step factor (default: 0.1, the paper’s value).

  • smoothing (bool, optional) – Per-iteration N-spline smoothing (default: True).

  • n_segments (int, optional) – Number of segments when knots=None.

Returns:

(spectrum, iterations, converged).

Return type:

Tuple[np.ndarray, int, bool]

bssunfold.core.unfold_nspline.solve_nspline_full(A: ndarray, b: ndarray, x0: ndarray | None = None, E_MeV: ndarray | None = None, knots: str | Sequence[float] | None = None, sigma_rel: ndarray | None = None, continuity: str = 'C0C1', max_iterations: int = 200, tol: float = 0.001, step_theta: float = 0.1, smoothing: bool = True, n_segments: int | None = None) → dict[str, Any][source]#

Full N-spline unfolding with diagnostics (Islamgulov & Lartsev, 2008).

Iteratively minimises the directed divergence H between the measured and calculated normalised activations, smoothing the spectrum through an N-spline fit at every iteration (the paper’s regularisation). Uses the paper’s stopping criteria (H at the measurement-error level or stalled relative decrease) and reports the nev residual statistic with the acceptance bound nev <= 1 + 2/sqrt(N).

Parameters:

  • A (np.ndarray) – Response matrix (m, n) – activation responses of the detectors on the energy grid.

  • b (np.ndarray) – Measured readings / activation integrals (m,).

  • x0 (np.ndarray, optional) – Initial spectrum guess (n,). None means a flat spectrum. In the spirit of the paper this should be a Monte-Carlo calculated spectrum when available.

  • E_MeV (np.ndarray, optional) – Energy grid (MeV). Required; must be positive.

  • knots (str / Sequence[float] / None, optional) – Knot preset name (see NSPLINE_KNOT_PRESETS), explicit knot sequence or None (automatic log-uniform grid).

  • sigma_rel (np.ndarray, optional) – Relative measurement uncertainties dQ_i/Q_i (m,) used by the stopping criteria and nev (default: 0.1 for every detector).

  • continuity (str, optional) – "C0C1" (default), "C0" or "none".

  • max_iterations (int, optional) – Iteration budget (default: 200).

  • tol (float, optional) – Relative-decrease stopping tolerance for H (default: 1e-3).

  • step_theta (float, optional) – Conservative initial step factor: dmu = step_theta / sup|R-Rbar| (the paper’s 0.1); backtracking halves it while H increases.

  • smoothing (bool, optional) – Re-fit the N-spline after every iteration (default: True, the paper’s procedure; False reduces to the plain MIRD loop).

  • n_segments (int, optional) – Number of spline segments when knots=None.

Returns:

Diagnostics dictionary with keys spectrum, iterations, converged, stop_reason, H, H_history,

H_target, nev, nev_limit, acceptable, Qr, relative_residuals, fluence, mean_energy, knots, knots_source, continuity and params.

Return type:

Dict[str, Any]

bssunfold.core.unfold_nspline.fit_nspline(E: ndarray, phi: ndarray, knots: str | Sequence[float] | None = None, rel_err: ndarray | None = None, continuity: str = 'C0C1', n_segments: int | None = None) → tuple[ndarray, dict[str, Any]][source]#

Approximate a pointwise spectrum by an N-spline (Eqs. 2, 5-7).

Solves the weighted log-domain least-squares problem with continuity constraints of the paper:

min_X sum_j w_j^2 (a_kj + u_j q_kj + E_j r_kj - ln phi_j)^2 s.t. D X = 0, w_j = 1 / eps_j,

via the KKT (Lagrange multiplier) system

[[G^T W G, D^T], [D, 0]] [X; lam] = [G^T W Y; 0].

Parameters:

  • E (np.ndarray) – Energy grid (MeV), positive values.

  • phi (np.ndarray) – Spectrum values on the grid (non-negative; zeros are floored and down-weighted).

  • knots (str / Sequence[float] / None, optional) – Knot preset name, explicit knots or None for an automatic log-uniform grid (default).

  • rel_err (np.ndarray, optional) – Relative pointwise errors eps_j; weights w_j = 1/eps_j (paper’s Eq. 7). None means unit weights.

  • continuity (str, optional) – "C0C1" (default), "C0" or "none" – see build_continuity_matrix().

  • n_segments (int, optional) – Number of segments when knots=None (default: adaptive).

Returns:

(N_E, info) where N_E is the fitted spline on E and info carries knots, a/q/r parameter arrays, the weighted RMS log-residual and the continuity option.

Return type:

Tuple[np.ndarray, Dict[str, Any]]

bssunfold.core.unfold_nspline.build_continuity_matrix(knots: Sequence[float], continuity: str = 'C0C1') → ndarray[source]#

Build the spline continuity matrix D of Eq. (5).

The N-spline parameter vector is X = (a, q, r)^T with a = (a_1..a_M), q = (q_1..q_M), r = (r_1..r_M). The continuity conditions at the interior knots read (Eqs. 3-4):

C0: a_k - a_{k+1} + u (q_k - q_{k+1}) + E (r_k - r_{k+1}) = 0 C1: (q_k - q_{k+1}) + E (r_k - r_{k+1}) = 0,

with u = ln E and E the knot value, and are assembled into

D = [[A, B, C], [0, A, C]], D X = 0.

Parameters:

  • knots (Sequence[float]) – Knot sequence (M = len(knots) - 1 segments).

  • continuity (str, optional) – "C0C1" (default) – continuous value and derivative; "C0" – continuous value only (only the first block row of D is kept, as noted in the paper); "none" – no continuity.

Returns:

Matrix of shape (rows, 3M) where rows is 0, M-1 or 2(M-1) depending on continuity.

Return type:

np.ndarray

bssunfold.core.unfold_nspline.nspline_eval(E: ndarray, a: Sequence[float], q: Sequence[float], r: Sequence[float], knots: Sequence[float]) → ndarray[source]#

Evaluate the N-spline N(E) = exp(a_k + q_k ln E + r_k E).

Parameters:

  • E (np.ndarray) – Evaluation energies (MeV), must be positive.

  • a (Sequence[float]) – Length-M parameter vectors.

  • q (Sequence[float]) – Length-M parameter vectors.

  • r (Sequence[float]) – Length-M parameter vectors.

  • knots (Sequence[float]) – The M+1 knot values.

Returns:

N-spline values (positive by construction).

Return type:

np.ndarray

bssunfold.core.unfold_nspline.auto_knots(E_MeV: ndarray, n_segments: int = 12) → tuple[float, ...][source]#

Build a log-uniform knot grid spanning the energy range of E_MeV.

Parameters:

  • E_MeV (np.ndarray) – Energy grid (MeV), must contain positive values.

  • n_segments (int, optional) – Number of spline segments (default: 12); the returned tuple has n_segments + 1 knots from min(E) to max(E).

Returns:

Strictly increasing knot sequence (MeV).

Return type:

Tuple[float, …]

Usage#

from bssunfold import Detector

det = Detector()
result = det.unfold_nspline(
    readings,
    knots="BARS5_channel",     # or None (auto), or explicit knots
    continuity="C0C1",         # spline continuity option
    relative_uncertainty=0.05,
    max_iterations=300,
)
print(result["nev"], result["acceptable"])   # paper's statistic
print(result["H_history"])                   # convergence trace

See examples/41-nspline.ipynb for a worked comparison with GRAVEL on a synthetic spectrum, and examples/42-nspline-iaea.ipynb for unfolding an IAEA Compendium Monte-Carlo BSA spectrum (t4-14-s.txt_1) from GSF Bonner-sphere readings.

Notes and limitations#

  • nev uses unweighted relative residuals: with readings spanning many orders of magnitude the weakest detector dominates the statistic (its contribution to \(H\) is negligible because \(H\) works with normalised activations). In the paper’s activation-foil applications all \(Q_i\) are comparable and \(nev \approx 1\)–2.

  • Energy regions where all detector responses vanish are shaped by the C0/C1 spline continuation (prior-driven); the pointwise smoothing weights follow the paper’s \(w = 1/\varepsilon\) philosophy with sensitivity-graded errors.

  • The method is verified by tests/test_nspline.py (32 tests): exact interpolation properties, continuity constraints, synthetic-spectrum recovery, the nev statistic and the Detector integration.