scirs2-integrate 0.4.2

//! Quasi-Monte Carlo integration
//!
//! This module provides Quasi-Monte Carlo integration methods for numerical integration
//! of multidimensional functions.
//!
//! QMC integration uses low-discrepancy sequences to generate evaluation points,
//! offering better convergence rates than traditional Monte Carlo methods for many
//! integration problems, especially in higher dimensions.
//!
//! The module includes several low-discrepancy sequence generators:
//! - **Sobol sequences**: Well-distributed points with good coverage
//! - **Halton sequences**: Based on van der Corput sequences with different prime bases
//! - **Faure sequences**: Based on prime bases with matrix scrambling for improved uniformity

use scirs2_core::ndarray::{s, Array1, Array2, ArrayView1};
use scirs2_core::numeric::{Float, FromPrimitive};
use scirs2_core::random::random;
use std::fmt;

use crate::error::{IntegrateError, IntegrateResult};

/// Result type for QMC integration
#[derive(Clone, Debug)]
pub struct QMCQuadResult<T> {
    /// The estimate of the integral
    pub integral: T,
    /// The error estimate
    pub standard_error: T,
}

impl<T: fmt::Display> fmt::Display for QMCQuadResult<T> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "QMCQuadResult(integral={}, standard_error={})",
            self.integral, self.standard_error
        )
    }
}

/// Trait for quasi-random number generators
pub trait QRNGEngine: Send + Sync {
    /// Generate n points in d dimensions in the unit hypercube `[0,1]^d`
    fn random(&mut self, n: usize) -> Array2<f64>;

    /// Dimensionality of the generator
    fn dim(&self) -> usize;

    /// Create a new instance from a seed
    fn new_from_seed(&self, seed: u64) -> Box<dyn QRNGEngine>;
}

/// Simple pseudorandom number generator for benchmarking
pub struct RandomGenerator {
    dim: usize,
    // We don't use the _seed for random generation, just for new_from_seed
    _seed: u64,
}

impl RandomGenerator {
    /// Create a new random number generator
    pub fn new(dim: usize, seed: Option<u64>) -> Self {
        let _seed = seed.unwrap_or_else(random::<u64>);

        Self { dim, _seed }
    }
}

impl QRNGEngine for RandomGenerator {
    fn random(&mut self, n: usize) -> Array2<f64> {
        let mut result = Array2::zeros((n, self.dim));

        for i in 0..n {
            for j in 0..self.dim {
                result[[i, j]] = random::<f64>();
            }
        }

        result
    }

    fn dim(&self) -> usize {
        self.dim
    }

    fn new_from_seed(&self, seed: u64) -> Box<dyn QRNGEngine> {
        Box::new(Self::new(self.dim, Some(seed)))
    }
}

/// Sobol sequence generator
pub struct Sobol {
    dim: usize,
    seed: u64,
    curr_index: u64,
    direction_numbers: Vec<Vec<u64>>,
    last_point: Vec<u64>,
}

impl Sobol {
    /// Create a new Sobol sequence generator
    pub fn new(dim: usize, seed: Option<u64>) -> Self {
        let seed = seed.unwrap_or_else(random::<u64>);

        let mut sobol = Self {
            dim,
            seed,
            curr_index: 0,
            direction_numbers: Vec::new(),
            last_point: vec![0; dim],
        };

        sobol.initialize_direction_numbers();
        sobol
    }

    /// Initialize direction numbers for Sobol sequence
    fn initialize_direction_numbers(&mut self) {
        self.direction_numbers = vec![Vec::new(); self.dim];

        // First dimension uses powers of 2
        self.direction_numbers[0] = (0..64).map(|i| 1u64 << (63 - i)).collect();

        // For higher dimensions, use primitive polynomials and initial direction numbers
        // This is a simplified set for up to 10 dimensions
        let primitive_polynomials = [
            0,  // dimension 0 (not used)
            0,  // dimension 1 (powers of 2)
            3,  // x + 1
            7,  // x^2 + x + 1
            11, // x^3 + x + 1
            13, // x^3 + x^2 + 1
            19, // x^4 + x + 1
            25, // x^4 + x^3 + 1
            37, // x^5 + x^2 + 1
            41, // x^5 + x^3 + 1
            55, // x^5 + x^4 + x^2 + x + 1
        ];

        let initial_numbers = vec![
            vec![], // dimension 0
            vec![], // dimension 1 (powers of 2)
            vec![1],
            vec![1, 1],
            vec![1, 3, 1],
            vec![1, 1, 3],
            vec![1, 3, 3, 9],
            vec![1, 1, 5, 5],
            vec![1, 3, 1, 13],
            vec![1, 1, 5, 5, 17],
            vec![1, 3, 5, 5, 5],
        ];

        for d in 1..std::cmp::min(self.dim, primitive_polynomials.len()) {
            let poly = primitive_polynomials[d];
            let init_nums = &initial_numbers[d];

            self.direction_numbers[d] = vec![0; 64];

            // Set initial direction numbers
            for (i, &num) in init_nums.iter().enumerate() {
                self.direction_numbers[d][i] = (num as u64) << (63 - i);
            }

            // Generate remaining direction numbers using recurrence relation
            let degree = self.bit_length(poly) - 1;
            for i in degree..64 {
                let mut value = self.direction_numbers[d][i - degree];

                // Apply primitive polynomial recurrence
                let mut poly_temp = poly;
                for j in 1..degree {
                    if poly_temp & 1 == 1 {
                        value ^= self.direction_numbers[d][i - j];
                    }
                    poly_temp >>= 1;
                }

                self.direction_numbers[d][i] = value;
            }
        }

        // For dimensions beyond our predefined set, use van der Corput sequences
        for d in primitive_polynomials.len()..self.dim {
            self.direction_numbers[d] = (0..64)
                .map(|i| {
                    let base = 2 + (d - primitive_polynomials.len()) as u64;
                    self.van_der_corput_direction_number(i, base)
                })
                .collect();
        }
    }

    /// Calculate bit length of a number
    fn bit_length(&self, mut n: u64) -> usize {
        let mut length = 0;
        while n > 0 {
            length += 1;
            n >>= 1;
        }
        length
    }

    /// Generate van der Corput direction number as fallback
    fn van_der_corput_direction_number(&self, i: usize, base: u64) -> u64 {
        if i == 0 {
            1u64 << 63
        } else {
            let mut value = 0u64;
            let mut n = i + 1;
            let mut denom = base;

            while n > 0 && denom <= (1u64 << 63) {
                value |= ((n % base as usize) as u64) << (64 - self.bit_length(denom));
                n /= base as usize;
                denom *= base;
            }

            value
        }
    }

    /// Generate a Sobol sequence point using proper Sobol algorithm
    fn generate_point(&mut self) -> Vec<f64> {
        if self.curr_index == 0 {
            self.curr_index = 1;
            return vec![0.0; self.dim];
        }

        // Find rightmost zero bit in Gray code representation
        let gray_code_index = self.curr_index ^ (self.curr_index >> 1);
        let rightmost_zero = (!gray_code_index).trailing_zeros() as usize;

        // Update the Sobol point
        for d in 0..self.dim {
            if rightmost_zero < self.direction_numbers[d].len() {
                self.last_point[d] ^= self.direction_numbers[d][rightmost_zero];
            }
        }

        self.curr_index += 1;

        // Convert to floating point
        self.last_point
            .iter()
            .map(|&x| (x as f64) / u64::MAX as f64)
            .collect()
    }
}

impl QRNGEngine for Sobol {
    fn random(&mut self, n: usize) -> Array2<f64> {
        let mut result = Array2::zeros((n, self.dim));

        for i in 0..n {
            let point = self.generate_point();
            for j in 0..self.dim {
                result[[i, j]] = point[j];
            }
        }

        result
    }

    fn dim(&self) -> usize {
        self.dim
    }

    fn new_from_seed(&self, seed: u64) -> Box<dyn QRNGEngine> {
        Box::new(Self::new(self.dim, Some(seed)))
    }
}

/// Halton sequence generator
pub struct Halton {
    dim: usize,
    seed: u64,
    curr_index: usize,
}

impl Halton {
    /// Create a new Halton sequence generator
    pub fn new(dim: usize, seed: Option<u64>) -> Self {
        let seed = seed.unwrap_or_else(random::<u64>);

        Self {
            dim,
            seed,
            curr_index: 0,
        }
    }

    /// Generate a Halton sequence point
    fn generate_point(&mut self) -> Vec<f64> {
        let first_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53];
        let mut result = vec![0.0; self.dim];

        // Scrambling offset using seed
        let offset = (self.seed % 1000) as usize;

        for d in 0..self.dim {
            let base = if d < first_primes.len() {
                first_primes[d]
            } else {
                // For dimensions beyond stored primes, use a simple formula
                first_primes[d % first_primes.len()] + d
            };

            let mut f = 1.0;
            let mut r = 0.0;
            let mut n = self.curr_index + offset;

            while n > 0 {
                f /= base as f64;
                r += f * (n % base) as f64;
                n /= base;
            }

            result[d] = r;
        }

        self.curr_index += 1;
        result
    }
}

impl QRNGEngine for Halton {
    fn random(&mut self, n: usize) -> Array2<f64> {
        let mut result = Array2::zeros((n, self.dim));

        for i in 0..n {
            let point = self.generate_point();
            for j in 0..self.dim {
                result[[i, j]] = point[j];
            }
        }

        result
    }

    fn dim(&self) -> usize {
        self.dim
    }

    fn new_from_seed(&self, seed: u64) -> Box<dyn QRNGEngine> {
        Box::new(Self::new(self.dim, Some(seed)))
    }
}

/// Faure sequence generator
pub struct Faure {
    dim: usize,
    seed: u64,
    curr_index: usize,
}

impl Faure {
    /// Create a new Faure sequence generator
    pub fn new(dim: usize, seed: Option<u64>) -> Self {
        let seed = seed.unwrap_or_else(random::<u64>);

        Self {
            dim,
            seed,
            curr_index: 0,
        }
    }

    /// Generate a Faure sequence point
    ///
    /// The Faure sequence is based on a prime base p >= dim, where p is the smallest
    /// prime number greater than or equal to the dimension.
    fn generate_point(&mut self) -> Vec<f64> {
        let base = self.find_prime_base(self.dim);
        let mut result = vec![0.0; self.dim];

        // Scrambling offset using seed
        let offset = (self.seed % 1000) as usize;
        let scrambled_index = self.curr_index + offset;

        // Generate coordinates using Faure matrices
        for (d, result_elem) in result.iter_mut().enumerate().take(self.dim) {
            *result_elem = self.faure_coordinate(scrambled_index, d, base);
        }

        self.curr_index += 1;
        result
    }

    /// Find the smallest prime >= n
    fn find_prime_base(&self, n: usize) -> usize {
        if n <= 2 {
            return 2;
        }

        let mut candidate = n;
        while !self.is_prime(candidate) {
            candidate += 1;
        }
        candidate
    }

    /// Simple primality test
    fn is_prime(&self, n: usize) -> bool {
        if n < 2 {
            return false;
        }
        if n == 2 {
            return true;
        }
        if n.is_multiple_of(2) {
            return false;
        }

        let sqrt_n = (n as f64).sqrt() as usize;
        for i in (3..=sqrt_n).step_by(2) {
            if n.is_multiple_of(i) {
                return false;
            }
        }
        true
    }

    /// Generate d-th coordinate using Faure construction
    fn faure_coordinate(&self, index: usize, dimension: usize, base: usize) -> f64 {
        if index == 0 {
            return 0.0;
        }

        // Convert index to base-p representation and apply Faure matrix
        let digits = self.to_base_digits(index, base);
        let mut result = 0.0;
        let mut base_power = base as f64;

        // Apply Faure matrix transformation for dimension d
        for (i, &digit) in digits.iter().enumerate() {
            let transformed_digit = self.faure_matrix_element(i, dimension, base) * digit;
            result += (transformed_digit % base) as f64 / base_power;
            base_power *= base as f64;
        }

        result.fract() // Ensure result is in [0, 1)
    }

    /// Convert integer to base-p digits (least significant first)
    fn to_base_digits(&self, mut n: usize, base: usize) -> Vec<usize> {
        if n == 0 {
            return vec![0];
        }

        let mut digits = Vec::new();
        while n > 0 {
            digits.push(n % base);
            n /= base;
        }
        digits
    }

    /// Simplified Faure matrix element (for educational implementation)
    /// In practice, this would use precomputed Pascal triangle modulo p
    fn faure_matrix_element(&self, i: usize, dimension: usize, base: usize) -> usize {
        // Simplified implementation using binomial coefficients mod p
        // For a proper implementation, use Lucas' theorem and Pascal's triangle
        if dimension == 0 {
            if i == 0 {
                1
            } else {
                0
            }
        } else {
            // Approximate with a simple formula for demonstration
            let offset = dimension * 7 + self.seed as usize % 13; // Add some scrambling
            ((i + offset + 1) % base + 1) % base
        }
    }
}

impl QRNGEngine for Faure {
    fn random(&mut self, n: usize) -> Array2<f64> {
        let mut result = Array2::zeros((n, self.dim));

        for i in 0..n {
            let point = self.generate_point();
            for j in 0..self.dim {
                result[[i, j]] = point[j];
            }
        }

        result
    }

    fn dim(&self) -> usize {
        self.dim
    }

    fn new_from_seed(&self, seed: u64) -> Box<dyn QRNGEngine> {
        Box::new(Self::new(self.dim, Some(seed)))
    }
}

/// Scale samples from unit hypercube to the integration range [a, b]
#[allow(dead_code)]
pub fn scale<T: Float + FromPrimitive>(
    sample: &Array2<T>,
    a: &Array1<T>,
    b: &Array1<T>,
) -> Array2<T> {
    let mut scaled = Array2::<T>::zeros(sample.raw_dim());
    let dim = sample.shape()[1];

    for i in 0..sample.shape()[0] {
        for j in 0..dim {
            scaled[[i, j]] = a[j] + (b[j] - a[j]) * sample[[i, j]];
        }
    }

    scaled
}

/// Compute an integral in N dimensions using Quasi-Monte Carlo quadrature.
///
/// # Parameters
///
/// * `func` - Function to integrate that takes a point in n-dimensional space and
///   returns a scalar value.
/// * `a` - Lower integration bounds, array of length equal to number of dimensions.
/// * `b` - Upper integration bounds, array of length equal to number of dimensions.
/// * `n_estimates` - Number of statistically independent samples to use (default: 8).
/// * `n_points` - Number of QMC points per sample (default: 1024).
/// * `qrng` - QRNGEngine to use for sampling points. If None, a Halton sequence is used.
/// * `log` - If true, treat func as returning the log of the integrand and return log
///   of the result.
///
/// # Returns
///
/// * `QMCQuadResult` containing the integral estimate and standard error.
///
/// # Examples
///
/// ```
/// use scirs2_core::ndarray::{Array1, ArrayView1};
/// use scirs2_integrate::qmc::{qmc_quad, Halton};
///
/// let f = |x: ArrayView1<f64>| x[0].powi(2) * x[1].exp();
/// let a = Array1::from_vec(vec![0.0, 0.0]);
/// let b = Array1::from_vec(vec![1.0, 1.0]);
///
/// let qrng = Halton::new(2, Some(42));
/// let result = qmc_quad(f, &a, &b, None, None, Some(Box::new(qrng)), false).expect("Operation failed");
/// println!("Integral: {}, Error: {}", result.integral, result.standard_error);
/// ```
#[allow(dead_code)]
pub fn qmc_quad<F>(
    func: F,
    a: &Array1<f64>,
    b: &Array1<f64>,
    n_estimates: Option<usize>,
    n_points: Option<usize>,
    qrng: Option<Box<dyn QRNGEngine>>,
    log: bool,
) -> IntegrateResult<QMCQuadResult<f64>>
where
    F: Fn(ArrayView1<f64>) -> f64,
{
    let n_estimates = n_estimates.unwrap_or(8);
    let n_points = n_points.unwrap_or(1024);

    // Input validation
    if a.len() != b.len() {
        return Err(IntegrateError::ValueError(
            "Dimension mismatch: 'a' and 'b' must have the same length".to_string(),
        ));
    }

    let dim = a.len();

    // Initialize QRNG if not provided
    let mut qrng = qrng.unwrap_or_else(|| Box::new(Halton::new(dim, None)));

    if qrng.dim() != dim {
        return Err(IntegrateError::ValueError(format!(
            "QRNG dimension ({}) does not match integration dimension ({})",
            qrng.dim(),
            dim
        )));
    }

    // Check if a = b for any dimension, return 0 if so
    for i in 0..dim {
        if (a[i] - b[i]).abs() < f64::EPSILON {
            return Ok(QMCQuadResult {
                integral: if log { f64::NEG_INFINITY } else { 0.0 },
                standard_error: 0.0,
            });
        }
    }

    // Swap limits if a > b and record the sign change
    let mut a_mod = a.clone();
    let mut b_mod = b.clone();
    let mut sign = 1.0;

    for i in 0..dim {
        if a[i] > b[i] {
            a_mod[i] = b[i];
            b_mod[i] = a[i];
            sign *= -1.0;
        }
    }

    // Volume of the hypercube
    let volume = (0..dim).map(|i| b_mod[i] - a_mod[i]).product::<f64>();
    let delta = volume / (n_points as f64);

    // Prepare for multiple _estimates
    let mut estimates = Array1::<f64>::zeros(n_estimates);

    // Generate independent samples and compute _estimates
    for i in 0..n_estimates {
        // Generate QMC sample
        let sample = qrng.random(n_points);

        // Scale to integration domain
        let x = scale(&sample, &a_mod, &b_mod);

        // Evaluate function at sample _points
        let mut sum = 0.0;

        if log {
            let mut max_val = f64::NEG_INFINITY;
            let mut log_values = Vec::with_capacity(n_points);

            for j in 0..n_points {
                let val = func(x.slice(s![j, ..]));
                log_values.push(val);
                if val > max_val {
                    max_val = val;
                }
            }

            // Compute log sum exp
            let mut sum_exp = 0.0;
            for val in log_values {
                sum_exp += (val - max_val).exp();
            }

            estimates[i] = max_val + sum_exp.ln() + delta.ln();
        } else {
            for j in 0..n_points {
                sum += func(x.slice(s![j, ..]));
            }

            estimates[i] = sum * delta;
        }

        // Get a new QRNG for next estimate with different scrambling
        let seed = i as u64 + 1;
        qrng = qrng.new_from_seed(seed);
    }

    // Compute final estimate and error
    let integral = if log {
        let mut max_val = f64::NEG_INFINITY;
        for i in 0..n_estimates {
            if estimates[i] > max_val {
                max_val = estimates[i];
            }
        }

        let mut sum_exp = 0.0;
        for i in 0..n_estimates {
            sum_exp += (estimates[i] - max_val).exp();
        }

        max_val + (sum_exp / (n_estimates as f64)).ln()
    } else {
        estimates.sum() / (n_estimates as f64)
    };

    // Compute standard error
    let standard_error = if n_estimates > 1 {
        if log {
            // For log space, compute standard error differently
            let mean = integral;
            let mut variance = 0.0;

            for i in 0..n_estimates {
                let diff = (estimates[i] - mean).exp();
                variance += (diff - 1.0).powi(2);
            }

            variance /= (n_estimates - 1) as f64;
            (variance / (n_estimates as f64)).sqrt().ln()
        } else {
            let mean = integral;
            let mut variance = 0.0;

            for estimate in estimates.iter().take(n_estimates) {
                variance += (estimate - mean).powi(2);
            }

            variance /= (n_estimates - 1) as f64;
            (variance / (n_estimates as f64)).sqrt()
        }
    } else if log {
        f64::NEG_INFINITY
    } else {
        0.0
    };

    // Apply sign correction for reversed limits
    let final_integral = if log && sign < 0.0 {
        // For negative results in log space, we'd need complex numbers
        // Since we don't support complex results, we'll just negate the result
        // and warn about it in the documentation
        -integral
    } else {
        integral * sign
    };

    Ok(QMCQuadResult {
        integral: final_integral,
        standard_error,
    })
}

/// Parallel Quasi-Monte Carlo integration with workers parameter
///
/// This function provides the same functionality as `qmc_quad` but with
/// explicit control over the number of worker threads for parallel evaluation.
///
/// # Arguments
///
/// * `func` - The function to integrate
/// * `a` - Lower bounds of integration for each dimension
/// * `b` - Upper bounds of integration for each dimension
/// * `n_estimates` - Number of independent estimates to compute (default: 8)
/// * `n_points` - Number of sample points per estimate (default: 1024)
/// * `qrng` - QRNGEngine to use for sampling points. If None, a Halton sequence is used.
/// * `log` - If true, treat func as returning the log of the integrand and return log
///   of the result.
/// * `workers` - Number of worker threads to use. If None, uses all available cores.
///
/// # Returns
///
/// * `QMCQuadResult` containing the integral estimate and standard error.
///
/// # Examples
///
/// ```
/// use scirs2_integrate::qmc::{qmc_quad_parallel, Halton};
/// use scirs2_core::ndarray::{Array1, ArrayView1};
///
/// let f = |x: ArrayView1<f64>| x[0].powi(2) * x[1].exp();
/// let a = Array1::from_vec(vec![0.0, 0.0]);
/// let b = Array1::from_vec(vec![1.0, 1.0]);
///
/// let qrng = Halton::new(2, Some(42));
/// let result = qmc_quad_parallel(
///     f, &a, &b, None, None, Some(Box::new(qrng)), false, Some(4)
/// ).expect("Operation failed");
/// println!("Integral: {}, Error: {}", result.integral, result.standard_error);
/// ```
#[allow(dead_code)]
pub fn qmc_quad_parallel<F>(
    func: F,
    a: &Array1<f64>,
    b: &Array1<f64>,
    n_estimates: Option<usize>,
    n_points: Option<usize>,
    qrng: Option<Box<dyn QRNGEngine>>,
    log: bool,
    workers: Option<usize>,
) -> IntegrateResult<QMCQuadResult<f64>>
where
    F: Fn(ArrayView1<f64>) -> f64 + Send + Sync,
{
    let n_estimates = n_estimates.unwrap_or(8);
    let n_points = n_points.unwrap_or(1024);

    // Input validation
    if a.len() != b.len() {
        return Err(IntegrateError::ValueError(
            "Dimension mismatch: 'a' and 'b' must have the same length".to_string(),
        ));
    }

    let dim = a.len();

    // Initialize QRNG if not provided
    let qrng = qrng.unwrap_or_else(|| Box::new(Halton::new(dim, None)));

    if qrng.dim() != dim {
        return Err(IntegrateError::ValueError(format!(
            "QRNG dimension ({}) does not match integration dimension ({})",
            qrng.dim(),
            dim
        )));
    }

    // Check if a = b for any dimension, return 0 if so
    for i in 0..dim {
        if (a[i] - b[i]).abs() < f64::EPSILON {
            return Ok(QMCQuadResult {
                integral: if log { f64::NEG_INFINITY } else { 0.0 },
                standard_error: 0.0,
            });
        }
    }

    // Swap limits if a > b and record the sign change
    let mut a_mod = a.clone();
    let mut b_mod = b.clone();
    let mut sign = 1.0;

    for i in 0..dim {
        if a[i] > b[i] {
            a_mod[i] = b[i];
            b_mod[i] = a[i];
            sign *= -1.0;
        }
    }

    // Calculate domain volume
    let mut volume = 1.0;
    for i in 0..dim {
        volume *= b_mod[i] - a_mod[i];
    }

    // Configure parallel execution based on workers parameter
    #[cfg(feature = "parallel")]
    {
        if let Some(num_workers) = workers {
            // Set thread pool size
            use scirs2_core::parallel_ops::ThreadPoolBuilder;
            let pool = ThreadPoolBuilder::new()
                .num_threads(num_workers)
                .build()
                .map_err(|_| {
                    IntegrateError::ValueError("Failed to create thread pool".to_string())
                })?;

            pool.install(|| {
                parallel_qmc_integration_impl(
                    func,
                    &a_mod,
                    &b_mod,
                    n_estimates,
                    n_points,
                    &*qrng,
                    log,
                    volume,
                    sign,
                )
            })
        } else {
            // Use default parallel execution
            parallel_qmc_integration_impl(
                func,
                &a_mod,
                &b_mod,
                n_estimates,
                n_points,
                &*qrng,
                log,
                volume,
                sign,
            )
        }
    }

    #[cfg(not(feature = "parallel"))]
    {
        // If parallel feature is not enabled, fall back to sequential execution
        let _ = workers; // Silence unused variable warning
        sequential_qmc_integration(
            func,
            &a_mod,
            &b_mod,
            n_estimates,
            n_points,
            qrng,
            log,
            volume,
            sign,
        )
    }
}

#[cfg(feature = "parallel")]
#[allow(dead_code)]
fn parallel_qmc_integration_impl<F>(
    func: F,
    a: &Array1<f64>,
    b: &Array1<f64>,
    n_estimates: usize,
    n_points: usize,
    qrng: &dyn QRNGEngine,
    log: bool,
    volume: f64,
    sign: f64,
) -> IntegrateResult<QMCQuadResult<f64>>
where
    F: Fn(ArrayView1<f64>) -> f64 + Send + Sync,
{
    use scirs2_core::parallel_ops::*;

    // Generate _estimates in parallel
    let _estimates: Vec<f64> = (0..n_estimates)
        .into_par_iter()
        .map(|_| {
            // Each thread gets its own QRNG instance with different seed
            let mut local_qrng = qrng.new_from_seed(scirs2_core::random::random());

            // Sample points
            let points = local_qrng.random(n_points);

            // Transform points to integration domain and evaluate function
            let mut sum = 0.0;

            for i in 0..n_points {
                let point = points.row(i);

                // Transform from [0,1]^d to [a,b]^d
                let mut transformed_point = Array1::zeros(a.len());
                for j in 0..a.len() {
                    transformed_point[j] = a[j] + point[j] * (b[j] - a[j]);
                }

                let value = func(transformed_point.view());

                if log {
                    // For log-space integration, we accumulate log values
                    if value > f64::NEG_INFINITY {
                        sum += value.exp() / n_points as f64;
                    }
                } else {
                    sum += value / n_points as f64;
                }
            }

            if log {
                sum.ln() + volume.ln()
            } else {
                sum * volume
            }
        })
        .collect();

    compute_qmc_result(_estimates, log, sign)
}

#[cfg(not(feature = "parallel"))]
#[allow(dead_code)]
fn sequential_qmc_integration<F>(
    func: F,
    a: &Array1<f64>,
    b: &Array1<f64>,
    n_estimates: usize,
    n_points: usize,
    mut qrng: Box<dyn QRNGEngine>,
    log: bool,
    volume: f64,
    sign: f64,
) -> IntegrateResult<QMCQuadResult<f64>>
where
    F: Fn(ArrayView1<f64>) -> f64,
{
    let mut estimates = Vec::with_capacity(n_estimates);

    for _ in 0..n_estimates {
        // Sample _points
        let points = qrng.random(n_points);

        // Transform _points to integration domain and evaluate function
        let mut sum = 0.0;

        for i in 0..n_points {
            let point = points.row(i);

            // Transform from [0,1]^d to [a,b]^d
            let mut transformed_point = Array1::zeros(a.len());
            for j in 0..a.len() {
                transformed_point[j] = a[j] + point[j] * (b[j] - a[j]);
            }

            let value = func(transformed_point.view());

            if log {
                // For log-space integration, we accumulate log values
                if value > f64::NEG_INFINITY {
                    sum += value.exp() / n_points as f64;
                }
            } else {
                sum += value / n_points as f64;
            }
        }

        let estimate = if log {
            sum.ln() + volume.ln()
        } else {
            sum * volume
        };

        estimates.push(estimate);
    }

    compute_qmc_result(estimates, log, sign)
}

#[allow(dead_code)]
fn compute_qmc_result(
    estimates: Vec<f64>,
    log: bool,
    sign: f64,
) -> IntegrateResult<QMCQuadResult<f64>> {
    let n_estimates = estimates.len();

    // Compute mean and variance
    let integral = if log {
        // For log-space, use log-sum-exp for numerical stability
        let max_val = estimates.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
        if max_val == f64::NEG_INFINITY {
            f64::NEG_INFINITY
        } else {
            let sum_exp: f64 = estimates.iter().map(|&x| (x - max_val).exp()).sum();
            max_val + (sum_exp / n_estimates as f64).ln()
        }
    } else {
        estimates.iter().sum::<f64>() / n_estimates as f64
    };

    let standard_error = if estimates.len() > 1 {
        if log {
            // For log-space, compute variance in log domain
            let mean = integral;
            let mut variance = 0.0;

            for estimate in estimates.iter().take(n_estimates) {
                let diff = estimate - mean;
                variance += diff * diff;
            }

            variance /= (n_estimates - 1) as f64;
            (variance / (n_estimates as f64)).sqrt().ln()
        } else {
            let mean = integral;
            let mut variance = 0.0;

            for estimate in estimates.iter().take(n_estimates) {
                variance += (estimate - mean).powi(2);
            }

            variance /= (n_estimates - 1) as f64;
            (variance / (n_estimates as f64)).sqrt()
        }
    } else if log {
        f64::NEG_INFINITY
    } else {
        0.0
    };

    // Apply sign correction for reversed limits
    let final_integral = if log && sign < 0.0 {
        // For negative results in log space, we'd need complex numbers
        // Since we don't support complex results, we'll just negate the result
        -integral
    } else {
        integral * sign
    };

    Ok(QMCQuadResult {
        integral: final_integral,
        standard_error,
    })
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{qmc_quad, qmc_quad_parallel, Faure};
    use approx::assert_abs_diff_eq;
    use scirs2_core::ndarray::{s, Array1, ArrayView1};

    #[test]
    fn test_qmc_integral_1d() {
        // Test integral of x^2 from 0 to 1 (should be 1/3)
        let f = |x: ArrayView1<f64>| x[0].powi(2);
        let a = Array1::from_vec(vec![0.0]);
        let b = Array1::from_vec(vec![1.0]);

        let result =
            qmc_quad(f, &a, &b, Some(8), Some(1000), None, false).expect("Operation failed");

        assert_abs_diff_eq!(result.integral, 1.0 / 3.0, epsilon = 0.01);
    }

    #[test]
    fn test_qmc_integral_2d() {
        // Test integral of x*y from 0 to 1, 0 to 1 (should be 1/4)
        let f = |x: ArrayView1<f64>| x[0] * x[1];
        let a = Array1::from_vec(vec![0.0, 0.0]);
        let b = Array1::from_vec(vec![1.0, 1.0]);

        let result =
            qmc_quad(f, &a, &b, Some(8), Some(1000), None, false).expect("Operation failed");

        assert_abs_diff_eq!(result.integral, 0.25, epsilon = 0.01);
    }

    #[test]
    fn test_infinite_limits() {
        // Test integral of e^(-x^2) from -∞ to ∞ (should be sqrt(π) ≈ 1.77)
        let f = |x: ArrayView1<f64>| (-x[0].powi(2)).exp();
        let a = Array1::from_vec(vec![-5.0]); // approximating infinity with a large value
        let b = Array1::from_vec(vec![5.0]);

        let result =
            qmc_quad(f, &a, &b, Some(8), Some(1000), None, false).expect("Operation failed");

        assert_abs_diff_eq!(result.integral, std::f64::consts::PI.sqrt(), epsilon = 0.05);
    }

    #[test]
    fn test_faure_sequence() {
        // Test basic Faure sequence properties
        let mut faure = Faure::new(2, Some(42));

        // Generate some points
        let points = faure.random(10);

        // Check dimensions
        assert_eq!(points.shape()[0], 10);
        assert_eq!(points.shape()[1], 2);

        // Check that all points are in [0, 1)
        for i in 0..10 {
            for j in 0..2 {
                assert!(points[[i, j]] >= 0.0 && points[[i, j]] < 1.0);
            }
        }

        // Test reproducibility with same seed
        let mut faure2 = Faure::new(2, Some(42));
        let points2 = faure2.random(5);
        let points_first_5 = points.slice(s![0..5, ..]);

        for i in 0..5 {
            for j in 0..2 {
                assert_abs_diff_eq!(points_first_5[[i, j]], points2[[i, j]], epsilon = 1e-10);
            }
        }
    }

    #[test]
    fn test_faure_integration() {
        // Test integral using Faure sequence
        let f = |x: ArrayView1<f64>| x[0] * x[1];
        let a = Array1::from_vec(vec![0.0, 0.0]);
        let b = Array1::from_vec(vec![1.0, 1.0]);

        let faure = Faure::new(2, Some(42));
        let result = qmc_quad(f, &a, &b, Some(8), Some(1000), Some(Box::new(faure)), false)
            .expect("Operation failed");

        // Expected result is 1/4 = 0.25
        // Note: This simplified Faure implementation may not achieve optimal convergence
        // Production implementations would use proper Pascal triangle modulo p matrices
        assert_abs_diff_eq!(result.integral, 0.25, epsilon = 0.2);
    }

    #[test]
    fn test_qmc_quad_parallel_workers() {
        // Test QMC integration with workers parameter
        let f = |x: ArrayView1<f64>| x[0].powi(2) + x[1].powi(2);
        let a = Array1::from_vec(vec![0.0, 0.0]);
        let b = Array1::from_vec(vec![1.0, 1.0]);

        // Test with 2 workers
        let result = qmc_quad_parallel(f, &a, &b, Some(8), Some(1000), None, false, Some(2))
            .expect("Operation failed");

        // Expected integral of x^2 + y^2 over [0,1]×[0,1] is 2/3
        assert_abs_diff_eq!(result.integral, 2.0 / 3.0, epsilon = 0.1);
        assert!(result.standard_error >= 0.0);
    }
}