Examples ======== Basic Usage ----------- .. code-block:: python import pandas as pd from bssunfold import Detector # Load response functions rf_df = pd.read_csv('response_functions.csv') # Initialize detector detector = Detector(rf_df) # Define readings readings = { 'sphere_1': 150.2, 'sphere_2': 120.5, 'sphere_3': 95.7 } # Unfold with cvxpy result_cvxpy = detector.unfold_cvxpy( readings, regularization=0.001, calculate_errors=True ) # Unfold with Landweber result_landweber = detector.unfold_landweber( readings, max_iterations=500, tolerance=1e-5, calculate_errors=True ) # Unfold with MLEM result_mlem = detector.unfold_mlem_odl( readings, max_iterations=500, calculate_errors=True ) # Unfold with GRAVEL (no extra deps) result_gravel = detector.unfold_gravel( readings, max_iterations=200, tolerance=1e-6 ) # Unfold with MAXED (maximum entropy, no extra deps) result_maxed = detector.unfold_maxed( readings, sigma_factor=0.1 ) # Unfold with Bayes (D'Agostini, no extra deps) result_bayes = detector.unfold_bayes( readings, max_iterations=200, tolerance=1e-3 ) # Unfold with TSVD (truncated SVD, no extra deps) result_tsvd = detector.unfold_tsvd( readings, k=5, method='l_curve' ) # Unfold with StatReg (Turchin, no extra deps) result_statreg = detector.unfold_statreg( readings, unfoldermethod='EmpiricalBayes' ) Parametric Unfolding -------------------- The parametric methods model the neutron spectrum as a weighted sum of thermal, epithermal, and fast components (FRUIT-style model). They are especially useful when the spectrum shape can be approximated by these three components. .. code-block:: python import pandas as pd from bssunfold import Detector detector = Detector(pd.read_csv('response_functions.csv')) readings = {"0in": 0.0003, "2in": 0.0099, "3in": 0.0536, "5in": 0.1841} # Simple lmfit-based parametric unfolding result_param = detector.unfold_parametric( readings, parametric_method='thermal+epithermal+fast', optimizer='lmfit', calculate_errors=True, ) # SQP via cvxpy (requires initial_params or auto-scan) result_cvxpy = detector.unfold_parametric( readings, parametric_method='thermal+epithermal+fast', optimizer='cvxpy', solver_backend='cvxpy:ECOS', calculate_errors=True, ) # Combined: lmfit first, then QP refinement result_combined = detector.unfold_parametric( readings, parametric_method='thermal+epithermal+fast', optimizer='combined', solver_backend='cvxpy', calculate_errors=True, ) # Standalone parametric solver with custom initial guess from bssunfold.core.unfold_parametric import solve_parametric_cvxpy result = solve_parametric_cvxpy( A_matrix=detector.response_functions.values, b_readings=np.array([readings[k] for k in detector.sphere_names]), E=detector.energy_bins, parametric_method='thermal+epithermal+fast', initial_params={'P_th': 1e5, 'P_epi': 1e5, 'P_f': 1e5, 'b': 0.5, 'beta_prime': 0.5, 'alpha': 2.0, 'beta': 0.5}, max_iter=20, tolerance=1e-6, solver_backend='auto', ) BON95 Parametric Unfolding -------------------------- The BON95 method models the lethargy spectrum E*Phi(E) as a linear combination of four components (thermal, epithermal, intermediate, fast) with shape parameters found by grid search and linear coefficients solved by weighted NLS. After parametric fitting, the result is refined by directed-divergence (I-divergence) iterations. .. code-block:: python import pandas as pd from bssunfold import Detector detector = Detector(pd.read_csv('response_functions.csv')) readings = {"0in": 0.0003, "2in": 0.0099, "3in": 0.0536, "5in": 0.1841} # BON95 parametric unfolding (grid search + directed-divergence) result_bon95 = detector.unfold_parametric2( readings, b_range=(0.5, 2.0, 5), # epithermal exponent grid Tf_range=(0.5, 10.0, 5), # fast peak energy grid (MeV) c_range=(0.5, 3.0, 4), # fast peak width grid noise_level=0.05, # 5% measurement uncertainty calculate_errors=True, ) # Compare with FRUIT parametric result_fruit = detector.unfold_parametric(readings, optimizer='lmfit') print("BON95 spectrum shape:", result_bon95['spectrum'].shape) print("FRUIT spectrum shape:", result_fruit['spectrum'].shape) # Standalone BON95 solver from bssunfold.core.unfold_parametric2 import solve_parametric2 import numpy as np E = detector.E_MeV ln_steps = np.zeros(len(E)) log_e = np.log10(E + 1e-15) ln_steps[0] = log_e[1] - log_e[0] ln_steps[-1] = log_e[-1] - log_e[-2] ln_steps[1:-1] = (log_e[2:] - log_e[:-2]) / 2.0 ln_steps *= np.log(10) A = np.array([detector.sensitivities[n] for n in readings]) b = np.array([readings[n] for n in readings]) spectrum, success, msg, nfev = solve_parametric2(A, b, E, ln_steps) print(f"Converged: {success}, message: {msg}") SQP Optimizers ^^^^^^^^^^^^^^ ``unfold_parametric2`` supports multiple optimizers via the ``optimizer`` parameter. The default is ``"grid"`` (exhaustive grid search + NLS). Three additional SQP-based solvers are available: - ``"cvxpy"`` — sequential quadratic programming via cvxpy - ``"qpsolvers"`` — sequential quadratic programming via qpsolvers - ``"combined"`` — grid search followed by SQP refinement .. code-block:: python # Grid search (default) — thorough but slow for fine grids result = detector.unfold_parametric2( readings, optimizer="grid", b_range=(0.5, 2.0, 5), Tf_range=(0.5, 10.0, 5), c_range=(0.5, 3.0, 4), ) # CVXPY SQP — fast, no grid needed result = detector.unfold_parametric2( readings, optimizer="cvxpy", initial_guess=(1.0, 2.0, 1.5), # (b, Tf, c) initial guess ) # QPSolvers SQP — alternative backend (requires OSQP, SCS, or similar) result = detector.unfold_parametric2( readings, optimizer="qpsolvers", solver_backend="osqp", initial_guess=(1.0, 2.0, 1.5), ) # Combined — grid search for coarse optimum, then SQP refinement result = detector.unfold_parametric2( readings, optimizer="combined", b_range=(0.5, 2.0, 5), Tf_range=(0.5, 10.0, 5), c_range=(0.5, 3.0, 4), ) Compressive Sensing (CS) Unfolding ---------------------------------- The ``unfold_cs`` method unfolds a neutron spectrum using **Compressive Sensing (CS)**. The spectrum ``x`` is represented sparsely in a learned dictionary ``D`` as ``x = D @ alpha``, where ``alpha`` is a sparse coefficient vector. The measurement equation ``b = A @ x`` becomes ``b = (A @ D) @ alpha``, which is solved for the sparse ``alpha`` using the **SL0** algorithm. The dictionary is learned with **K-SVD** and sparse coding is performed with **OMP**. This approach is well suited for the highly underdetermined problem where the number of energy groups greatly exceeds the number of detector readings. .. code-block:: python import pandas as pd from bssunfold import Detector detector = Detector(pd.read_csv('response_functions.csv')) readings = {"0in": 0.0003, "2in": 0.0099, "3in": 0.0536, "5in": 0.1841} # Compressive sensing unfolding (K-SVD dictionary + OMP + SL0) result_cs = detector.unfold_cs( readings, n_atoms=80, # number of dictionary atoms sparsity=6, # target sparsity of the coefficient vector max_iterations=200, # SL0 outer iterations random_state=0, # reproducibility calculate_errors=True, ) print("CS spectrum shape:", result_cs['spectrum'].shape) print("CS method:", result_cs['method']) # Standalone CS solver from bssunfold.core.unfold_cs import solve_cs import numpy as np A = np.array([detector.sensitivities[n] for n in readings]) b = np.array([readings[n] for n in readings]) spectrum, iterations, converged = solve_cs(A, b, n_atoms=80, sparsity=6) print(f"Converged: {converged}, iterations: {iterations}") Non-negative K-SVD (NN-KSVD) Unfolding --------------------------------------- The ``unfold_nnksvd`` method implements the BNCT epithermal neutron spectrum unfolding approach of Xu et al. (NIMA 2026). It combines non-negative K-SVD dictionary learning with three sparse-coding strategies: ``nnls_topk`` (proposed), ``omp``, and ``nn_omp``. .. code-block:: python import pandas as pd from bssunfold import Detector detector = Detector(pd.read_csv('response_functions.csv')) readings = {"0in": 0.0003, "2in": 0.0099, "3in": 0.0536, "5in": 0.1841} # NN-KSVD unfolding with the proposed NNLS-TopK coder result = detector.unfold_nnksvd( readings, n_atoms=15, # number of dictionary atoms sparsity=4, # sparsity level K n_dictionary_iterations=80, lambda_tik=0.01, # Tikhonov regularization prior_wt=0.5, # training-sample prior weight sparse_coder='nnls_topk', random_state=42, calculate_errors=True, ) print("NN-KSVD spectrum shape:", result['spectrum'].shape) # Standalone solver from bssunfold.core.unfold_nnksvd import solve_nnksvd_unfold import numpy as np A = np.array([detector.sensitivities[n] for n in readings]) b = np.array([readings[n] for n in readings]) spectrum = solve_nnksvd_unfold(A, b, detector.E_MeV, n_atoms=15, sparsity=4) print("Spectrum shape:", spectrum.shape) Maximum Neutron Energy Cutoff ----------------------------- All ``unfold_*`` methods accept a ``max_neutron_energy`` parameter (in MeV) that forces the reconstructed fluence to zero above the specified energy. This is useful when you know a priori that the neutron field does not contain neutrons above a certain energy, or when you want to isolate a particular energy region. Two internal strategies are used depending on the solver: - **UB array** (QP solvers — ``cvxpy``, ``qpsolvers``, ``docplex``, ``scip``, ``mystic``): the full response matrix is passed to the solver; a per-bin upper bound vector is set to 0 for bins above the cutoff. - **Trimming** (iterative / matrix solvers): the response matrix is sliced to active energy bins, the reduced system is solved, and the result is expanded back to the full grid with zeros above the cutoff. .. code-block:: python import pandas as pd from bssunfold import Detector, RF_LANL df = pd.DataFrame.from_dict(RF_LANL, orient='columns') det = Detector(df) reference = pd.read_csv('MonteCarlo_Calculated_spectra_from_IAEA_Comp_for_comparison.csv') readings = det.get_effective_readings_for_spectra(reference[['E_MeV', 'ISO_ref_Cf252']]) # Restrict unfolding to energies below 10 MeV result = det.unfold_cvxpy(readings, max_neutron_energy=10.0) # Iterative method — same parameter result2 = det.unfold_landweber(readings, max_neutron_energy=10.0) # Verify zero fluence above cutoff above = result['spectrum'][det.E_MeV > 10.0] print(f"Max fluence above 10 MeV: {above.max():.2e}") # 0.00e+00 IAEA Compendium Benchmark Notebooks ------------------------------------ Every public ``unfold_*`` method has a dedicated example notebook that exercises it on the `IAEA Compendium of Monte-Carlo calculated neutron spectra `_. The benchmark CSV ships with the repository at ``tests/MonteCarlo_Calculated_spectra_from_IAEA_Comp_for_comparison.csv`` and contains 20 reference spectra (ISO reference fields, BNCT-like BSA spectra, accelerator-based neutron sources, etc.). Each notebook follows the same template: 1. Build a :class:`~bssunfold.Detector` from the built-in ``RF_GSF`` response functions. 2. Fold a chosen reference spectrum into detector readings via :meth:`~bssunfold.Detector.get_effective_readings_for_spectra`. 3. Run the unfolder with method-specific kwargs. 4. Compute spectral-distance metrics with :func:`~bssunfold.utils.comparison.compare_spectra`. 5. Plot the unfolded spectrum against the ground truth and the per-sphere residuals. 6. Sweep across several IAEA spectra and report the ``comprehensive_score`` for each. Notebooks for the 23 methods that previously lacked a dedicated example (files ``56-…`` through ``78-…`` in ``examples/``), the LOUHI78 (``79-…``) and OSEM-ANLM (``80-…``) benchmark notebooks, and the CUQIpy quality-analysis notebook (``81-…``): ============================================================== ============================ Notebook Method ============================================================== ============================ ``56-gnowee-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_gnowee` ``57-maeo-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_maeo` ``58-nnqp-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_nnqp` ``59-qpmad-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_qpmad` ``60-zfit-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_zfit` ``61-iterative-refinement-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_iterative_refinement` ``62-odl-pdhg-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_odl_pdhg` ``63-amaxed-regularization-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_amaxed_regularization` ``64-crystal-ball-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_crystal_ball` ``65-directed-divergence-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_directed_divergence` ``66-express-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_express` ``67-ferdor-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_ferdor` ``68-imaxed-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_imaxed` ``69-mystic-hybrid-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_mystic_hybrid` ``70-nsduaz-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_nsduaz` ``71-odl-douglas-rachford-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_odl_douglas_rachford` ``72-qubo-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_qubo` ``73-rebunki-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_rebunki` ``74-rfsp-jul-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_rfsp_jul` ``75-scipy-direct-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_scipy_direct_method` ``76-staysl-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_staysl` ``77-tikhonov-legendre-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_tikhonov_legendre` ``78-ensemble-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_ensemble` ``79-louhi-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_louhi` ``80-osem-anlm-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_osem_anlm` ``81-cuqi-quality-analysis-iaea.ipynb`` :meth:`~bssunfold.Detector.unfold_cuqi` ============================================================== ============================ Notebook ``81`` additionally demonstrates the CUQIpy quality-analysis workflow: hierarchical Hybrid Gibbs sampling with Gamma hyperprior, ESS / Gelman-Rubin convergence diagnostics, 95% HPD credible intervals, residual analysis, EURADOS-style energy-group fluences and a comparison against CVXPY / MLEM / Landweber / GRAVEL baselines. The full API surface (imports, signatures, parameter assignment and end-to-end execution) is covered by the test module ``tests/test_all_unfold_methods_api.py``.