use ndarray::{Array1, Array2};

/// 计算1D高斯卷积核，支持高阶微分(order)
pub fn gaussian_kernel1d(sigma: f64, order: usize, radius: usize) -> Result<Array1<f64>, String> {
    let sigma2 = sigma * sigma;

    // x范围[-radius, radius]
    let x: Array1<f64> = Array1::from_iter((-(radius as isize)..=radius as isize).map(|v| v as f64));
    // phi_x = exp(-0.5 / sigma2 * x^2)
    let mut phi_x = x.mapv(|v| (-0.5 / sigma2 * v * v).exp());
    let sum_phi = phi_x.sum();
    phi_x /= sum_phi;

    if order == 0 {
        return Ok(phi_x);
    }

    // 计算q(x)系数向量 q = [1, 0, 0, ..., 0]
    let mut q = Array1::<f64>::zeros(order + 1);
    q[0] = 1.0;

    // exponent_range = [0, 1, ..., order]
    let exponent_range: Vec<usize> = (0..=order).collect();

    // 构造矩阵 D (order+1 x order+1)，
    // D[i, i+1] = exponent_range[i+1]
    let mut d = Array2::<f64>::zeros((order + 1, order + 1));
    for i in 0..order {
        d[[i, i + 1]] = (exponent_range[i + 1]) as f64;
    }

    // 构造矩阵 P (order+1 x order+1),
    // P[i, i-1] = 1 / -sigma2
    let mut p = Array2::<f64>::zeros((order + 1, order + 1));
    for i in 1..=order {
        p[[i, i - 1]] = -1.0 / sigma2;
    }

    // Q_deriv = D + P
    let q_deriv = d + p;

    // q = (Q_deriv)^order * q 初始q=[1,0,...]
    for _ in 0..order {
        q = q_deriv.dot(&q);
    }

    // 计算 (x^exponent_range) dot q
    // 先构造 x的幂次矩阵 shape=(x.len(), order+1)
    let x_len = x.len();
    let mut x_matrix = Array2::<f64>::zeros((x_len, order + 1));
    for i in 0..x_len {
        for j in 0..=order {
            x_matrix[[i, j]] = x[i].powi(j as i32);
        }
    }

    // x_matrix dot q
    let q_poly = x_matrix.dot(&q);

    // 返回 q_poly * phi_x
    Ok(&q_poly * &phi_x)
}