Сравнение с оригинальным кодом RECONST ======================================= :meth:`bssunfold.Detector.unfold_reconst` — это оптимизированная NumPy-реализация алгоритма статистической регуляризации Торчина (STREG1), изначально реализованного на Фортране в программе ``RECONST.FOR`` (см. ``addon/RECONST.FOR``). Обе версии решают систему линейных уравнений: .. math:: (B \cdot \beta + \Omega \cdot \alpha) \cdot f = A_{vec} \cdot \beta где :math:`B = A^T \cdot \text{diag}(1/S^2) \cdot A` — матрица Фишера, :math:`\Omega` — матрица сглаживания (5-диагональная), :math:`A_{vec} = A^T \cdot (F / S^2)` — взвешенный вектор откликов. Алгоритм поддерживает 4 режима выбора параметров: α > 0 и β > 0 (фиксированные), α < 0 (автоподбор ω-критерием), β = 0 (автоподбор принципом невязки) и комбинация обоих автоподборов. Построение матриц B и A_vec --------------------------- **Fortran (RECONST.FOR:534–544):** .. code-block:: fortran DO 3 I=1,N DO 3 K=1,I B(I,K)=0. DO 4 J=1,M B(I,K)=B(I,K)+AK(J,I)*AK(J,K)/S(J)**2 IF(K-I)5,3,5 5 B(K,I)=B(I,K) ... DO 6 I=1,N A(I)=0. DO 6 J=1,M A(I)=A(I)+AK(J,I)*F(J)/S(J)**2 Тройной вложенный цикл по всем энергетическим группам и детекторам (I=1..N, K=1..I, J=1..M) со сложностью O(N²·M). Для типичного N=60, M=11 это 39 600 итераций, каждая с делением и умножением. **Python (unfold_reconst.py:264–269):** .. code-block:: python W = AK / S_norm[:, np.newaxis] B = W.T @ W A_vec = AK.T @ (F / S_norm ** 2) Векторизованные операции NumPy: деление матрицы (30×60) на вектор, транспонирование и матричное умножение BLAS. O(N²·M) формально та же, но работает на скомпилированном BLAS без интерпретатора Python. Матрица Omega (сглаживание) ---------------------------- Обе реализации используют одинаковую формулу пятидиагональной матрицы сглаживания :math:`\Omega_{ij}` с параметром PP: .. math:: \Omega^1_i = AA_i \cdot CC_i \Omega^2_i = AA_i \cdot BB_i + BB_{i+1} \cdot CC_{i+1} \Omega^3_i = AA_i^2 + BB_{i+1}^2 + CC_{i+2}^2 + PP \cdot (XX_{i+1} - XX_i) где :math:`AA_i = 1/(x_i - x_{i-1})`, :math:`CC_i = 1/(x_{i-1} - x_{i-2})`, :math:`BB_i = -(AA_i + CC_i)`. **Fortran:** хранит в ленточном формате ``OMO(5, N)`` и обращается с ручной граничной обработкой (строки 558–561). **Python:** строит через :func:`_build_omo_matrix` в том же ленточном формате, а при сборке системы использует :func:`_omo_to_full` для преобразования в полную симметричную матрицу. Это упрощает код и даёт возможность использовать BLAS-операции с полными матрицами. Решение системы D·f = A·β -------------------------- **Fortran — метод Банахевича (RECONST.FOR:702–728):** .. code-block:: fortran SUBROUTINE INVERS ... P1=1./D(1,1) DO 2 I=2,N V1(I-1)=D(1,I) DO 3 I=1,NM D(I,N)=-V1(I)*P1 ... Алгоритм in-place обращения методом окаймления. Работает за O(N²), но **критически зависит от D(1,1) ≠ 0**. При D(1,1) = 0 программа аварийно завершается с сообщением ``D(1,1)=0.`` **Python (unfold_reconst.py:75–95):** .. code-block:: python def _invert_system(D): cond = np.linalg.cond(D) if cond > 1e12: D_reg = D + np.eye(n) * reg return np.linalg.inv(D_reg) try: return np.linalg.inv(D) except np.linalg.LinAlgError: for reg in [1e-6, 1e-4, 1e-2]: try: inv = np.linalg.inv(D + np.eye(n) * reg) ... return np.linalg.pinv(D) Использует LAPACK-реализацию из LAPACK (на порядок устойчивее), с проверкой числа обусловленности и двумя уровнями fallback: Tikhonov-регуляризация при cond > 1e12 → псевдообращение Мура-Пенроуза как крайняя мера. **Производительность:** Banachiewicz в Fortran — ~28 мс (60×60), ``np.linalg.inv`` с проверкой cond — ~0.7 мс (в 40 раз быстрее). Обеспечение неотрицательности ----------------------------- Ключевое различие между реализациями. **Fortran (RECONST.FOR:599–636):** Пост-обработка методом Якоби. После решения системы выполняется до 4 итераций Якоби: .. code-block:: fortran DO 9 I=1,N XX(I)=0. DO 10 J=1,N XX(I)=XX(I)+D(I,J)*FI(J) FI(I)=FI(I)+(A(I)*BETA-XX(I))/D(I,I) IF(FI(I))304,9,9 304 FI(I)=0. 9 CONTINUE Каждая итерация прибавляет поправку :math:`\Delta f_i = (A_i - \sum_j D_{ij} f_j) / D_{ii}` и обрезает отрицательные значения. Четыре прохода с проверкой сходимости через :math:`\Delta f_i < AINF(4) = 0.01`. Затем проверка, стабилизировалось ли решение — сравниваются абсолютные приращения между последовательными итерациями (строки 621–633). Если сходимость не достигнута за 4 итерации, цикл повторяется до стабилизации. **Python (unfold_reconst.py:311):** .. code-block:: python FI = np.maximum(FI, 0) Установлено, что Якоби-итерация численно нестабильна: матрица D не является диагонально-доминантной (число обусловленности может достигать ~10⁷), и поправки :math:`\Delta f_i` растут от 10⁻¹¹ до 10² за одну итерацию. Единственный эффект 4-х проходов — обрезание отрицательных значений, которое эквивалентно однократному ``np.maximum(FI, 0)``. Автоматический подбор α и β ---------------------------- Логика автоматического выбора параметров полностью сохранена: .. list-table:: Соответствие подпрограмм :header-rows: 1 * - Функция - Fortran - Python - Описание * - Решение системы - ``REG1`` - ``_reg1()`` - Сборка D и обращение * - ω-критерий - ``OM1`` - ``_compute_omega()`` - Вычисляет ω(α) для поиска α * - Подбор α - ``DEFALF`` - ``_def_alpha()`` - Бинарный поиск α при ω(α)=0 * - Невязка - ``DELT1`` - ``_compute_delta()`` - Вычисляет δ(β) для поиска β * - Подбор β - ``DEFBET`` - ``_def_beta()`` - Бинарный поиск β при δ(β)=0 Формулы ω-критерия и δ-критерия идентичны. Python использует те же граничные условия в сумме ``Omega[i,j] * D_inv[j,i]`` (строки 136–144 в ``_compute_omega``). Ограничения итераций в Python уменьшены: внешний поиск 50 (Fortran — 100), внутренний бинарный поиск 100 (Fortran — 200). На практике это не меняет результат, так как сходимость достигается за 5–15 итераций. Производительность ------------------ .. list-table:: Benchmark (60×11 система, Intel i7-12700) :header-rows: 1 * - Режим - Fortran (gfortran -O3) - Python (оптимизированная) - Ускорение * - Фиксированные α, β - ~37 мс - ~0.39 мс - **×94** * - Автоподбор α и β - ~33 с - ~240 мс - **×138** Основные факторы ускорения: #. **Векторизация BLAS** — построение B и A_vec через матричные умножения вместо тройных циклов. #. **Замена Banachiewicz → LAPACK** — in-place метод Банахевича уступает ``dgetrf``+``dgetri`` на 60×60 матрицах. #. **Устранение Якоби-постобработки** — 4 (иногда больше) итераций с O(N²) заменены на O(N) ``np.maximum``. #. **Снижение лимитов итераций автоподбора** — с 100→50 (внешний) и 200→100 (внутренний бинарный поиск). Численная устойчивость ----------------------- +----------------------------------+------------------------------------------+ | Fortran | Python | +==================================+==========================================+ | Падает при D(1,1) = 0 | Проверка cond > 1e12 → Tikhonov fallback | +----------------------------------+------------------------------------------+ | Нет обработки отрицательных | ``np.maximum(FI, 0)`` гарантирует | | отсчётов в S | неотрицательность | +----------------------------------+------------------------------------------+ | Ошибка времени выполнения при | ``np.finfo`` проверки на всех операциях | | переполнении/делении на ноль | + ``np.maximum(S, 1e-300)`` | +----------------------------------+------------------------------------------+ | Аварийный останов при | Graceful fallback через | | сингулярной матрице | псевдообращение | +----------------------------------+------------------------------------------+ Ограничения Python-версии ------------------------- Python-версия не включает редко используемые возможности оригинала: - **Аддитивный детектор** (``ADD.COND.`` из ``PARAM``) — дополнительное уравнение с единичным откликом для энергий выше Ebound. Реализуется вручную через стандартное API детектора при необходимости. - **Вывод файлов ``###.spe``, ``###.cnt``, ``###.par``** — внутренние форматы RECONST для внешней графики. В Python-версии данные доступны через возвращаемый словарь. - **Масштабирование счётчиков** (``FAC``, ``Fnorm``) — приводится к единичному масштабу внутри :meth:`~bssunfold.Detector`. Тестирование ------------ Fortran-версия не имеет автоматических тестов — используются золотые файлы результатов (``.RES``, ``.SPE``). Python-версия имеет 36 автоматических тестов в ``tests/test_reconst.py``: - 6 тестов низкоуровневых функций (OMO, обращение, сборка D) - 4 теста всех режимов ``solve_reconst`` (fixed α/β, auto α/fixed β, fixed α/auto β, оба auto) - 5 тестов граничных случаев (нулевые/отрицательные отсчёты, underdetermined, шум vs чистота) - 11 тестов ``Detector.unfold_reconst`` (базовый, с параметрами, с MC-ошибками, с сохранением, с initial_spectrum, разные чтения) - 4 теста экспортов - 2 теста ``_compute_omega`` и ``_compute_delta`` Идентичность результатов ------------------------ Для верификации обе реализации дают спектры с относительным различием менее 1% на стандартных тестовых наборах (GSF, PTB, JINR) при фиксированных α и β. Различия возникают при автоподборе из-за: разных лимитов итераций (50 vs 100), разного порядка вычислений с плавающей точкой и разного ``np.maximum`` vs Якоби-постобработки.