BSSunfold - Neutron Spectrum Unfolding Package for Bonner Sphere Spectrometers

🔍 Overview

BSSUnfold is a Python package for neutron spectrum unfolding from measurements obtained with Bonner Sphere Spectrometers (BSS). The package implements several mathematical algorithms for solving the inverse problem of unfolding neutron energy spectra from detector readings, with applications in radiation protection, nuclear physics research, and accelerator facilities. Iterative solvers are accelerated with Numba JIT compilation (Landweber, D'Agostini Bayes, Doroshenko, Kaczmarz, MLEM, GRAVEL) for 3–50x speedups.

Contents

• Features
• Installation
• Quick start
• Available Unfolding Methods
• Spectrum Comparison
• Project structure
• Technical requirements
• Authors
• Citing
• Documentation
• Publications

📦 Features

• Multiple Unfolding Algorithms (70+ methods):

  • Tikhonov-type: CVXPY, qpsolvers, Legendre basis, TSVD (truncated SVD), EPIC (Equal Posterior Information Condition)
  • Krylov/hybrid: Lanczos, GKS (Golub-Kahan bidiagonalization + projected GCV/DP/L-curve), CGLS, FISTA (accelerated proximal gradient), Hybrid GMRES
  • Iterative: Landweber, MLEM (pure NumPy + ODL), MLEM-STOP (J-factor stopping), GRAVEL, Doroshenko, Kaczmarz, SART
  • EM family: OSEM (ordered subsets), MAP-EM (penalised one-step-late EM), BSREM (block-sequential regularised EM)
  • Multi-sphere ratio methods: SAND-II (geometric-mean ratios), BUNKI / BUNKI-UT (SPUNIT and BON31G)
  • Classic codes (independent reimplementations): CRYSTAL BALL (direct delta-operator), RFSP-JUL (damped least squares), STAY'SL (single-step Bayesian least squares) — reimplemented from their published algorithmic descriptions; the original codes are proprietary
  • Bayesian: D'Agostini iterative (Bayes), Bayes with spline regularization, zfit likelihood-based inference
  • Maximum Entropy: MAXED (primal log-space dual minimisation), IMAXED, AMAXED, AMAXED-Regularization (Wong 2024 PhD thesis methods)
  • Statistical Regularization: Turchin's method (StatReg, reimplementation of Reconst)
  • Dictionary / Sparse: NN-KSVD (non-negative K-SVD with NNLS-TopK/OMP/NN-OMP sparse coders, Xu et al. NIMA 2026), CS (compressive sensing with K-SVD + SL0)
  • Optimization-based: lmfit (L1/L2/Elastic Net), Scipy direct solvers (CG, GMRES, LSQR), Mystic (direct-search: fmin, Powell, diffev), SMT (exact solving via Z3), Genetic (meta-heuristic: PSO, GA, DE, ES, EP, ABC, GWO, CMA-ES via MEALPY), SCIP (pyscipopt), CPLEX (docplex), QUBO (quantum-inspired annealing)
  • Evolutionary: MAEO (Multi-Algorithm Evolutionary Optimization with NSGA-III, C-TAEA, AGE-MOEA-II, SPEA2)
  • Advanced Proximal: ODL PDHG, ODL Douglas-Rachford (Total Variation regularization)
  • Pipeline: Combined approach for chaining multiple methods
  • Ensemble & Refinement: Ensemble (robust combination of base solvers via weighted average, median, trimmed mean, or best residual), Iterative refinement (two-pass unfold with an auto-selected blending factor), and Bin-wise adaptive (per-bin best method selection from pre-computed benchmark lookup)
  • Parametric: FRUIT-style thermal/epithermal/fast model (lmfit, cvxpy SQP, qpsolvers SQP, combined); BON95 4-component model with directed-divergence iterations; N-spline directed-divergence unfolding (Islamgulov & Lartsev, Atomic Energy 104(5) 2008) with C0/C1 knot continuity and BARS-5/IGRIK/YAGUAR knot presets

• Maximum Energy Cutoff: max_neutron_energy parameter on all unfold_* methods — forces zero fluence above a user-specified energy. QP solvers receive the full response matrix with a ub array; iterative solvers use matrix trimming with automatic result expansion.

• Numba JIT-Accelerated Iterative Solvers:

  • @njit(cache=True) compiled inner loops for Landweber, D'Agostini Bayes, Doroshenko, Kaczmarz, MLEM, GRAVEL
  • 3–50x speedup on iterative solvers (see Performance)
  • Automatic disk caching of compiled code; graceful fallback when numba is not installed

• Radiation Dose Calculations:

  • Effective dose calculations for different irradiation types based on conversion coefficients from 116 publication of International commission on radiological protection (ICRP)

• Comprehensive Data Management:

  • Automatic response function processing
  • Uncertainty quantification via Monte Carlo methods

• Advanced Visualization:

  • Spectrum plotting with uncertainty bands
  • Detector reading comparison

📥 Installation

Using uv (recommended)

uv add bssunfold

Using pip

pip install bssunfold

Using conda

conda install conda-forge::bssunfold

From Source

git clone https://github.com/radiationsafety/bssunfold.git
cd bssunfold
pip install -e .

Optional dependencies

# Basic installation (without additional solvers)
uv add bssunfold

# all methods
uv add bssunfold[all]

# With numba JIT acceleration (recommended for iterative solvers)
uv add "bssunfold[numba]"

# With additional cross-platform solvers (recommended)
uv add "bssunfold[solvers-core]"

# All solvers (Unix/Linux/macOS)
uv add "bssunfold[all-solvers]"

# Windows (all except proxsuite)
uv add "bssunfold[windows]"

# With QP interpretation via pyoptexplain
uv add "bssunfold[interpret]"

# With Bayesian MCMC unfolding (PyMC + ArviZ)
uv add "bssunfold[mcmc]"

Install with all solvers (Unix/Linux/Mac):

uv add bssunfold[all-solvers]

For Windows is recommended to use the following command because of the problem with proxsuite:

uv add bssunfold[windows]

🎯 Quick Start

Open in interactive notebooks:

import pandas as pd
from bssunfold import Detector

# Load response functions
rf_df = pd.read_csv("../data/response_functions/rf_GSF.csv")

# Initialize detector
detector = Detector(rf_df)

# Provide detector readings [reading per second]
readings = {
    "0in": 0.0003,
    "2in": 0.0099,
    "3in": 0.0536,
    "5in": 0.1841,
    "6in": 0.2196,
    "8in": 0.2200,
    "10in": 0.172,
    "12in": 0.120,
    "15in": 0.066,
    "18in": 0.034,
}

# Unfold spectrum using convex optimization
result = detector.unfold_cvxpy(
    readings,
    regularization=1e-4,
    calculate_errors=True
)

# Visualize results
detector.plot_with_uncertainty(result, plot_style='errorbar')

# Calculate and display dose rates
print("Dose rates [pcSv/s]:", result['doserates'])

# Restrict unfolding to energies below 10 MeV (zero fluence above cutoff)
result_limited = detector.unfold_cvxpy(
    readings,
    regularization=1e-4,
    max_neutron_energy=10.0
)

📊 Input Data Structure

Response Functions

Response functions must be provided as a CSV file with the following format:

E_MeV,0in,2in,3in,5in,8in,10in,12in
1.00E-09,0.001,0.005,0.01,0.02,0.03,0.04,0.05
1.00E-08,0.002,0.006,0.012,0.022,0.032,0.042,0.052
...

Detector Readings

Readings should be provided as a dictionary mapping sphere names to measured values:

readings = {
    'sphere_0in': 150.2,   # Bare detector
    'sphere_2in': 120.5,   # 2-inch polyethylene sphere
    'sphere_3in': 95.7,    # 3-inch polyethylene sphere
    # ... additional spheres
}

📦 Built-in Response Functions

The package includes 7 built-in response function datasets for immediate use:

⚠️ Note: RF_EURADOS has a narrower energy range (max 20 MeV) and RF_IHEP has a wider range (max 2000 MeV) compared to the standard 631 MeV used by GSF/PTB/LANL/JINR/Fermilab. Use caution when comparing results across datasets.

from bssunfold import Detector, RF_JINR

# Use built-in response functions directly
detector = Detector(RF_JINR)
result = detector.unfold_cvxpy(readings, regularization=1e-4)

🔢 Dose Conversion Coefficients

The package includes 4 dose conversion coefficient datasets for flexible dose rate calculations:

⚠️ Note: NRB99_2009_effective covers a limited energy range (25 eV – 20 MeV). Values outside this range are set to zero.

from bssunfold import Detector, get_coefficients

# Method 1: Set on Detector (affects all subsequent unfolds)
detector = Detector(cc_type="ICRP74_effective")
result = detector.unfold_cvxpy(readings)

# Method 2: Change after creation
detector.set_dose_coefficients("ICRP74_operational")

# Method 3: Get coefficients directly for custom use
cc = get_coefficients("NRB99_2009_effective")
from bssunfold import interpolate_coefficients
cc_interp = interpolate_coefficients(cc, detector.E_MeV)

⚙️ Available Unfolding Methods

graph TD
    A[Unfolding Methods] --> B[Tikhonov-type]
    A --> J[Krylov/hybrid]
    A --> C[Iterative]
    A --> D[Bayesian]
    A --> E[Maximum Entropy]
    A --> F[Statistical Regularization]
    A --> G[Optimization-based]
    A --> L[Dictionary / Sparse]
    A --> H[Pipeline]
    A --> I[Parametric]
    A --> K[Classic codes]

    B --> B1[unfold_cvxpy]
    B --> B2[unfold_qpsolvers]
    B --> B3[unfold_tsvd]
    B --> B4[unfold_tikhonov_legendre]

    J --> J1[unfold_lanczos]
    J --> J2[unfold_gks]
    J --> J3[unfold_cgls]
    J --> J4[unfold_hybrid_gmres]
    J --> J5[unfold_fista]

    C --> C1[unfold_landweber]
    C --> C2[unfold_mlem]
    C --> C3[unfold_mlem_stop]
    C --> C4[unfold_mlem_odl]
    C --> C5[unfold_gravel]
    C --> C6[unfold_doroshenko]
    C --> C7[unfold_kaczmarz]

    D --> D1[unfold_bayes]
    D --> D2[unfold_bayes_spline_regularization]

    E --> E1[unfold_maxed]
    E --> E2[unfold_imaxed]
    E --> E3[unfold_amaxed]
    E --> E4[unfold_amaxed_regularization]
    F --> F1[unfold_statreg]
    F --> F2[unfold_reconst]

    G --> G1[unfold_lmfit]
    G --> G2[unfold_scipy_direct_method]
    G --> G3[unfold_mystic]
    G --> GH[unfold_mystic_hybrid]
    G --> G4[unfold_smt]
    G --> G5[unfold_genetic]
    G --> G6[unfold_scip]
    G --> G7[unfold_docplex]
    G --> G8[unfold_epic]

    L --> L1[unfold_cs]
    L --> L2[unfold_nnksvd]

    H --> H1[unfold_combined]
    H --> H2[unfold_cascade]
    H --> H3[unfold_composite]
    H --> H4[unfold_ensemble]
    H --> H5[unfold_iterative_refinement]
    H --> H6[unfold_binned]

    I --> I1[unfold_parametric]
    I --> I2[unfold_parametric_cvxpy]
    I --> I3[unfold_parametric_qpsolvers]
    I --> I4[unfold_parametric_combined]
    I --> I5[unfold_parametric2]
    I --> I6[unfold_fruit_like]
    I --> I7[unfold_hybrid_parametric]
    I --> I8[unfold_bayesian_parametric]
    I --> I9[unfold_nspline]

    K --> K1[unfold_crystal_ball]
    K --> K2[unfold_rfsp_jul]
    K --> K3[unfold_staysl]

    style A fill:#4a90d9,color:#fff
    style B fill:#e8f0fe
    style C fill:#e8f0fe
    style D fill:#e8f0fe
    style E fill:#e8f0fe
    style F fill:#e8f0fe
    style G fill:#e8f0fe
    style H fill:#e8f0fe
    style I fill:#e8f0fe
    style J fill:#e8f0fe
    style K fill:#e8f0fe
    style L fill:#e8f0fe

Method Reference Table

Common parameters (shared by most methods): readings, initial_spectrum, calculate_errors, noise_level, n_montecarlo, save_result, random_state.

Pipeline Example

result = detector.unfold_combined(
    readings=readings,
    pipeline=[
        {"method": "cvxpy", "params": {"regularization": 1e-4}},
        {"method": "landweber", "params": {"max_iterations": 2000}},
    ],
    calculate_errors=True,
)

Parametric Example

# FRUIT-style parametric model (thermal + epithermal + fast)
result = detector.unfold_parametric(
    readings=readings,
    parametric_method='thermal+epithermal+fast',
    optimizer='cvxpy',           # or 'lmfit', 'qpsolvers', 'combined'
    solver_backend='cvxpy:ECOS', # or 'qpsolvers:osqp'
    calculate_errors=True,
)

# The parametric model fit yields spectrum components
print(result['doserates'])

N-spline Example

# N-spline unfolding (Islamgulov & Lartsev, Atomic Energy 2008)
# spectrum = exp(a_k + q_k ln E + r_k E) per segment, MIRD directed-divergence loop
result = detector.unfold_nspline(
    readings=readings,
    knots="BARS5_channel",       # or None (auto log-uniform), or explicit MeV knots
    continuity="C0C1",           # spline continuity at interior knots
    relative_uncertainty=0.05,
    max_iterations=300,
)

# Paper's quality control: nev statistic and acceptability bound
print(result["nev"], result["acceptable"])   # nev <= 1 + 2/sqrt(N)
print(result["H_history"])                   # directed-divergence convergence trace

📊 5-Detector Comparison

The dose rate evaluation scripts compare results across 5 detector configurations:

ISO scatter plots with per-detector color coding are generated in tests/iso_plots/ and tests/iaea_compendium_iso_plots/.

📊 Spectrum Comparison

Compare two or more unfolded spectra using a comprehensive set of 41 metrics: 27 simple pairwise diagnostics, plus EURADOS-style integral quantities and spectral diagnostics that are computed automatically when an energy grid is available.

import numpy as np
from bssunfold import Detector

detector = Detector()

r1 = detector.unfold_qpsolvers(readings, save_result=False)
r2 = detector.unfold_cvxpy(readings, save_result=False)

# Compare two results (all simple metrics)
result = detector.compare(r1, r2)
print(result['cosine_similarity'], result['mean_squared_error'])

# Compare with specific metrics
detector.compare(r1, r2, metrics=['cosine_similarity', 'kl_divergence'])

# Compare raw spectra
df = detector.compare(
    np.ones(detector.n_energy_bins),
    np.ones(detector.n_energy_bins) * 2,
    np.ones(detector.n_energy_bins) * 3,
    labels=['Ref', 'A', 'B'],
)
print(df)

# Visual comparison
detector.compare(r1, r2, plot=True, save_to='comparison.png')

# Independent usage
from bssunfold.utils.comparison import compare_spectra, kl_divergence
all_metrics = compare_spectra(s1, s2)
print(kl_divergence(s1, s2))

graph TD
    A[Comparison Metrics<br/>44 total] --> B[Entropy]
    A --> C[Distribution]
    A --> D[Correlation]
    A --> E[Error]
    A --> F[Similarity]
    A --> G[Chi-squared]
    A --> H[Statistical]
    A --> I[EURADOS Integral]
    A --> J[Spectral Diagnostics]
    A --> K[Xu 2026 BNCT]

    B --> B1[kl_divergence]
    B --> B2[cross_entropy]
    B --> B3[entropy_difference_percent]

    C --> C1[wasserstein_dist]
    C --> C2[energy_dist]
    C --> C3[kolmogorov_smirnov_stat]

    D --> D1[pearson_r]
    D --> D2[spearman_r]

    E --> E1[mean_squared_error]
    E --> E2[root_mean_squared_error]
    E --> E3[mean_absolute_error]
    E --> E4[mape]
    E --> E5[r2_score]
    E --> E6[max_error]
    E --> E7[median_absolute_error]

    F --> F1[cosine_similarity]
    F --> F2[mmd_rbf]
    F --> F3[total_flux_ratio]
    F --> F4[spectral_shape_similarity]

    G --> G1[chi_squared]
    G --> G2[g_test]
    G --> G3[freeman_tukey]
    G --> G4[cressie_read]

    H --> H1[anderson_darling]
    H --> H2[wilcoxon_test]
    H --> H3[mannwhitneyu_test]
    H --> H4[standardized_mean_difference]

    I --> I1[fluence_averaged_energy]
    I --> I2[energy_group_fluence]
    I --> I3[dose_averaged_energy]
    I --> I4[ambient_dose_equivalent_rate]

    J --> J1[fluence_difference_percent]
    J --> J2[energy_group_fluence_diff]
    J --> J3[dose_difference_percent]
    J --> J4[fluence_averaged_energy_diff]
    J --> J5[dose_averaged_energy_diff]
    J --> J6[log_lethargy_correlation]
    J --> J7[peak_location_error]
    J --> J8[peak_width_error]
    J --> J9[dose_weighted_error]
    J --> J10[response_matrix_consistency]

    K --> K1[relative_flux_error]
    K --> K2[flux_correlation_coefficient]
    K --> K3[comprehensive_score]

    style A fill:#4a90d9,color:#fff

All 44 Metrics

The 30 simple metrics below are always computed. The EURADOS Integral and Spectral Diagnostics groups require an energy grid (pass energy= or use unfolded result dicts, which carry it) and follow Gomez-Ros et al. 2022. The Xu 2026 (BNCT) group is computed by default alongside the simple metrics and follows Xu et al. (NIMA 2026, https://doi.org/10.1016/j.nima.2026.172070).

The simple metrics are implemented with pure NumPy/SciPy — no extra dependencies required. The integral quantities additionally use the packaged ICRP-116 dose conversion coefficients.

📈 Output Data

The package provides comprehensive output in standardized formats:

Spectrum Results

• Energy grid in MeV
• Unfolded neutron spectrum for the grid of energy bins
• Uncertainty estimates (if calculated)

Dose Calculations

• Effective dose rates for different geometries:
  • AP (Anterior-Posterior)
  • PA (Posterior-Anterior)
  • LLAT (Left Lateral)
  • RLAT (Right Lateral)
  • ROT (Rotational)
  • ISO (Isotropic)

Quality Metrics

• Residual norm
• Iteration counts

📝 Application Areas

Nuclear Research Facilities

• Neutron spectroscopy at particle accelerators
• Reactor neutron field characterization
• Fusion device diagnostics

Radiation Protection

• Workplace monitoring at nuclear power plants
• Medical accelerator facilities
• Industrial radiography installations

Scientific Research

• Space radiation studies
• Cosmic ray neutron measurements
• Nuclear physics experiments

🔬 Advanced Features

Result Management

# List all saved results
results = detector.list_results()
print(f"Available results: {results}")

# Retrieve specific result
result = detector.get_result('20240115_143022_cvxpy')

# Create comprehensive report
report = detector.create_summary_report(
    save_path='unfolding_report.json'
)

# Clear results history
detector.clear_results()

Custom Uncertainty Analysis

# Custom Monte Carlo parameters
result = detector.unfold_cvxpy(
    readings,
    calculate_errors=True,
    n_montecarlo=500,      # Number of samples
    noise_level=0.02       # 2% measurement noise
)

# Access uncertainty data
uncert_mean = result['spectrum_uncert_mean']

📂 Project Structure

bssunfold/
├── CHANGELOG.md
├── CITATION.cff
├── CODE_OF_CONDUCT.md
├── CONTRIBUTING.md
├── LICENSE
├── README.md
├── SECURITY.md
├── TESTS_AND_DOCS.md
├── pyproject.toml
├── uv.lock
├── environment.yml
├── assets/                      # Logos, images
├── conda.recipe/                # Conda build recipe
├── docs/                        # Sphinx documentation
│   ├── index.rst
│   ├── overview.rst             # Methods & metrics overview
│   ├── detector.rst             # Full API reference
│   ├── examples.rst
│   ├── conf.py
│   └── requirements.txt
├── examples/                    # Jupyter notebooks (40 notebooks)
├── scripts/                     # Standalone analysis / benchmark scripts
│   ├── rank_methods.py
│   ├── optimize_defaults_and_new_methods.py
│   ├── dose_rate_evaluation.py
│   ├── dose_rate_iaea_compendium.py
│   ├── dose_rate_qpsolvers_cvxpy.py
│   ├── analyze_iaea_compendium.py
│   └── summarize_iaea.py
├── tools/                       # Dev / CI analysis tools
│   ├── run_dynapyt.py           # DynaPyt dynamic analysis
│   └── regenerate_paper_figures.py
├── tests/                       # Test suite (48 files, 2314 tests)
│   ├── test_all.py
│   ├── test_boost_part1.py
│   ├── test_boost_part2.py
│   ├── test_boost_part3.py
│   ├── test_boost_part4.py
│   ├── test_cascade.py
│   ├── test_classic_unfolders.py
│   ├── test_comparison.py
│   ├── test_composite.py
│   ├── test_coverage.py         # Edge-case & fallback tests
│   ├── test_coverage_boost.py
│   ├── test_cs.py
│   ├── test_detector.py
│   ├── test_docplex.py
│   ├── test_dose_coefficients.py
│   ├── test_em_methods.py
│   ├── test_epic.py
│   ├── test_ferdor.py
│   ├── test_genetic.py
│   ├── test_genetic_improvements.py
│   ├── test_iaea_validation.py
│   ├── test_improvements.py
│   ├── test_interpret.py
│   ├── test_krylov_tv.py
│   ├── test_lanczos.py
│   ├── test_maeo.py
│   ├── test_mcmc.py
│   ├── test_methods2.py
│   ├── test_mlem.py
│   ├── test_mlem_stop.py
│   ├── test_mystic.py
│   ├── test_new_ensemble_refinement.py
│   ├── test_new_methods.py
│   ├── test_new_methods_fixed.py
│   ├── test_new_metrics.py
│   ├── test_new_unfold_methods.py
│   ├── test_nsduaz.py
│   ├── test_randomized_kaczmarz_eki.py
│   ├── test_readings.py
│   ├── test_rebunki.py
│   ├── test_reconst.py
│   ├── test_refactored_fixed.py
│   ├── test_regularization_new_criteria.py
│   ├── test_response_functions.py
│   ├── test_scip.py
│   ├── test_security.py
│   ├── test_smt.py
│   ├── test_unfold_parametric.py
│   ├── test_unfold_parametric2.py
│   └── test_wong2024_methods.py
└── src/
    └── bssunfold/
        ├── __init__.py              # Public API: Detector
        ├── constants.py             # ICRP-116 dose coefficients
        ├── logging_config.py
        ├── platform_check.py        # Solver availability checks
        └── core/
            ├── __init__.py
            ├── _base_unfolder.py
            ├── _bon95.py            # BON95 parametric family
            ├── _em_priors.py        # EM prior functions
            ├── _fruit.py            # FRUIT parametric model + NLS fit
            ├── _interpret_pyopt.py  # Interpretation QP build/solve
            ├── _interpret_report.py # Interpretation report dataclass
            ├── _matrix_utils.py     # SVD, derivative matrix
            ├── _max_energy.py       # max_neutron_energy helpers (UB, mask)
            ├── _montecarlo.py       # MC uncertainty
            ├── _multires.py         # Coarse-to-fine helpers
            ├── _numba_jit.py        # Numba JIT inner loops
            ├── _parametric_shared.py # Shared parametric constants
            ├── _solver_backends.py  # Shared solver-backend resolution
            ├── detector.py          # Main Detector class (70 unfold_* methods)
            ├── dose_calculation.py
            ├── regularization.py    # L-curve, GCV, DP, cosine, NCP, etc.
            ├── unfold_amaxed.py
            ├── unfold_amaxed_regularization.py
            ├── unfold_bayes.py
            ├── unfold_bayes_spline_regularization.py
            ├── unfold_bayesian_parametric.py
            ├── unfold_bsrem.py
            ├── unfold_bunki.py
            ├── unfold_bunkiut.py
            ├── unfold_cascade.py
            ├── unfold_cgls.py
            ├── unfold_combined.py
            ├── unfold_composite.py
            ├── unfold_crystal_ball.py
            ├── unfold_cs.py
            ├── unfold_cvxpy.py
            ├── unfold_docplex.py
            ├── unfold_doroshenko.py
            ├── unfold_eki.py
            ├── unfold_ensemble.py
            ├── unfold_epic.py
            ├── unfold_ferdor.py
            ├── unfold_fista.py
            ├── unfold_fruit_like.py
            ├── unfold_genetic.py
            ├── unfold_gks.py
            ├── unfold_gravel.py
            ├── unfold_hybrid_gmres.py
            ├── unfold_hybrid_parametric.py
            ├── unfold_imaxed.py
            ├── unfold_interpret.py
            ├── unfold_iterative_refinement.py
            ├── unfold_kaczmarz.py
            ├── unfold_lanczos.py
            ├── unfold_landweber.py
            ├── unfold_lmfit.py
            ├── unfold_maeo.py
            ├── unfold_mapem.py
            ├── unfold_maxed.py
            ├── unfold_mcmc.py
            ├── unfold_mlem.py
            ├── unfold_mlem_odl.py
            ├── unfold_mlem_stop.py
            ├── unfold_mystic.py
            ├── unfold_nnksvd.py
            ├── unfold_nsduaz.py
            ├── unfold_nspline.py
            ├── unfold_odl_advanced.py
            ├── unfold_osem.py
            ├── unfold_parametric.py
            ├── unfold_parametric2.py
            ├── unfold_qpsolvers.py
            ├── unfold_qubo.py
            ├── unfold_randomized_kaczmarz.py
            ├── unfold_rebunki.py
            ├── unfold_reconst.py
            ├── unfold_rfsp_jul.py
            ├── unfold_sandii.py
            ├── unfold_sart.py
            ├── unfold_scip.py
            ├── unfold_scipy_direct_method.py
            ├── unfold_smt.py
            ├── unfold_statreg.py
            ├── unfold_staysl.py
            ├── unfold_tikhonov_legendre.py
            ├── unfold_tikhonov_tv.py
            ├── unfold_tsvd.py
            └── unfold_zfit.py

🔧 Technical Requirements

Core Requirements

• Python 3.11+
• NumPy, SciPy, Pandas, Matplotlib
• cvxpy[ecos] — convex optimisation framework (CVXPY-based methods)

Optional Backends

• numba — JIT compilation for iterative solvers (3–50x speedup)
• pytikhonov — L-curve / GCV / DP regularisation (Tikhonov-type methods)
• qpsolvers[solvers-core] — QP solvers (unfold_qpsolvers)
• mystic — constrained/direct-search optimization (unfold_mystic)
• z3-solver — SMT exact solving (unfold_smt)
• pyscipopt — SCIP Optimization Suite interface (unfold_scip)
• docplex + cplex — IBM CPLEX modeling & engine (unfold_docplex)
• mealpy — population-based meta-heuristic optimization (unfold_genetic)
• lmfit — L1/L2/Elastic Net regularisation (unfold_lmfit)
• odl — Operator Discretization Library (unfold_mlem_odl)
• pymc + arviz — Bayesian MCMC/NUTS sampling (unfold_mcmc)

All other methods (GRAVEL, MAXED, Bayes, StatReg, Reconst, TSVD, ScipyDirect, Landweber, Kaczmarz, Doroshenko, MLEM, TikhonovLegendre, NSpline) have no extra dependencies beyond NumPy/SciPy.

See pyproject.toml for version constraints.

Performance

All iterative solvers use Numba JIT-compiled inner loops when numba is installed, with automatic fallback to pure Python.

Benchmarks on 60-bin energy grid, 500 iterations, macOS arm64.

Install numba for the best performance:

uv add bssunfold[numba]

📖 Citation

If you use BSSUnfold in your research, please cite paper:

@article{chizhov2024neutron,
  title={Neutron spectra unfolding from Bonner spectrometer readings by the regularization method using the Legendre polynomials},
  author={Chizhov, K and Beskrovnaya, L and Chizhov, A},
  journal={Physics of Particles and Nuclei},
  volume={55},
  number={3},
  pages={532--534},
  year={2024},
  publisher={Springer}
}

or software:

@misc{konstantin_radiationsafetybssunfold_2025,
        title = {Radiationsafety/bssunfold},
        copyright = {GNU General Public License v3.0 only},
        shorttitle = {Radiationsafety/bssunfold},
        url = {https://zenodo.org/doi/10.5281/zenodo.18056376},
        abstract = {first published version of package},
        urldate = {2026-01-12},
        publisher = {Zenodo},
        author = {Chizhov, Konstantin},
        month = dec,
        year = {2025},
        doi = {10.5281/ZENODO.18056376},
}

💬 Contributing

Contributions are welcome! Please see CONTRIBUTING.md for guidelines.

1. Fork the repository
2. Create a feature branch
3. Add tests for new functionality
4. Ensure all tests pass
5. Submit a pull request

📘 Documentation

Documentation and API reference is available in /docs folder. Theory and methodology in the research paper, examples of usage in /examples folder. Check the https://bssunfold.readthedocs.io/en/latest/

📄 License

This project is licensed under the GNU General Public License v3.0 - see the LICENSE file for details.

💬 Support

For questions, bug reports, or feature requests:

• Open an issue on GitHub
• Contact: kchizhov@jinr.ru

💻 Authors

• Konstantin Chizhov
• Alexei Chizhov
• Dmitry Borschev
• Maria Akimochkina

🌐 Acknowledgments

• ICRP and IAEA for data
• Contributors and testers
• Joint Institute for Nuclear Research (JINR)
• University "Dubna", School of Big Data Analytics

🎓 Publications

1. Chizhov A. V., Chizhov K. A. TSVD-Based Iterative Algorithm of Landweber for Neutron Spectra Unfolding by Bonner Multi-Sphere Spectrometer Readings // Phys. Part. Nuclei. 2026. Т. 57. № 4. С. 750–752. https://doi.org/10.1134/S1063779626700735
2. Чижов К.А., Чижов А.В., Борщев Д.С., Акимочкина М.А. Методы решения обратных задач для обработки результатов измерений на примере восстановления спектра нейтронов, Тридцать третья международная конференция "Математика. Компьютер. Образование, г. Дубна, 26 – 31 января 2026 г., https://mce.su
3. Chizhov, K., Chizhov, A. Optimization of the Neutron Spectrum Unfolding Algorithm Using Shifted Legendre Polynomials Based on Weighted Tikhonov Regularization. Phys. Part. Nuclei 56, 1395–1399 (2025). https://doi.org/10.1134/S106377962570056X
4. Chizhov K., Beskrovnaya L., Chizhov A. Neutron spectrum unfolding method based on shifted Legendre polynomials, its application to the IREN facility // Phys. Part. Nucl. Lett. — 2025. — V. 22, no. 2. — P. 337–340. — DOI: https://doi.org/10.1134/S154747712470239X
5. Chizhov K., Beskrovnaya L., Chizhov A. Neutron spectra unfolding from Bonner spectrometer readings by the regularization method using the Legendre polynomials // Phys. Part. Nucl. — 2024. — V. 55. — P. 532–534. — DOI: https://doi.org/10.1134/S1063779624030298
6. Chizhov K., Chizhov A. Optimization approach to neutron spectra unfolding with Bonner multi-sphere spectrometer // Math. Model. — 2024. — V. 7. — P. 89–90.
7. Чижов А. В., Чижов К. А. Восстановление спектров опорных нейтронных полей на Фазотроне (ОИЯИ) на основе показаний многошарового спектрометра Боннера методом усеченного сингулярного разложения Тезисы Трудов LXI Всероссийской конференции по физике РУДН 19 - 23 мая 2025.
8. Chizhov, K., Chizhov, A., TSVD-based neutron spectra unfolding by Bonner multi-sphere spectrometer readings with iteration procedure, proceedings of the International Conference "Distributed Computing and Grid-technologies in Science and Education".
9. Белый А.А., Стариковская М.Д., Чижов К.А. Разработка веб-приложения для эксперимента по восстановлению спектра нейтронов с применением алгоритмов нейронный сетей. Системный анализ в науке и образовании. 2025;(2):49–57.
10. Starikovskaya MD, Chizhov KA. Neutron spectrum unfolding based on random forest algorithm and generated training sample. In Российский университет дружбы народов им. П. Лумумбы; 2025 [cited 2025 Dec 25]. p. 389–94. Available from: https://www.elibrary.ru/item.asp?id=83014906
11. Chizhov KA, Bely AA, Starikovskaia MD, Volkov EN. Восстановление энергетического спектра потока нейтронного излучения с помощью алгоритма машинного обучения «случайный лес». Современные информационные технологии и ИТ-образование. 2024 Dec 15 [cited 2025 Apr 9]; 20(4). Available from: http://sitito.cs.msu.ru/index.php/SITITO/article/view/1167

📘 References

1. Compendium of neutron spectra and detector responses for radiation protection purposes: supplement to technical reports series no. 318. — Vienna: International Atomic Energy Agency, 2001. — Technical reports series no. 403. — STI/DOC/010/403. — ISBN 92-0-102201-8.
2. Diamond, S. and Boyd, S., 2016. CVXPY: A Python-embedded modeling language for convex optimization. Journal of Machine Learning Research, 17(83), pp.1-5.

BSSUnfold - Professional neutron spectrum unfolding for radiation science and nuclear applications.