#ifndef ADAPTIVESIMPSONINTEGRATE_H
#define ADAPTIVESIMPSONINTEGRATE_H

namespace AdaptiveSimpsonIntegrate {

#include <armadillo>
#include <cmath>
#include <iostream>
//#include <iomanip>

// 定义被积函数：y = A*exp(-(x-P)²/(2*delt²)) + H/(1+exp((x-P)/W)) + C
struct FitFunction {
    // 拟合公式的参数
    double A;    // 高斯项振幅
    double P;    // 中心位置
    double delt; // 高斯项宽度
    double H;    // Sigmoid项高度
    double W;    // Sigmoid项宽度
    double C;    // 常数项

    // 重载()运算符，计算给定点x的函数值
    double operator()(double x) const {
        // 高斯项：A*exp(-(x-P)²/(2*delt²))
        double gaussian_term = A * ::std::exp(-::std::pow(x - P, 2) / (2 * ::std::pow(delt, 2)));

        // Sigmoid项：H/(1+exp((x-P)/W))
        double sigmoid_term = H / (1 + ::std::exp((x - P) / W));

        // 常数项 + 两项之和
        return gaussian_term + sigmoid_term + C;
    }
};

// 辛普森积分核心计算（单区间）
template<typename Func>
double simpson(double a, double b, const Func& f) {
    double h = (b - a) / 2.0;
    double fa = f(a), fb = f(b), fc = f(a + h);
    return h * (fa + 4 * fc + fb) / 3.0;
}

// 自适应辛普森积分（递归实现，控制精度）
template<typename Func>
double adaptive_simpson(double a, double b, double eps, const Func& f,
                        double whole, double fa, double fb, double fc, int depth) {
    // 防止递归过深（避免栈溢出）
    if (depth > 20) return whole;

    double c = (a + b) / 2.0;
    double h = (b - a) / 4.0;
    double fd = f(a + h), fe = f(b - h);

    // 计算左右子区间的辛普森值
    double left = simpson(a, c, f);
    double right = simpson(c, b, f);

    // 精度判断：满足则返回，不满足则递归细分
    if (::std::abs(left + right - whole) <= 15 * eps) {
        return left + right + (left + right - whole) / 15.0;
    }

    return adaptive_simpson(a, c, eps/2, f, left, fa, fc, fd, depth+1) +
           adaptive_simpson(c, b, eps/2, f, right, fc, fb, fe, depth+1);
}

// 积分入口函数（对外接口）
template<typename Func>
double integrate(Func f, double a, double b, double eps = 1e-8) {
    double c = (a + b) / 2.0;
    double fa = f(a), fb = f(b), fc = f(c);
    double whole = simpson(a, b, f);
    return adaptive_simpson(a, b, eps, f, whole, fa, fb, fc, 0);
}

}

#endif // ADAPTIVESIMPSONINTEGRATE_H