use crate::numerical::elementary::eval_expr;
use crate::symbolic::core::Expr;
use std::collections::HashMap;
use std::convert::TryInto;
use std::sync::Arc;
/// Computes the numerical Taylor series expansion of a function around a point.
///
/// This function calculates the coefficients of the Taylor series by numerically
/// evaluating the function and its derivatives at the `at_point`. It then returns
/// a closure that can evaluate the resulting Taylor polynomial at any `x`.
///
/// # Arguments
/// * `f` - The symbolic expression for the function `f(x)`.
/// * `var` - The variable `x`.
/// * `at_point` - The point `a` around which to expand the series.
/// * `order` - The maximum order `N` of the series to compute.
///
/// # Returns
/// A `Result` containing a closure `Arc<dyn Fn(f64) -> f64>` that evaluates the Taylor polynomial,
/// or an error string if evaluation fails.
pub fn taylor_series_numerical(
    f: &Expr,
    var: &str,
    at_point: f64,
    order: usize,
) -> Result<Arc<dyn Fn(f64) -> f64>, String> {
    let mut coeffs = Vec::with_capacity(order + 1);
    let mut current_f = f.clone();
    let mut factorial = 1.0;
    let mut vars_map = HashMap::new();
    vars_map.insert(var.to_string(), at_point);
    coeffs.push(eval_expr(&current_f, &vars_map)?);
    for i in 1..=order {
        current_f = crate::symbolic::calculus::differentiate(&current_f, var);
        factorial *= i as f64;
        let coeff_val = eval_expr(&current_f, &vars_map)? / factorial;
        coeffs.push(coeff_val);
    }
    let a = at_point;
    let taylor_poly = move |x: f64| -> f64 {
        let mut sum = 0.0;
        for (i, coeff) in coeffs.iter().enumerate() {
            let power: i32 = i
                .try_into()
                .expect("Series index 'i' is too large to fit in i32");
            sum += coeff * (x - a).powi(power);
        }
        sum
    };
    Ok(Arc::new(taylor_poly))
}
