LOUHI78 unfolding: constrained weighted least squares with generalized smoothing#

The unfold_louhi method is the Python port of the LOUHI78 general purpose unfolding program of J. T. Routti and V. Sandberg (General purpose unfolding program LOUHI78 with linear and nonlinear regressions, Computer Physics Communications 21 (1980) 119-135, doi:10.1016/0010-4655(80)90021-4). LOUHI is one of the classical Bonner-sphere unfolding codes: it formulates the unfolding problem as a constrained weighted least-squares fit with a generalized smoothing term and determines the solution by quadratic programming.

Mathematical formulation#

Given the response matrix \(A\), the measurement vector \(b\), the per-detector uncertainties \(\sigma_i\) and the default (a-priori) spectrum \(\phi_0\), LOUHI minimizes the total chi-square functional

\[\chi^2(\phi) \;=\; \sum_i \left( \frac{b_i - (A\phi)_i}{\sigma_i} \right)^2 \;+\; \lambda^2 \, \big\| L (\phi - \phi_0) \big\|_2^2\]

subject to the physical non-negativity constraints \(\phi_j \ge 0\). The smoothing operator \(L\) supports three orders:

  • smooth_order=0 — identity: the solution shrinks toward the a-priori spectrum;

  • smooth_order=1 (default) — first differences of the deviation from the a-priori spectrum;

  • smooth_order=2 — second differences (curvature) of the deviation.

Assembling the normal equations turns the problem into the quadratic program

\[\min_{\phi \ge 0} \; \tfrac12 \phi^T H \phi - g^T \phi, \qquad H = 2\left( A^T W A + \lambda^2 L^T L \right), \quad g = 2\left( A^T W b + \lambda^2 L^T L \phi_0 \right),\]

with \(W = \operatorname{diag}(1/\sigma_i^2)\), which is solved by Hildreth’s iterative coordinate algorithm (the LSI step of LOUHI78): every sweep minimizes one coordinate exactly and projects it onto its non-negativity constraint, \(\phi_j \leftarrow \max(0,\, \phi_j - \nabla_j / H_{jj})\). The sweep stops when the relative change of the quadratic objective drops below tolerance.

Linear and nonlinear modes#

In the linear mode the smoothing weight smoothness (\(\lambda\)) is fixed by the user. In the nonlinear mode (auto_smooth=True) LOUHI adjusts \(\lambda\) automatically by a nonlinear regression: a golden-section search on \(\log_{10}\lambda\) drives the data chi-square to its expected value \(\chi^2_{\text{target}}\) (default: the number of detectors, the expectation of the chi-square distribution). This reproduces the automatic smoothing-parameter search of the original program.

Error propagation#

The statistical error report of LOUHI78 propagates the measurement uncertainties through the inverse of the Hessian restricted to the free (strictly positive) bins of the active set; constrained bins at zero carry no variance in the first-order propagation. The helper bssunfold.core.unfold_louhi.louhi_covariance() returns the resulting per-bin standard deviations.

Usage#

from bssunfold import Detector

detector = Detector()
result = detector.unfold_louhi(
    readings={"3in": 0.053, "5in": 0.184, "10in": 0.172, "18in": 0.034},
    smoothness=1.0,        # fixed smoothing weight (linear mode)
    smooth_order=1,        # first-difference smoothing operator
)

# nonlinear mode: automatic smoothing-weight regression
result = detector.unfold_louhi(readings, auto_smooth=True)

# explicit a-priori spectrum (anchored smoothing)
import numpy as np
phi0 = ...  # physically informed default spectrum
result = detector.unfold_louhi(readings, initial_spectrum=phi0)

API reference#

See also Detector Class for the full Detector API reference.

bssunfold.core.unfold_louhi.unfold_louhi(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, smoothness: float = 1.0, smooth_order: int = 1, auto_smooth: bool = False, chi2_target: float | None = None, max_iterations: int = 500, tolerance: float = 1e-06, relative_uncertainty: float = 0.1, calculate_errors: bool = False, noise_level: float = 0.01, n_montecarlo: int = 100, variance_reduction: str = 'none', save_result: bool = False, random_state: int | None = None) → dict[str, Any][source]

Unfold neutron spectrum using the LOUHI78 algorithm.

Parameters:

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

  • n_energy_bins (int) – Number of energy bins.

  • E_MeV (np.ndarray) – Energy grid.

  • 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) – Default (a-priori) spectrum. If None, a flat spectrum is used.

  • smoothness (float, optional) – Smoothing weight lambda (default: 1.0).

  • smooth_order (int, optional) – Smoothing operator order 0/1/2 (default: 1, first differences).

  • auto_smooth (bool, optional) – Adjust the smoothing weight automatically to reach chi2_target (default: False).

  • chi2_target (Optional[float], optional) – Target data chi-square for auto_smooth (default: number of detectors).

  • max_iterations (int, optional) – Maximum number of Hildreth sweeps (default: 500).

  • tolerance (float, optional) – Relative objective change per sweep for convergence (default: 1e-6).

  • relative_uncertainty (float, optional) – Relative measurement uncertainty (default: 0.1).

  • 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).

  • variance_reduction (str, optional) – MC variance reduction: ‘none’, ‘antithetic’, ‘control’, ‘both’ (default: ‘none’).

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

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

Returns:

Unfolding results dictionary.

Return type:

Dict[str, Any]

bssunfold.core.unfold_louhi.solve_louhi(A: ndarray, b: ndarray, x0: ndarray | None = None, smoothness: float = 1.0, smooth_order: int = 1, max_iterations: int = 500, tolerance: float = 1e-06, relative_uncertainty: float = 0.1, sigma: ndarray | None = None, auto_smooth: bool = False, chi2_target: float | None = None) → tuple[ndarray, int, bool][source]

Solve the unfolding problem with the LOUHI78 algorithm.

Parameters:

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

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

  • x0 (np.ndarray, optional) – Default (a-priori) spectrum (n,) used as the starting point and as the reference of the generalized smoothing term. When None, a flat unit spectrum is used.

  • smoothness (float, optional) – Smoothing weight lambda of the generalized smoothing term (default: 1.0). Ignored when auto_smooth is True.

  • smooth_order (int, optional) – Order of the smoothing operator: 0 (identity), 1 (first differences) or 2 (second differences); default 1.

  • max_iterations (int, optional) – Maximum number of Hildreth coordinate sweeps (default: 500).

  • tolerance (float, optional) – Maximum relative change between sweeps for convergence (default: 1e-6).

  • relative_uncertainty (float, optional) – Relative measurement uncertainty used to derive detector sigma values when sigma is not supplied (default: 0.1).

  • sigma (np.ndarray, optional) – Explicit per-detector measurement uncertainties (m,). When given, overrides relative_uncertainty.

  • auto_smooth (bool, optional) – Nonlinear regression mode of LOUHI78: adjust the smoothing weight by a golden-section search on log10(lambda) so that the data chi-square reaches chi2_target (default False).

  • chi2_target (float, optional) – Target data chi-square for auto_smooth. Defaults to the number of detectors (the expected value of the chi-square).

Returns:

(solution spectrum, sweeps used, converged flag).

Return type:

Tuple[np.ndarray, int, bool]

bssunfold.core.unfold_louhi.louhi_smoothing_matrix(n: int, smooth_order: int = 1) → ndarray[source]

Build the generalized smoothing operator L of LOUHI.

Parameters:

  • n (int) – Number of energy bins.

  • smooth_order (int, optional) – Order of the smoothing functional: 0 shrinks the solution toward the default spectrum (identity operator), 1 penalizes first differences of the deviation from the default spectrum and 2 penalizes second differences (default: 1).

Returns:

The (n, n) smoothing matrix L.

Return type:

np.ndarray

bssunfold.core.unfold_louhi.louhi_covariance(A: ndarray, sigma: ndarray, smoothness: float, smooth_order: int, x0: ndarray, x: ndarray) → ndarray[source]

Propagate measurement uncertainties for the LOUHI solution.

Follows the statistical error analysis of LOUHI78: the covariance of the free (strictly positive) spectrum bins is the inverse of the Hessian restricted to the active set, while constrained bins at zero carry no variance in the first-order propagation.

Parameters:

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

  • sigma (np.ndarray) – Per-detector measurement uncertainties (m,).

  • smoothness (float) – Smoothing weight used for the solution.

  • smooth_order (int) – Order of the smoothing operator used for the solution.

  • x0 (np.ndarray) – Default spectrum used in the smoothing term.

  • x (np.ndarray) – LOUHI solution spectrum.

Returns:

Standard deviations of the spectrum bins (n,).

Return type:

np.ndarray