quantrs2-device 0.1.3

//! Statistical analysis components

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

use super::super::core::SciRS2ProcessTomographer;
use super::super::results::*;
use crate::DeviceResult;

// Conditional imports
#[cfg(feature = "scirs2")]
use scirs2_stats::{mean, std, var};

#[cfg(feature = "scirs2")]
use scirs2_linalg::eigvals;

#[cfg(not(feature = "scirs2"))]
use super::super::fallback::{mean, std, var};

impl SciRS2ProcessTomographer {
    /// Fit statistical distribution to data
    pub(crate) fn fit_distribution(
        &self,
        data: &[f64],
        name: &str,
    ) -> DeviceResult<ElementDistribution> {
        #[cfg(feature = "scirs2")]
        {
            let data_array = Array1::from_vec(data.to_vec());
            let data_view = data_array.view();

            let mean_val = mean(&data_view).unwrap_or(0.0);
            let std_val = std(&data_view, 0, None).unwrap_or(1.0);

            // Test goodness of fit for normal distribution
            let mut goodness_of_fit = 0.0;
            let mut distribution_type = "normal".to_string();
            let mut parameters = vec![mean_val, std_val];

            // Try fitting different distributions and select best fit
            if data.iter().all(|&x| x >= 0.0) {
                // Try gamma distribution for positive data
                let gamma_alpha = mean_val * mean_val / (std_val * std_val);
                let gamma_beta = mean_val / (std_val * std_val);

                if gamma_alpha > 0.0 && gamma_beta > 0.0 {
                    distribution_type = "gamma".to_string();
                    parameters = vec![gamma_alpha, gamma_beta];
                    goodness_of_fit = 0.85; // Placeholder
                }
            }

            // Calculate confidence interval
            let confidence_interval = (
                mean_val - 1.96 * std_val / (data.len() as f64).sqrt(),
                mean_val + 1.96 * std_val / (data.len() as f64).sqrt(),
            );

            Ok(ElementDistribution {
                distribution_type,
                parameters,
                goodness_of_fit,
                confidence_interval,
            })
        }

        #[cfg(not(feature = "scirs2"))]
        {
            // Fallback implementation
            let mean_val = data.iter().sum::<f64>() / data.len() as f64;
            let var_val =
                data.iter().map(|x| (x - mean_val).powi(2)).sum::<f64>() / data.len() as f64;
            let std_val = var_val.sqrt();

            Ok(ElementDistribution {
                distribution_type: "normal".to_string(),
                parameters: vec![mean_val, std_val],
                goodness_of_fit: 0.9,
                confidence_interval: (mean_val - std_val, mean_val + std_val),
            })
        }
    }

    /// Calculate skewness of data
    pub(crate) fn calculate_skewness(&self, data: &[f64]) -> f64 {
        if data.len() < 3 {
            return 0.0;
        }

        let n = data.len() as f64;
        let mean = data.iter().sum::<f64>() / n;
        let variance = data.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / n;
        let std_dev = variance.sqrt();

        if std_dev < 1e-12 {
            return 0.0;
        }

        let skewness = data
            .iter()
            .map(|&x| ((x - mean) / std_dev).powi(3))
            .sum::<f64>()
            / n;

        skewness
    }

    /// Calculate kurtosis of data
    pub(crate) fn calculate_kurtosis(&self, data: &[f64]) -> f64 {
        if data.len() < 4 {
            return 0.0;
        }

        let n = data.len() as f64;
        let mean = data.iter().sum::<f64>() / n;
        let variance = data.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / n;
        let std_dev = variance.sqrt();

        if std_dev < 1e-12 {
            return 0.0;
        }

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

        kurtosis
    }

    /// Calculate entropy of data
    pub(crate) fn calculate_entropy(&self, data: &[f64]) -> f64 {
        if data.is_empty() {
            return 0.0;
        }

        // Create histogram for discrete entropy calculation
        let num_bins = (data.len() as f64).sqrt() as usize + 1;
        let min_val = data.iter().copied().fold(f64::INFINITY, f64::min);
        let max_val = data.iter().copied().fold(f64::NEG_INFINITY, f64::max);

        if (max_val - min_val).abs() < 1e-12 {
            return 0.0;
        }

        let bin_width = (max_val - min_val) / num_bins as f64;
        let mut histogram = vec![0; num_bins];

        for &value in data {
            let bin_idx = ((value - min_val) / bin_width).floor() as usize;
            let bin_idx = bin_idx.min(num_bins - 1);
            histogram[bin_idx] += 1;
        }

        // Calculate entropy
        let total_count = data.len() as f64;
        let mut entropy = 0.0;

        for &count in &histogram {
            if count > 0 {
                let probability = count as f64 / total_count;
                entropy -= probability * probability.ln();
            }
        }

        entropy
    }

    /// Analyze eigenvalue distribution of the process matrix
    pub(crate) fn analyze_eigenvalue_distribution(
        &self,
        process_matrix: &Array4<Complex64>,
    ) -> DeviceResult<ElementDistribution> {
        // Convert process matrix to Choi representation for eigenvalue analysis
        let choi_matrix = self.convert_to_choi_matrix(process_matrix)?;

        #[cfg(feature = "scirs2")]
        {
            use scirs2_linalg::eig;

            // Convert complex matrix to real parts for eigenvalue calculation
            let real_matrix = choi_matrix.mapv(|x| x.re);

            // Compute eigenvalues using SciRS2
            if let Ok(eigenvalues) = eigvals(&real_matrix.view(), None) {
                let real_eigenvalues: Vec<f64> = eigenvalues.iter().map(|x| x.re).collect();

                return self.fit_distribution(&real_eigenvalues, "eigenvalues");
            }
        }

        // Fallback
        Ok(ElementDistribution {
            distribution_type: "uniform".to_string(),
            parameters: vec![0.0, 1.0],
            goodness_of_fit: 0.8,
            confidence_interval: (0.0, 1.0),
        })
    }

    /// Convert process matrix to Choi representation
    pub(crate) fn convert_to_choi_matrix(
        &self,
        process_matrix: &Array4<Complex64>,
    ) -> DeviceResult<Array2<Complex64>> {
        let dim = process_matrix.dim().0;
        let choi_dim = dim * dim;
        let mut choi_matrix = Array2::zeros((choi_dim, choi_dim));

        // Convert Chi matrix to Choi matrix
        // Choi[i*d+j, k*d+l] = Chi[i,k,j,l]
        for i in 0..dim {
            for j in 0..dim {
                for k in 0..dim {
                    for l in 0..dim {
                        let choi_row = i * dim + j;
                        let choi_col = k * dim + l;
                        if choi_row < choi_dim && choi_col < choi_dim {
                            choi_matrix[[choi_row, choi_col]] = process_matrix[[i, k, j, l]];
                        }
                    }
                }
            }
        }

        Ok(choi_matrix)
    }

    /// Calculate process coefficient for measurement matrix
    pub(crate) fn compute_process_coefficient(
        &self,
        input_state: &Array2<Complex64>,
        measurement: &Array2<Complex64>,
        i: usize,
        j: usize,
        k: usize,
        l: usize,
    ) -> DeviceResult<f64> {
        // Simplified coefficient calculation
        // In practice, this would involve computing Tr(M * E_{ij} * rho * E_{kl}^†)

        let dim = input_state.dim().0;
        if i >= dim || j >= dim || k >= dim || l >= dim {
            return Ok(0.0);
        }

        // Create basis matrices E_{ij} and E_{kl}
        let mut e_ij = Array2::zeros((dim, dim));
        e_ij[[i, j]] = Complex64::new(1.0, 0.0);

        let mut e_kl = Array2::zeros((dim, dim));
        e_kl[[k, l]] = Complex64::new(1.0, 0.0);

        // Compute the trace: Tr(M * E_ij * rho * E_kl^†)
        let mut result = Complex64::new(0.0, 0.0);
        for p in 0..dim {
            for q in 0..dim {
                for r in 0..dim {
                    for s in 0..dim {
                        result += measurement[[p, q]]
                            * e_ij[[q, r]]
                            * input_state[[r, s]]
                            * e_kl[[s, p]].conj();
                    }
                }
            }
        }

        Ok(result.re)
    }

    /// Generate state combinations for multi-qubit systems
    pub(crate) fn generate_state_combinations(
        &self,
        single_qubit_states: &[Array2<Complex64>],
        num_qubits: usize,
    ) -> DeviceResult<Vec<Array2<Complex64>>> {
        if num_qubits == 1 {
            return Ok(single_qubit_states.to_vec());
        }

        let mut combinations = Vec::new();
        let num_states = single_qubit_states.len();
        let total_combinations = num_states.pow(num_qubits as u32);

        for combo_idx in 0..total_combinations {
            let mut combination = single_qubit_states[0].clone();
            let mut temp_idx = combo_idx;

            for qubit_idx in 1..num_qubits {
                temp_idx /= num_states;
                let state_idx = temp_idx % num_states;
                combination = self.tensor_product(&combination, &single_qubit_states[state_idx])?;
            }

            combinations.push(combination);
        }

        Ok(combinations)
    }

    /// Compute tensor product of two matrices
    pub(crate) fn tensor_product(
        &self,
        a: &Array2<Complex64>,
        b: &Array2<Complex64>,
    ) -> DeviceResult<Array2<Complex64>> {
        let a_shape = a.dim();
        let b_shape = b.dim();
        let result_shape = (a_shape.0 * b_shape.0, a_shape.1 * b_shape.1);

        let mut result = Array2::zeros(result_shape);

        for i in 0..a_shape.0 {
            for j in 0..a_shape.1 {
                for k in 0..b_shape.0 {
                    for l in 0..b_shape.1 {
                        let result_i = i * b_shape.0 + k;
                        let result_j = j * b_shape.1 + l;
                        result[[result_i, result_j]] = a[[i, j]] * b[[k, l]];
                    }
                }
            }
        }

        Ok(result)
    }

    /// Generate Pauli tensor products
    pub(crate) fn generate_pauli_tensor_products(
        &self,
        single_qubit_paulis: &[Array2<Complex64>],
        num_qubits: usize,
    ) -> DeviceResult<Vec<Array2<Complex64>>> {
        if num_qubits == 1 {
            return Ok(single_qubit_paulis.to_vec());
        }

        let mut tensor_products = Vec::new();
        let num_paulis = single_qubit_paulis.len();
        let total_products = num_paulis.pow(num_qubits as u32);

        for product_idx in 0..total_products {
            let mut product = single_qubit_paulis[0].clone();
            let mut temp_idx = product_idx;

            for qubit_idx in 1..num_qubits {
                temp_idx /= num_paulis;
                let pauli_idx = temp_idx % num_paulis;
                product = self.tensor_product(&product, &single_qubit_paulis[pauli_idx])?;
            }

            tensor_products.push(product);
        }

        Ok(tensor_products)
    }
}