#define EIGEN_MPL2_ONLY // Ограничиваем подключаемые части Eigen лицензией MPL-2.0.
#include <Eigen/QR> // Подключаем QR с выбором ведущего столбца для полиномиального МНК.
#include <cmath> // Проверяем конечность вычисленных коэффициентов перед передачей в C.
#include <new> // Перехватываем ошибку выделения памяти внутри границы C ABI.
#include "set_signal.h" // Используем общую подготовку точек, вычисление кривой и стабильный C ABI.

static int eigen_fit(const double *x,const double *y,size_t count,unsigned degree,double *coefficients) { // Вычисляем только коэффициенты по уже нормированным уникальным точкам.
    try { // Исключения C++ не должны пересекать границу ctypes или C-ядра.
        Eigen::MatrixXd design(count,degree+1); // Матрица имеет одну строку на точку и до шести степеней, включая свободный член.
        for(size_t i=0;i<count;++i) { // Заполняем нормированную матрицу Вандермонда в том же базисе, что исходный МНК.
            double power=1.0; // Нулевая степень задаёт столбец свободного члена без вызова pow.
            for(unsigned j=0;j<=degree;++j) { // Последовательно строим степени от нулевой до выбранной пользователем.
                design(i,j)=power; // Сохраняем коэффициент при соответствующей степени в строке измерения.
                power*=2.0*x[i]-1.0; // Переводим нормированное время из [0,1] в [-1,1], сохраняя обусловленность исходного алгоритма.
            } // Заканчиваем заполнение строки матрицы по одной временной точке.
        } // Все точки участвуют в МНК с одинаковым весом после общего усреднения повторов X.
        Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr(design); // Используем отражения Хаусхолдера с перестановкой столбцов вместо собственной ортогонализации.
        qr.setThreshold(1e-12); // Относительный порог ранга соответствует порядку допуска исходного решателя; он не является абсолютной погрешностью сигнала.
        if(qr.rank()!=static_cast<Eigen::Index>(degree+1)) return 3; // Отклоняем вырожденный базис, не выдавая неоднозначные коэффициенты за корректную кривую.
        Eigen::Map<const Eigen::VectorXd> values(y,count); // Читаем нормированные амплитуды из памяти C без дополнительного копирования входного вектора.
        Eigen::VectorXd solution=qr.solve(values); // Находим минимум суммы квадратов без формирования нормальных уравнений AᵀA.
        for(unsigned j=0;j<=degree;++j) { // Возвращаем степени в исходном порядке: Eigen отменяет перестановку столбцов при solve.
            if(!std::isfinite(solution[j])) return 2; // Не разрешаем NaN или бесконечности попасть в оценку кривой и RMSE.
            coefficients[j]=solution[j]; // Передаём коэффициент в буфер общего C-алгоритма вычисления полинома.
        } // Весь выходной полином проверен до передачи управления C-ядру.
        return 0; // Сообщаем общий для двух реализаций код успешного расчёта.
    } catch(const std::bad_alloc &) { // Преобразуем нехватку памяти необязательного настольного решателя в детерминированный код.
        return 5; // Код 5 обозначает невозможность выполнить расчёт, включая выделение матриц.
    } catch(...) { // Защищаем потребителей C ABI от остальных исключений библиотеки C++.
        return 5; // Возвращаем явную ошибку без подмены выбранного пользователем решателя.
    } // Завершаем защиту границы между C++ и C.
} // Решатель не хранит глобальных данных и допускает независимые параллельные вызовы.

extern "C" PCAN_ABI_API int set_signal_reconstruct_eigen( // Отдельный экспорт сохраняет возможность загрузить модуль без изменения основной DLL.
    const double *x,const double *y,size_t count,int method,unsigned degree, // Формат точек и предел степени совпадают с исходным API.
    size_t output_count,int endpoint,double *out_x,double *out_y, // Временная сетка и включение конечной точки рассчитываются общим C-ядром.
    double *meta,double *work,size_t work_count) { // Метаданные и рабочая память сохраняют формат аннотаций формулы GUI.
    if(method!=SET_SIGNAL_POLYNOMIAL) return 1; // Eigen доступен только для МНК; остальные методы не маскируются под альтернативный решатель.
    return set_signal_reconstruct_with_fit(x,y,count,method,degree,output_count,endpoint,out_x,out_y,meta,work,work_count,eigen_fit); // Повторно используем C99-подготовку, обработку повторов, вычисление кривой и RMSE.
} // Завершаем экспорт с соглашением вызова C для x86 и x64.