Сравнение с оригинальным кодом RECONST#

bssunfold.Detector.unfold_reconst() — это оптимизированная NumPy-реализация алгоритма статистической регуляризации Торчина (STREG1), изначально реализованного на Фортране в программе RECONST.FOR (см. addon/RECONST.FOR).

Обе версии решают систему линейных уравнений:

\[(B \cdot \beta + \Omega \cdot \alpha) \cdot f = A_{vec} \cdot \beta\]

где \(B = A^T \cdot \text{diag}(1/S^2) \cdot A\) — матрица Фишера, \(\Omega\) — матрица сглаживания (5-диагональная), \(A_{vec} = A^T \cdot (F / S^2)\) — взвешенный вектор откликов.

Алгоритм поддерживает 4 режима выбора параметров: α > 0 и β > 0 (фиксированные), α < 0 (автоподбор ω-критерием), β = 0 (автоподбор принципом невязки) и комбинация обоих автоподборов.

Построение матриц B и A_vec#

Fortran (RECONST.FOR:534–544):

 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):

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 (сглаживание)#

Обе реализации используют одинаковую формулу пятидиагональной матрицы сглаживания \(\Omega_{ij}\) с параметром PP:

\[\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)\]

где \(AA_i = 1/(x_i - x_{i-1})\), \(CC_i = 1/(x_{i-1} - x_{i-2})\), \(BB_i = -(AA_i + CC_i)\).

Fortran: хранит в ленточном формате OMO(5, N) и обращается с ручной граничной обработкой (строки 558–561).

Python: строит через _build_omo_matrix() в том же ленточном формате, а при сборке системы использует _omo_to_full() для преобразования в полную симметричную матрицу. Это упрощает код и даёт возможность использовать BLAS-операции с полными матрицами.

Решение системы D·f = A·β#

Fortran — метод Банахевича (RECONST.FOR:702–728):

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):

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 итераций Якоби:

 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

Каждая итерация прибавляет поправку \(\Delta f_i = (A_i - \sum_j D_{ij} f_j) / D_{ii}\) и обрезает отрицательные значения. Четыре прохода с проверкой сходимости через \(\Delta f_i < AINF(4) = 0.01\). Затем проверка, стабилизировалось ли решение — сравниваются абсолютные приращения между последовательными итерациями (строки 621–633). Если сходимость не достигнута за 4 итерации, цикл повторяется до стабилизации.

Python (unfold_reconst.py:311):

FI = np.maximum(FI, 0)

Установлено, что Якоби-итерация численно нестабильна: матрица D не является диагонально-доминантной (число обусловленности может достигать ~10⁷), и поправки \(\Delta f_i\) растут от 10⁻¹¹ до 10² за одну итерацию. Единственный эффект 4-х проходов — обрезание отрицательных значений, которое эквивалентно однократному np.maximum(FI, 0).

Автоматический подбор α и β#

Логика автоматического выбора параметров полностью сохранена:

Соответствие подпрограмм#

Функция

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 итераций.

Производительность#

Benchmark (60×11 система, Intel i7-12700)#

Режим

Fortran (gfortran -O3)

Python (оптимизированная)

Ускорение

Фиксированные α, β

~37 мс

~0.39 мс

×94

Автоподбор α и β

~33 с

~240 мс

×138

Основные факторы ускорения:

  1. Векторизация BLAS — построение B и A_vec через матричные умножения вместо тройных циклов.

  2. Замена Banachiewicz → LAPACK — in-place метод Банахевича уступает dgetrf``+``dgetri на 60×60 матрицах.

  3. Устранение Якоби-постобработки — 4 (иногда больше) итераций с O(N²) заменены на O(N) np.maximum.

  4. Снижение лимитов итераций автоподбора — с 100→50 (внешний) и 200→100 (внутренний бинарный поиск).

Численная устойчивость#

Ограничения Python-версии#

Python-версия не включает редко используемые возможности оригинала:

  • Аддитивный детектор (ADD.COND. из PARAM) — дополнительное уравнение с единичным откликом для энергий выше Ebound. Реализуется вручную через стандартное API детектора при необходимости.

  • Вывод файлов ``###.spe``, ``###.cnt``, ``###.par`` — внутренние форматы RECONST для внешней графики. В Python-версии данные доступны через возвращаемый словарь.

  • Масштабирование счётчиков (FAC, Fnorm) — приводится к единичному масштабу внутри 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 Якоби-постобработки.