BSSUnfold

Read full text
Authors: Konstantin Chizhov
Cited by: 0

Abstract

Added License-required commercial QP solver methods — five new unfolding methods unfold_gurobi, unfold_mosek, unfold_cplex, unfold_copt and unfold_xpress (plus solve_gurobi, solve_mosek, solve_cplex, solve_copt, solve_xpress), driving the proprietary engines of Gurobi, MOSEK, IBM CPLEX, Cardinal Optimizer (COPT) and FICO Xpress through their cvxpy interfaces. Each method solves the same Tikhonov-regularized non-negative least-squares QP as the open-source path (0.5*||A x - b||^2 plus an identity or derivative penalty, L1/L2 norm options, Monte-Carlo uncertainty, max_neutron_energy support). License required: bssunfold never ships or checks a license; the engine package (gurobipy, mosek, cplex, coptpy, xpress) must be installed and licensed on the machine. Unlike unfold_cvxpy there is no open-source fallback — a missing engine or license warns and yields a zero spectrum, so results always come from the licensed engine that was requested. New general entry points solve_commercial / unfold_commercial taking a solver= alias, and metadata helpers commercial_solver_info / is_commercial_solver_available; the five public wrappers are pre-bound partials of these. New file core/unfold_commercial.py with the shared QP backend core/_commercial_qp.py (cvxpy quad_form build + solve, kwargs retry, non-optimal/license-error handling). New Detector methods unfold_gurobi, unfold_mosek, unfold_cplex, unfold_copt, unfold_xpress (full canonical signatures); all names registered in core/__init__.py. New optional extras in pyproject.toml: [gurobi], [mosek], [cplex], [copt], [xpress] and the aggregate [commercial] — deliberately not included in [all] because they are license-required proprietary packages (pip install "bssunfold[commercial]" after obtaining licenses). platform_check: new check_commercial_solvers_availability() (probes cvxpy installed_solvers() for the five engines, cached in COMMERCIAL_SOLVERS_AVAILABLE) and per-solver entries in get_available_solvers(). New test file tests/test_commercial.py (63 tests): solver metadata and alias validation, graceful zero-spectrum + warning behaviour for unavailable engines (each missing engine exercised via block_import), ill-formed input and norm guards (norm=1 requires nonneg), end-to-end numerics with the licensed CPLEX engine where installed (L2/L1/smoothness/ub/warm start, agreement with a tight-tolerance qpsolvers reference, objective equality), solver failure paths (mocked Problem.solve errors and non-optimal status), Monte-Carlo/save_result/regularization_method integration and max_neutron_energy masking. The commercial methods are wired into tests/test_all_unfold_methods_api.py with auto-skip when the engine package is not installed (CPLEX runs for real in the dev environment). OSEM-ANLM unfolding method (unfold_osem_anlm / solve_osem_anlm): port of the OSEM-ANLM algorithm of Jamaati et al. (2026), "Enhanced sparse view CT reconstruction using ordered subset expectation maximization and asymptotic non-local means algorithms", Scientific Reports (https://doi.org/10.1038/s41598-026-70607-1), adapted to neutron spectrum unfolding. Ordered-subset EM updates interleaved with the two-stage asymptotic non-local means (ANLM) filter applied after every subset update (article pseudo-code steps 4-5), or once to the OSEM result (anlm_mode='post'). The ANLM filter implements the article's two-stage scheme: stage 1 with the uniform parameter h1 = 0.5 * sigma, stage 2 with the point-wise parameter h2(i) = sqrt(sum_j w(i,j)^2 * sigma^2) (article eq. 6); the NLM weights use the Gaussian-weighted patch distances over the search window N (default 11) and similarity window nu (default 3, article optima). The noise level sigma is either user-supplied (h) or estimated automatically with a robust MAD estimator on second differences (estimate_noise_1d); by default the filter operates in log space (log_space=True), making it scale-free for spectra spanning orders of magnitude (log_space=False reproduces the raw-unit CT formulation). New file core/unfold_osem_anlm.py with the standalone 1D filter anlm_filter_1d. New Detector.unfold_osem_anlm() method; solve_osem_anlm, unfold_osem_anlm, anlm_filter_1d and estimate_noise_1d registered in core/__init__.py. New test file tests/test_osem_anlm.py (38 tests): filter identity and smoothing properties, noise estimator, solver equivalence with plain OSEM for the identity filter, post-mode composition solve_osem_anlm(post) == anlm_filter_1d(solve_osem(...)), subset variants, validation errors and Detector wrapper coverage. LOUHI78 unfolding (unfold_louhi / solve_louhi, Routti & Sandberg 1980, Computer Physics Communications 21, doi:10.1016/0010-4655(80)90021-4): constrained weighted least squares with generalized smoothing — the classic Bonner-sphere unfolding program, ported as min ||(b - A phi)/sigma||^2 + lambda^2 ||L (phi - phi0)||^2 s.t. phi >= 0: The quadratic program is solved by Hildreth's iterative coordinate algorithm (the LSI step of LOUHI78) with a relative-objective-change stopping rule; the default spectrum phi0 anchors both the smoothing term and the starting point. Generalized smoothing operators (new file core/unfold_louhi.py, helper louhi_smoothing_matrix): order 0 shrinks the solution toward the a-priori spectrum, orders 1/2 penalize first/second differences of the deviation from the a-priori. Nonlinear regression mode (auto_smooth=True): the smoothing weight is adjusted automatically by a golden-section search on log10(lambda) (reusing the core/_line_search.py building blocks) until the data chi-square reaches chi2_target (default: the number of detectors, i.e. the expected chi-square value). Statistical error propagation (louhi_covariance): inverse Hessian on the free (strictly positive) bins of the active set, as in the LOUHI78 error report. New Detector.unfold_louhi() method; registered in core/__init__.py. New tests in tests/test_louhi.py (42 tests covering the smoothing operators, the Hildreth QP core (interior-point equivalence and projection behavior), the automatic smoothing regression, the error propagation, parameter validation and the Detector integration). New example notebook examples/79-louhi-iaea.ipynb (LOUHI unfolding of an IAEA Compendium benchmark spectrum with the linear and nonlinear smoothing modes) and examples/80-osem-anlm-iaea.ipynb (OSEM-ANLM in the subset/post modes, the exact post == OSEM + anlm_filter_1d composition property and the h smoothing-strength knob); every public unfold_* method now has a dedicated IAEA benchmark notebook. New example notebook examples/81-cuqi-quality-analysis-iaea.ipynb (CUQIpy quality analysis, optional bssunfold[cuqi] backend): hierarchical Hybrid Gibbs unfolding with a data-driven flat prior center, ESS / Gelman-Rubin convergence diagnostics, 95% HPD credible intervals and coverage, residual analysis, EURADOS-style energy-group fluences and a comparison against CVXPY / MLEM / Landweber / GRAVEL. IAEA-Compendium benchmark notebook suite — 23 new example notebooks (examples/56-gnowee-iaea.ipynb through examples/78-ensemble-iaea.ipynb) for every public unfold_* method that previously lacked a dedicated example (gnowee, maeo, nnqp, qpmad, zfit, iterative_refinement, odl_pdhg, amaxed_regularization, crystal_ball, directed_divergence, express, ferdor, imaxed, mystic_hybrid, nsduaz, odl_douglas_rachford, qubo, rebunki, rfsp_jul, scipy_direct_method, staysl, tikhonov_legendre, ensemble). Each notebook follows the same template: build a Detector from the built-in RF_GSF response functions, fold an IAEA Compendium reference spectrum into readings, unfold, compute spectral-distance metrics (utils.comparison.compare_spectra), plot and sweep across several IAEA spectra. Documented in docs/examples.rst. New API test module tests/test_all_unfold_methods_api.py (~250 tests): every public unfold_* must be importable from bssunfold.core, accept the canonical (readings, initial_spectrum, ...) signature plus the common Monte-Carlo/housekeeping kwargs, return the standardized result keys on a smoke reading set, accept every documented kwarg non-default and reject unknown kwargs, run end-to-end on the IAEA benchmark, plus headless nbconvert execution of the new notebooks (auto-skipped when an optional backend or jupyter nbconvert is unavailable). Test-infrastructure fixes: tests/test_smt.py defers import z3 behind pytest.importorskip("z3") (collection no longer breaks when the optional z3-solver package is not installed) and the qubo / zfit pytest markers are registered in pyproject.toml. Changed Readme.md - shorten, additional info is in docs.

Publication details

Journal / source
Zenodo (CERN European Organization for Nuclear Research)
Publication date
23 September 2026