quantrs2_device/dynamical_decoupling/

fallback_scirs2.rs

//! Fallback implementations for SciRS2 functions when the feature is not available

use scirs2_core::ndarray::{Array1, Array2, ArrayView1, ArrayView2};
use std::collections::HashMap;

/// Fallback optimization result
#[derive(Debug, Clone)]
pub struct OptimizeResult {
    /// Optimal parameters
    pub x: Array1<f64>,
    /// Optimal function value
    pub fun: f64,
    /// Number of iterations
    pub nit: usize,
    /// Number of function evaluations
    pub nfev: usize,
    /// Success flag
    pub success: bool,
    /// Status message
    pub message: String,
}

/// Fallback minimize function
pub fn minimize<F>(
    objective: F,
    initial: &Array1<f64>,
    _method: &str,
) -> Result<OptimizeResult, String>
where
    F: Fn(&Array1<f64>) -> f64,
{
    // Simple gradient-free optimization using Nelder-Mead-like approach
    let mut current_x = initial.clone();
    let mut current_f = objective(&current_x);
    let n_params = initial.len();

    let mut step_size = 0.1;
    let max_iterations = 1000;
    let tolerance = 1e-6;

    let mut nfev = 1;

    for _iteration in 0..max_iterations {
        let mut improved = false;

        // Try step in each parameter direction
        for i in 0..n_params {
            // Positive step
            let mut test_x = current_x.clone();
            test_x[i] += step_size;
            let test_f = objective(&test_x);
            nfev += 1;

            if test_f < current_f {
                current_x = test_x;
                current_f = test_f;
                improved = true;
                continue;
            }

            // Negative step
            let mut test_x = current_x.clone();
            test_x[i] -= step_size;
            let test_f = objective(&test_x);
            nfev += 1;

            if test_f < current_f {
                current_x = test_x;
                current_f = test_f;
                improved = true;
            }
        }

        if !improved {
            // Reduce step size
            step_size *= 0.5;
            if step_size < tolerance {
                break;
            }
        }
    }

    Ok(OptimizeResult {
        x: current_x,
        fun: current_f,
        nit: max_iterations,
        nfev,
        success: true,
        message: "Optimization completed".to_string(),
    })
}

/// Fallback statistical functions
pub fn mean(data: &ArrayView1<f64>) -> Result<f64, String> {
    if data.is_empty() {
        return Err("Cannot compute mean of empty array".to_string());
    }
    Ok(data.sum() / data.len() as f64)
}

pub fn std(data: &ArrayView1<f64>, ddof: i32, _workers: Option<usize>) -> Result<f64, String> {
    if data.len() <= ddof as usize {
        return Err("Insufficient data for standard deviation calculation".to_string());
    }

    let mean_val = mean(data)?;
    let variance = data.iter().map(|x| (x - mean_val).powi(2)).sum::<f64>()
        / (data.len() as f64 - ddof as f64);

    Ok(variance.sqrt())
}

pub fn var(data: &ArrayView1<f64>, ddof: i32, _workers: Option<usize>) -> Result<f64, String> {
    if data.len() <= ddof as usize {
        return Err("Insufficient data for variance calculation".to_string());
    }

    let mean_val = mean(data)?;
    let variance = data.iter().map(|x| (x - mean_val).powi(2)).sum::<f64>()
        / (data.len() as f64 - ddof as f64);

    Ok(variance)
}

pub fn pearsonr(x: &ArrayView1<f64>, y: &ArrayView1<f64>) -> Result<f64, String> {
    if x.len() != y.len() {
        return Err("Arrays must have same length".to_string());
    }

    if x.len() < 2 {
        return Err("Need at least 2 data points".to_string());
    }

    let mean_x = mean(x)?;
    let mean_y = mean(y)?;

    let numerator: f64 = x
        .iter()
        .zip(y.iter())
        .map(|(&xi, &yi)| (xi - mean_x) * (yi - mean_y))
        .sum();

    let sum_sq_x: f64 = x.iter().map(|&xi| (xi - mean_x).powi(2)).sum();
    let sum_sq_y: f64 = y.iter().map(|&yi| (yi - mean_y).powi(2)).sum();

    let denominator = (sum_sq_x * sum_sq_y).sqrt();

    if denominator == 0.0 {
        Ok(0.0)
    } else {
        Ok(numerator / denominator)
    }
}

pub fn spearmanr(x: &ArrayView1<f64>, y: &ArrayView1<f64>) -> Result<f64, String> {
    if x.len() != y.len() {
        return Err("Arrays must have same length".to_string());
    }

    // Convert to ranks and compute Pearson correlation of ranks
    let x_ranks = rank_array(x);
    let y_ranks = rank_array(y);

    let x_ranks_view = x_ranks.view();
    let y_ranks_view = y_ranks.view();

    pearsonr(&x_ranks_view, &y_ranks_view)
}

fn rank_array(data: &ArrayView1<f64>) -> Array1<f64> {
    let mut indexed_data: Vec<(usize, f64)> =
        data.iter().enumerate().map(|(i, &x)| (i, x)).collect();
    indexed_data.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));

    let mut ranks = Array1::zeros(data.len());
    for (rank, &(index, _)) in indexed_data.iter().enumerate() {
        ranks[index] = rank as f64 + 1.0;
    }

    ranks
}

pub fn ttest_1samp(data: &ArrayView1<f64>, pop_mean: f64) -> Result<(f64, f64), String> {
    if data.len() < 2 {
        return Err("Need at least 2 data points for t-test".to_string());
    }

    let sample_mean = mean(data)?;
    let sample_std = std(data, 1, None)?;
    let n = data.len() as f64;

    let t_statistic = (sample_mean - pop_mean) / (sample_std / n.sqrt());

    // Simplified p-value calculation (assuming normal distribution)
    let df = n - 1.0;
    let p_value = 2.0 * (1.0 - normal_cdf(t_statistic.abs()));

    Ok((t_statistic, p_value))
}

pub fn ks_2samp(x: &ArrayView1<f64>, y: &ArrayView1<f64>) -> Result<(f64, f64), String> {
    if x.is_empty() || y.is_empty() {
        return Err("Both samples must be non-empty".to_string());
    }

    let mut x_sorted = x.to_vec();
    let mut y_sorted = y.to_vec();
    x_sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    y_sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

    let mut all_values = x_sorted.clone();
    all_values.extend_from_slice(&y_sorted);
    all_values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    all_values.dedup();

    let mut max_diff = 0.0f64;

    for &value in &all_values {
        let cdf_x = x_sorted.iter().filter(|&&x| x <= value).count() as f64 / x_sorted.len() as f64;
        let cdf_y = y_sorted.iter().filter(|&&y| y <= value).count() as f64 / y_sorted.len() as f64;
        let diff = (cdf_x - cdf_y).abs();
        max_diff = max_diff.max(diff);
    }

    // Simplified p-value calculation
    let n_x = x.len() as f64;
    let n_y = y.len() as f64;
    let sqrt_term = ((n_x + n_y) / (n_x * n_y)).sqrt();
    let ks_statistic = max_diff;
    let p_value = 2.0f64 * (-2.0f64 * ks_statistic.powi(2) / sqrt_term.powi(2)).exp();

    Ok((ks_statistic, p_value.min(1.0)))
}

pub fn shapiro_wilk(data: &ArrayView1<f64>) -> Result<(f64, f64), String> {
    if data.len() < 3 || data.len() > 5000 {
        return Err("Shapiro-Wilk test requires 3-5000 observations".to_string());
    }

    // Simplified implementation - just check if data looks roughly normal
    let mean_val = mean(data)?;
    let std_val = std(data, 1, None)?;

    // Calculate skewness and kurtosis as rough normality indicators
    let n = data.len() as f64;
    let skewness = data
        .iter()
        .map(|&x| ((x - mean_val) / std_val).powi(3))
        .sum::<f64>()
        / n;

    let kurtosis = data
        .iter()
        .map(|&x| ((x - mean_val) / std_val).powi(4))
        .sum::<f64>()
        / n
        - 3.0;

    // Simple heuristic for W statistic
    let w_statistic = 1.0 - (skewness.abs() + kurtosis.abs()) / 10.0;
    let w_statistic = w_statistic.clamp(0.0, 1.0);

    // Simple p-value based on W statistic
    let p_value = if w_statistic > 0.95 {
        0.5
    } else if w_statistic > 0.9 {
        0.1
    } else {
        0.01
    };

    Ok((w_statistic, p_value))
}

pub fn percentile(data: &ArrayView1<f64>, percentile: f64) -> Result<f64, String> {
    if data.is_empty() {
        return Err("Cannot compute percentile of empty array".to_string());
    }

    if !(0.0..=100.0).contains(&percentile) {
        return Err("Percentile must be between 0 and 100".to_string());
    }

    let mut sorted_data = data.to_vec();
    sorted_data.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

    let index = (percentile / 100.0) * (sorted_data.len() - 1) as f64;
    let lower_index = index.floor() as usize;
    let upper_index = index.ceil() as usize;

    if lower_index == upper_index {
        Ok(sorted_data[lower_index])
    } else {
        let weight = index - lower_index as f64;
        Ok(sorted_data[lower_index].mul_add(1.0 - weight, sorted_data[upper_index] * weight))
    }
}

pub fn trace(matrix: &ArrayView2<f64>) -> Result<f64, String> {
    let (rows, cols) = matrix.dim();
    if rows != cols {
        return Err("Matrix must be square for trace calculation".to_string());
    }

    let mut trace_val = 0.0;
    for i in 0..rows {
        trace_val += matrix[(i, i)];
    }

    Ok(trace_val)
}

/// Real matrix inverse via Gauss-Jordan elimination with partial pivoting
/// (pure Rust fallback for when the `scirs2` feature is disabled). Returns an
/// honest error on non-square or singular matrices.
pub fn inv(matrix: &ArrayView2<f64>) -> Result<Array2<f64>, String> {
    let (rows, cols) = matrix.dim();
    if rows != cols {
        return Err("Matrix must be square for inversion".to_string());
    }
    let n = rows;
    if n == 0 {
        return Ok(Array2::zeros((0, 0)));
    }

    // Working copy of A and the identity that becomes A^-1.
    let mut a = matrix.to_owned();
    let mut inverse = Array2::<f64>::eye(n);

    for col in 0..n {
        // Partial pivot: pick the row with the largest magnitude in this column.
        let mut pivot_row = col;
        let mut pivot_val = a[(col, col)].abs();
        for r in (col + 1)..n {
            let v = a[(r, col)].abs();
            if v > pivot_val {
                pivot_val = v;
                pivot_row = r;
            }
        }
        if pivot_val <= f64::MIN_POSITIVE {
            return Err("singular matrix".to_string());
        }
        if pivot_row != col {
            for c in 0..n {
                a.swap((col, c), (pivot_row, c));
                inverse.swap((col, c), (pivot_row, c));
            }
        }

        // Scale the pivot row so the pivot becomes 1.
        let pivot = a[(col, col)];
        for c in 0..n {
            a[(col, c)] /= pivot;
            inverse[(col, c)] /= pivot;
        }

        // Eliminate the pivot column from every other row.
        for r in 0..n {
            if r == col {
                continue;
            }
            let factor = a[(r, col)];
            if factor == 0.0 {
                continue;
            }
            for c in 0..n {
                let a_val = a[(col, c)];
                let inv_val = inverse[(col, c)];
                a[(r, c)] -= factor * a_val;
                inverse[(r, c)] -= factor * inv_val;
            }
        }
    }

    Ok(inverse)
}

// Helper function for normal CDF approximation
fn normal_cdf(x: f64) -> f64 {
    // Simplified normal CDF approximation
    0.5 * (1.0 + erf(x / 2.0_f64.sqrt()))
}

// Simplified error function approximation
fn erf(x: f64) -> f64 {
    // Abramowitz and Stegun approximation
    let a1 = 0.254_829_592;
    let a2 = -0.284_496_736;
    let a3 = 1.421_413_741;
    let a4 = -1.453_152_027;
    let a5 = 1.061_405_429;
    let p = 0.327_591_1;

    let sign = if x < 0.0 { -1.0 } else { 1.0 };
    let x = x.abs();

    let t = 1.0 / (1.0 + p * x);
    let y = ((a5 * t + a4).mul_add(t, a3).mul_add(t, a2).mul_add(t, a1) * t)
        .mul_add(-(-x * x).exp(), 1.0);

    sign * y
}

#[cfg(test)]
mod tests {
    use super::*;
    use scirs2_core::ndarray::Array1;

    #[test]
    fn test_mean() {
        let data = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
        let result = mean(&data.view()).expect("Mean calculation should succeed");
        assert!((result - 3.0).abs() < 1e-10);
    }

    #[test]
    fn test_std() {
        let data = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
        let result = std(&data.view(), 1, None).expect("Std calculation should succeed");
        // Standard deviation of [1,2,3,4,5] with ddof=1 is sqrt(2.5) ≈ 1.58
        assert!((result - 1.5811388300841898).abs() < 1e-10);
    }

    #[test]
    fn test_pearsonr() {
        let x = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
        let y = Array1::from_vec(vec![2.0, 4.0, 6.0, 8.0, 10.0]);
        let result = pearsonr(&x.view(), &y.view()).expect("Pearson correlation should succeed");
        // Perfect correlation should be 1.0
        assert!((result - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_minimize() {
        let initial = Array1::from_vec(vec![0.0]);
        let objective = |x: &Array1<f64>| (x[0] - 2.0).powi(2);

        let result =
            minimize(objective, &initial, "nelder-mead").expect("Minimization should succeed");

        // Should find minimum near x = 2.0
        assert!((result.x[0] - 2.0).abs() < 0.5);
        assert!(result.success);
    }

    #[test]
    fn test_inv_identity_property() {
        use scirs2_core::ndarray::array;
        let a = array![[4.0, 7.0], [2.0, 6.0]];
        let a_inv = inv(&a.view()).expect("inverse should succeed");
        for r in 0..2 {
            for c in 0..2 {
                let mut acc = 0.0;
                for k in 0..2 {
                    acc += a[(r, k)] * a_inv[(k, c)];
                }
                let expected = if r == c { 1.0 } else { 0.0 };
                assert!((acc - expected).abs() < 1e-6, "[{r},{c}]={acc}");
            }
        }
    }

    #[test]
    fn test_inv_3x3_identity_property() {
        use scirs2_core::ndarray::array;
        let a = array![[2.0, 1.0, 1.0], [1.0, 3.0, 2.0], [1.0, 0.0, 0.0]];
        let a_inv = inv(&a.view()).expect("inverse should succeed");
        for r in 0..3 {
            for c in 0..3 {
                let mut acc = 0.0;
                for k in 0..3 {
                    acc += a[(r, k)] * a_inv[(k, c)];
                }
                let expected = if r == c { 1.0 } else { 0.0 };
                assert!((acc - expected).abs() < 1e-6, "[{r},{c}]={acc}");
            }
        }
    }

    #[test]
    fn test_inv_singular_errors() {
        use scirs2_core::ndarray::array;
        let a = array![[1.0, 2.0], [2.0, 4.0]];
        assert!(inv(&a.view()).is_err());
    }

    #[test]
    fn test_inv_non_square_errors() {
        use scirs2_core::ndarray::array;
        let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]];
        assert!(inv(&a.view()).is_err());
    }
}