pub mod error;
pub mod types;

use crate::error::FilterError;
use crate::types::BoundaryMode;
use num_traits::{Float, FloatConst, NumAssign};

// --- KERNEL GENERATION ---

pub fn gaussian_kernel<T>(sigma: T, truncate: T) -> Vec<T>
where
    T: Float + FloatConst + NumAssign,
{
    // SciPy Logic: int(truncate * sigma + 0.5)
    // We add 0.5 and floor to simulate "Round Half Up" behavior perfectly
    let radius_f = (truncate * sigma + T::from(0.5).unwrap()).floor();
    let radius = radius_f.to_isize().unwrap();
    let size = (2 * radius + 1) as usize;

    if sigma <= T::zero() {
        let mut k = vec![T::zero(); size];
        k[radius as usize] = T::one();
        return k;
    }

    let mut kernel = Vec::with_capacity(size);
    let two = T::from(2.0).unwrap();
    let denom = two * sigma * sigma;

    for i in -radius..=radius {
        let x = T::from(i).unwrap();
        kernel.push((-x * x / denom).exp());
    }

    // Normalize
    let sum: T = kernel.iter().fold(T::zero(), |acc, &x| acc + x);
    for v in &mut kernel {
        *v /= sum;
    }
    kernel
}

// --- PADDING LOGIC ---

fn pad_data<T>(data: &[T], radius: usize, mode: BoundaryMode<T>) -> Vec<T>
where
    T: Float + Copy,
{
    let len = data.len();
    let new_len = len + 2 * radius;
    let mut padded = Vec::with_capacity(new_len);

    // 1. Left Padding
    match mode {
        BoundaryMode::Constant(val) => {
            padded.extend(std::iter::repeat(val).take(radius));
        }
        BoundaryMode::Nearest => {
            let val = data[0];
            padded.extend(std::iter::repeat(val).take(radius));
        }
        BoundaryMode::Reflect => {
            // SciPy-style reflect: d c b a | a b c ...
            for i in (1..=radius).rev() {
                 let src = reflect_index_helper(-(i as isize), len);
                 padded.push(data[src]);
            }
        }
    }

    // 2. Middle (Copy Data)
    padded.extend_from_slice(data);

    // 3. Right Padding
    match mode {
        BoundaryMode::Constant(val) => {
            padded.extend(std::iter::repeat(val).take(radius));
        }
        BoundaryMode::Nearest => {
            let val = data[len - 1];
            padded.extend(std::iter::repeat(val).take(radius));
        }
        BoundaryMode::Reflect => {
             for i in 0..radius {
                 let src = reflect_index_helper((len + i) as isize, len);
                 padded.push(data[src]);
             }
        }
    }

    padded
}

#[inline(always)]
fn reflect_index_helper(idx: isize, len: usize) -> usize {
    if len == 0 { return 0; }
    let n_i = len as isize;
    let mut i = idx;
    if i < 0 { i = -i - 1; }
    let period = 2 * n_i;
    i = ((i % period) + period) % period;
    if i >= n_i { (period - 1 - i) as usize } else { i as usize }
}

// --- CONVOLUTION HELPER (HOT LOOP) ---

/// Inner loop broken out for LLVM optimization.
/// Marked inline(always) to ensure it dissolves into the caller for SIMD.
#[inline(always)]
fn convolve_sample<T>(window: &[T], kernel: &[T]) -> T
where 
    T: Float + NumAssign
{
    let mut acc = T::zero();
    let len = kernel.len(); 
    
    // SAFETY: The caller ensures window and kernel are the same length.
    // This removes bounds checks to allow aggressive SIMD auto-vectorization.
    unsafe {
        for j in 0..len {
            acc += *window.get_unchecked(j) * *kernel.get_unchecked(j);
        }
    }
    acc
}

// --- MAIN ENTRY POINT ---

pub fn gaussian_filter1d<T>(
    data: &[T],
    sigma: T,
    truncate: T,
    mode: BoundaryMode<T>,
) -> Result<Vec<T>, FilterError>
where
    T: Float + FloatConst + NumAssign + std::fmt::Debug,
{
    if data.is_empty() {
        return Err(FilterError::EmptyInput);
    }
    if sigma < T::zero() {
        return Err(FilterError::InvalidSigma(sigma.to_f64().unwrap()));
    }
    // Optimization: If sigma is effectively zero, return copy
    if sigma.abs() < T::epsilon() {
        return Ok(data.to_vec());
    }

    // 1. Prepare Kernel
    let kernel = gaussian_kernel(sigma, truncate);
    let radius = (kernel.len() - 1) / 2;

    // 2. Pad Data
    // Pre-calculating padding allows the inner loop to be branchless
    let padded = pad_data(data, radius, mode);
    
    let mut out = Vec::with_capacity(data.len());
    let k_len = kernel.len();

    // 3. Convolution
    // Iterate only over the valid data region
    for i in 0..data.len() {
        let window = &padded[i..i + k_len];
        out.push(convolve_sample(window, &kernel));
    }

    Ok(out)
}