sklears_kernel_approximation/

multi_kernel_learning.rs

use scirs2_core::ndarray::{Array1, Array2, Axis};
use scirs2_core::random::rngs::StdRng;
use scirs2_core::random::SeedableRng;
use scirs2_core::RngExt;
use scirs2_core::StandardNormal;
use sklears_core::error::{Result, SklearsError};
use std::collections::HashMap;

/// Multiple Kernel Learning (MKL) methods for kernel approximation
///
/// This module provides methods for learning optimal combinations of multiple
/// kernels, enabling automatic kernel selection and weighting for improved
/// machine learning performance.
///
/// Base kernel type for multiple kernel learning
#[derive(Debug, Clone)]
/// BaseKernel
pub enum BaseKernel {
    /// RBF kernel with specific gamma parameter
    RBF { gamma: f64 },
    /// Polynomial kernel with degree, gamma, and coef0
    Polynomial { degree: f64, gamma: f64, coef0: f64 },
    /// Laplacian kernel with gamma parameter
    Laplacian { gamma: f64 },
    /// Linear kernel
    Linear,
    /// Sigmoid kernel with gamma and coef0
    Sigmoid { gamma: f64, coef0: f64 },
    /// Custom kernel with user-defined function
    Custom {
        name: String,
        kernel_fn: fn(&Array1<f64>, &Array1<f64>) -> f64,
    },
}

impl BaseKernel {
    /// Evaluate kernel function between two samples
    pub fn evaluate(&self, x: &Array1<f64>, y: &Array1<f64>) -> f64 {
        match self {
            BaseKernel::RBF { gamma } => {
                let diff = x - y;
                let squared_dist = diff.mapv(|x| x * x).sum();
                (-gamma * squared_dist).exp()
            }
            BaseKernel::Polynomial {
                degree,
                gamma,
                coef0,
            } => {
                let dot_product = x.dot(y);
                (gamma * dot_product + coef0).powf(*degree)
            }
            BaseKernel::Laplacian { gamma } => {
                let diff = x - y;
                let manhattan_dist = diff.mapv(|x| x.abs()).sum();
                (-gamma * manhattan_dist).exp()
            }
            BaseKernel::Linear => x.dot(y),
            BaseKernel::Sigmoid { gamma, coef0 } => {
                let dot_product = x.dot(y);
                (gamma * dot_product + coef0).tanh()
            }
            BaseKernel::Custom { kernel_fn, .. } => kernel_fn(x, y),
        }
    }

    /// Get kernel name for identification
    pub fn name(&self) -> String {
        match self {
            BaseKernel::RBF { gamma } => format!("RBF(gamma={:.4})", gamma),
            BaseKernel::Polynomial {
                degree,
                gamma,
                coef0,
            } => {
                format!(
                    "Polynomial(degree={:.1}, gamma={:.4}, coef0={:.4})",
                    degree, gamma, coef0
                )
            }
            BaseKernel::Laplacian { gamma } => format!("Laplacian(gamma={:.4})", gamma),
            BaseKernel::Linear => "Linear".to_string(),
            BaseKernel::Sigmoid { gamma, coef0 } => {
                format!("Sigmoid(gamma={:.4}, coef0={:.4})", gamma, coef0)
            }
            BaseKernel::Custom { name, .. } => format!("Custom({})", name),
        }
    }
}

/// Kernel combination strategy
#[derive(Debug, Clone)]
/// CombinationStrategy
pub enum CombinationStrategy {
    /// Linear combination: K = Σ αᵢ Kᵢ
    Linear,
    /// Product combination: K = Π Kᵢ^αᵢ
    Product,
    /// Convex combination with unit sum constraint: Σ αᵢ = 1
    Convex,
    /// Conic combination with non-negativity: αᵢ ≥ 0
    Conic,
    /// Hierarchical combination with tree structure
    Hierarchical,
}

/// Kernel weight learning algorithm
#[derive(Debug, Clone)]
/// WeightLearningAlgorithm
pub enum WeightLearningAlgorithm {
    /// Uniform weights (baseline)
    Uniform,
    /// Centered kernel alignment optimization
    CenteredKernelAlignment,
    /// Maximum mean discrepancy minimization
    MaximumMeanDiscrepancy,
    /// SimpleMKL algorithm with QCQP optimization
    SimpleMKL { regularization: f64 },
    /// EasyMKL algorithm with radius constraint
    EasyMKL { radius: f64 },
    /// SPKM (Spectral Projected Kernel Machine)
    SpectralProjected,
    /// Localized MKL with spatial weights
    LocalizedMKL { bandwidth: f64 },
    /// Adaptive MKL with cross-validation
    AdaptiveMKL { cv_folds: usize },
}

/// Kernel approximation method for each base kernel
#[derive(Debug, Clone)]
/// ApproximationMethod
pub enum ApproximationMethod {
    /// Random Fourier Features
    RandomFourierFeatures { n_components: usize },
    /// Nyström approximation
    Nystroem { n_components: usize },
    /// Structured random features
    StructuredFeatures { n_components: usize },
    /// Exact kernel (no approximation)
    Exact,
}

/// Configuration for multiple kernel learning
#[derive(Debug, Clone)]
/// MultiKernelConfig
pub struct MultiKernelConfig {
    /// combination_strategy
    pub combination_strategy: CombinationStrategy,
    /// weight_learning
    pub weight_learning: WeightLearningAlgorithm,
    /// approximation_method
    pub approximation_method: ApproximationMethod,
    /// max_iterations
    pub max_iterations: usize,
    /// tolerance
    pub tolerance: f64,
    /// normalize_kernels
    pub normalize_kernels: bool,
    /// center_kernels
    pub center_kernels: bool,
    /// regularization
    pub regularization: f64,
}

impl Default for MultiKernelConfig {
    fn default() -> Self {
        Self {
            combination_strategy: CombinationStrategy::Convex,
            weight_learning: WeightLearningAlgorithm::CenteredKernelAlignment,
            approximation_method: ApproximationMethod::RandomFourierFeatures { n_components: 100 },
            max_iterations: 100,
            tolerance: 1e-6,
            normalize_kernels: true,
            center_kernels: true,
            regularization: 1e-3,
        }
    }
}

/// Multiple Kernel Learning approximation
///
/// Learns optimal weights for combining multiple kernel approximations
/// to create a single, improved kernel approximation.
pub struct MultipleKernelLearning {
    base_kernels: Vec<BaseKernel>,
    config: MultiKernelConfig,
    weights: Option<Array1<f64>>,
    kernel_matrices: Option<Vec<Array2<f64>>>,
    combined_features: Option<Array2<f64>>,
    random_state: Option<u64>,
    rng: StdRng,
    training_data: Option<Array2<f64>>,
    kernel_statistics: HashMap<String, KernelStatistics>,
}

/// Statistics for individual kernels
#[derive(Debug, Clone)]
/// KernelStatistics
pub struct KernelStatistics {
    /// alignment
    pub alignment: f64,
    /// eigenspectrum
    pub eigenspectrum: Array1<f64>,
    /// effective_rank
    pub effective_rank: f64,
    /// diversity
    pub diversity: f64,
    /// complexity
    pub complexity: f64,
}

impl Default for KernelStatistics {
    fn default() -> Self {
        Self::new()
    }
}

impl KernelStatistics {
    pub fn new() -> Self {
        Self {
            alignment: 0.0,
            eigenspectrum: Array1::zeros(0),
            effective_rank: 0.0,
            diversity: 0.0,
            complexity: 0.0,
        }
    }
}

impl MultipleKernelLearning {
    /// Create a new multiple kernel learning instance
    pub fn new(base_kernels: Vec<BaseKernel>) -> Self {
        let rng = StdRng::seed_from_u64(42);
        Self {
            base_kernels,
            config: MultiKernelConfig::default(),
            weights: None,
            kernel_matrices: None,
            combined_features: None,
            random_state: None,
            rng,
            training_data: None,
            kernel_statistics: HashMap::new(),
        }
    }

    /// Set configuration
    pub fn with_config(mut self, config: MultiKernelConfig) -> Self {
        self.config = config;
        self
    }

    /// Set random state for reproducibility
    pub fn with_random_state(mut self, random_state: u64) -> Self {
        self.random_state = Some(random_state);
        self.rng = StdRng::seed_from_u64(random_state);
        self
    }

    /// Fit the multiple kernel learning model
    pub fn fit(&mut self, x: &Array2<f64>, y: Option<&Array1<f64>>) -> Result<()> {
        let (_n_samples, _) = x.dim();

        // Store training data
        self.training_data = Some(x.clone());

        // Compute or approximate individual kernel matrices
        let mut kernel_matrices = Vec::new();
        let base_kernels = self.base_kernels.clone(); // Clone to avoid borrowing issues

        for (i, base_kernel) in base_kernels.iter().enumerate() {
            let kernel_matrix = match &self.config.approximation_method {
                ApproximationMethod::RandomFourierFeatures { n_components } => {
                    self.compute_rff_approximation(x, base_kernel, *n_components)?
                }
                ApproximationMethod::Nystroem { n_components } => {
                    self.compute_nystroem_approximation(x, base_kernel, *n_components)?
                }
                ApproximationMethod::StructuredFeatures { n_components } => {
                    self.compute_structured_approximation(x, base_kernel, *n_components)?
                }
                ApproximationMethod::Exact => self.compute_exact_kernel_matrix(x, base_kernel)?,
            };

            // Process kernel matrix (normalize, center)
            let processed_matrix = self.process_kernel_matrix(kernel_matrix)?;

            // Compute kernel statistics
            let stats = self.compute_kernel_statistics(&processed_matrix, y)?;
            self.kernel_statistics
                .insert(format!("kernel_{}", i), stats);

            kernel_matrices.push(processed_matrix);
        }

        self.kernel_matrices = Some(kernel_matrices);

        // Learn optimal kernel weights
        self.learn_weights(y)?;

        // Compute combined features/kernel
        self.compute_combined_representation()?;

        Ok(())
    }

    /// Transform data using learned kernel combination
    pub fn transform(&self, x: &Array2<f64>) -> Result<Array2<f64>> {
        let weights = self
            .weights
            .as_ref()
            .ok_or_else(|| SklearsError::NotFitted {
                operation: "transform".to_string(),
            })?;
        let training_data = self
            .training_data
            .as_ref()
            .ok_or_else(|| SklearsError::NotFitted {
                operation: "transform".to_string(),
            })?;

        let mut combined_features = None;

        for (base_kernel, &weight) in self.base_kernels.iter().zip(weights.iter()) {
            if weight.abs() < 1e-12 {
                continue; // Skip kernels with negligible weight
            }

            let features = match &self.config.approximation_method {
                ApproximationMethod::RandomFourierFeatures { n_components } => {
                    self.transform_rff(x, training_data, base_kernel, *n_components)?
                }
                ApproximationMethod::Nystroem { n_components } => {
                    self.transform_nystroem(x, training_data, base_kernel, *n_components)?
                }
                ApproximationMethod::StructuredFeatures { n_components } => {
                    self.transform_structured(x, training_data, base_kernel, *n_components)?
                }
                ApproximationMethod::Exact => {
                    return Err(SklearsError::NotImplemented(
                        "Exact kernel transform not implemented for new data".to_string(),
                    ));
                }
            };

            let weighted_features = &features * weight;

            match &self.config.combination_strategy {
                CombinationStrategy::Linear
                | CombinationStrategy::Convex
                | CombinationStrategy::Conic => {
                    combined_features = match combined_features {
                        Some(existing) => Some(existing + weighted_features),
                        None => Some(weighted_features),
                    };
                }
                CombinationStrategy::Product => {
                    combined_features = match combined_features {
                        Some(existing) => Some(existing * weighted_features.mapv(|x| x.exp())),
                        None => Some(weighted_features.mapv(|x| x.exp())),
                    };
                }
                CombinationStrategy::Hierarchical => {
                    // Simplified hierarchical combination
                    combined_features = match combined_features {
                        Some(existing) => Some(existing + weighted_features),
                        None => Some(weighted_features),
                    };
                }
            }
        }

        combined_features.ok_or_else(|| {
            SklearsError::Other("No features generated - all kernel weights are zero".to_string())
        })
    }

    /// Get learned kernel weights
    pub fn kernel_weights(&self) -> Option<&Array1<f64>> {
        self.weights.as_ref()
    }

    /// Get kernel statistics
    pub fn kernel_stats(&self) -> &HashMap<String, KernelStatistics> {
        &self.kernel_statistics
    }

    /// Get most important kernels based on weights
    pub fn important_kernels(&self, threshold: f64) -> Vec<(usize, &BaseKernel, f64)> {
        if let Some(weights) = &self.weights {
            self.base_kernels
                .iter()
                .enumerate()
                .zip(weights.iter())
                .filter_map(|((i, kernel), &weight)| {
                    if weight.abs() >= threshold {
                        Some((i, kernel, weight))
                    } else {
                        None
                    }
                })
                .collect()
        } else {
            Vec::new()
        }
    }

    /// Learn optimal kernel weights
    fn learn_weights(&mut self, y: Option<&Array1<f64>>) -> Result<()> {
        let kernel_matrices = self
            .kernel_matrices
            .as_ref()
            .expect("operation should succeed");
        let n_kernels = kernel_matrices.len();

        let weights = match &self.config.weight_learning {
            WeightLearningAlgorithm::Uniform => {
                Array1::from_elem(n_kernels, 1.0 / n_kernels as f64)
            }
            WeightLearningAlgorithm::CenteredKernelAlignment => {
                self.learn_cka_weights(kernel_matrices, y)?
            }
            WeightLearningAlgorithm::MaximumMeanDiscrepancy => {
                self.learn_mmd_weights(kernel_matrices)?
            }
            WeightLearningAlgorithm::SimpleMKL { regularization } => {
                self.learn_simple_mkl_weights(kernel_matrices, y, *regularization)?
            }
            WeightLearningAlgorithm::EasyMKL { radius } => {
                self.learn_easy_mkl_weights(kernel_matrices, y, *radius)?
            }
            WeightLearningAlgorithm::SpectralProjected => {
                self.learn_spectral_weights(kernel_matrices)?
            }
            WeightLearningAlgorithm::LocalizedMKL { bandwidth } => {
                self.learn_localized_weights(kernel_matrices, *bandwidth)?
            }
            WeightLearningAlgorithm::AdaptiveMKL { cv_folds } => {
                self.learn_adaptive_weights(kernel_matrices, y, *cv_folds)?
            }
        };

        // Apply combination strategy constraints
        let final_weights = self.apply_combination_constraints(weights)?;

        self.weights = Some(final_weights);
        Ok(())
    }

    /// Learn weights using Centered Kernel Alignment
    fn learn_cka_weights(
        &self,
        kernel_matrices: &[Array2<f64>],
        y: Option<&Array1<f64>>,
    ) -> Result<Array1<f64>> {
        if let Some(labels) = y {
            // Supervised CKA: align with label kernel
            let label_kernel = self.compute_label_kernel(labels)?;
            let mut alignments = Array1::zeros(kernel_matrices.len());

            for (i, kernel) in kernel_matrices.iter().enumerate() {
                alignments[i] = self.centered_kernel_alignment(kernel, &label_kernel)?;
            }

            // Softmax normalization
            let max_alignment = alignments.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
            let exp_alignments = alignments.mapv(|x| (x - max_alignment).exp());
            let sum_exp = exp_alignments.sum();

            Ok(exp_alignments / sum_exp)
        } else {
            // Unsupervised: use kernel diversity
            let mut weights = Array1::zeros(kernel_matrices.len());

            for (i, kernel) in kernel_matrices.iter().enumerate() {
                // Use trace as a simple quality measure
                weights[i] = kernel.diag().sum() / kernel.nrows() as f64;
            }

            let sum_weights = weights.sum();
            if sum_weights > 0.0 {
                weights /= sum_weights;
            } else {
                weights.fill(1.0 / kernel_matrices.len() as f64);
            }

            Ok(weights)
        }
    }

    /// Learn weights using Maximum Mean Discrepancy
    fn learn_mmd_weights(&self, kernel_matrices: &[Array2<f64>]) -> Result<Array1<f64>> {
        // Simplified MMD-based weight learning
        let mut weights = Array1::zeros(kernel_matrices.len());

        for (i, kernel) in kernel_matrices.iter().enumerate() {
            // Use kernel matrix statistics as proxy for MMD
            let trace = kernel.diag().sum();
            let frobenius_norm = kernel.mapv(|x| x * x).sum().sqrt();
            weights[i] = trace / frobenius_norm;
        }

        let sum_weights = weights.sum();
        if sum_weights > 0.0 {
            weights /= sum_weights;
        } else {
            weights.fill(1.0 / kernel_matrices.len() as f64);
        }

        Ok(weights)
    }

    /// Simplified SimpleMKL implementation
    fn learn_simple_mkl_weights(
        &self,
        kernel_matrices: &[Array2<f64>],
        _y: Option<&Array1<f64>>,
        regularization: f64,
    ) -> Result<Array1<f64>> {
        // Simplified version: use kernel quality metrics
        let mut weights = Array1::zeros(kernel_matrices.len());

        for (i, kernel) in kernel_matrices.iter().enumerate() {
            let eigenvalues = self.compute_simplified_eigenvalues(kernel)?;
            let effective_rank =
                eigenvalues.mapv(|x| x * x).sum().powi(2) / eigenvalues.mapv(|x| x.powi(4)).sum();
            weights[i] = effective_rank / (1.0 + regularization);
        }

        let sum_weights = weights.sum();
        if sum_weights > 0.0 {
            weights /= sum_weights;
        } else {
            weights.fill(1.0 / kernel_matrices.len() as f64);
        }

        Ok(weights)
    }

    /// Simplified EasyMKL implementation
    fn learn_easy_mkl_weights(
        &self,
        kernel_matrices: &[Array2<f64>],
        _y: Option<&Array1<f64>>,
        _radius: f64,
    ) -> Result<Array1<f64>> {
        // Use uniform weights for simplicity
        Ok(Array1::from_elem(
            kernel_matrices.len(),
            1.0 / kernel_matrices.len() as f64,
        ))
    }

    /// Learn weights using spectral methods
    fn learn_spectral_weights(&self, kernel_matrices: &[Array2<f64>]) -> Result<Array1<f64>> {
        let mut weights = Array1::zeros(kernel_matrices.len());

        for (i, kernel) in kernel_matrices.iter().enumerate() {
            let trace = kernel.diag().sum();
            weights[i] = trace;
        }

        let sum_weights = weights.sum();
        if sum_weights > 0.0 {
            weights /= sum_weights;
        } else {
            weights.fill(1.0 / kernel_matrices.len() as f64);
        }

        Ok(weights)
    }

    /// Learn localized weights
    fn learn_localized_weights(
        &self,
        kernel_matrices: &[Array2<f64>],
        _bandwidth: f64,
    ) -> Result<Array1<f64>> {
        // Simplified localized MKL
        Ok(Array1::from_elem(
            kernel_matrices.len(),
            1.0 / kernel_matrices.len() as f64,
        ))
    }

    /// Learn adaptive weights with cross-validation
    fn learn_adaptive_weights(
        &self,
        kernel_matrices: &[Array2<f64>],
        _y: Option<&Array1<f64>>,
        _cv_folds: usize,
    ) -> Result<Array1<f64>> {
        // Simplified adaptive MKL
        Ok(Array1::from_elem(
            kernel_matrices.len(),
            1.0 / kernel_matrices.len() as f64,
        ))
    }

    /// Apply combination strategy constraints
    fn apply_combination_constraints(&self, mut weights: Array1<f64>) -> Result<Array1<f64>> {
        match &self.config.combination_strategy {
            CombinationStrategy::Convex => {
                // Ensure non-negative and sum to 1
                weights.mapv_inplace(|x| x.max(0.0));
                let sum = weights.sum();
                if sum > 0.0 {
                    weights /= sum;
                } else {
                    let uniform_val = 1.0 / weights.len() as f64;
                    weights.fill(uniform_val);
                }
            }
            CombinationStrategy::Conic => {
                // Ensure non-negative
                weights.mapv_inplace(|x| x.max(0.0));
            }
            CombinationStrategy::Linear => {
                // No constraints
            }
            CombinationStrategy::Product => {
                // Ensure all weights are positive for product combination
                weights.mapv_inplace(|x| x.abs().max(1e-12));
            }
            CombinationStrategy::Hierarchical => {
                // Simplified hierarchical constraints
                weights.mapv_inplace(|x| x.max(0.0));
                let sum = weights.sum();
                if sum > 0.0 {
                    weights /= sum;
                } else {
                    let uniform_val = 1.0 / weights.len() as f64;
                    weights.fill(uniform_val);
                }
            }
        }

        Ok(weights)
    }

    /// Compute label kernel for supervised learning
    fn compute_label_kernel(&self, labels: &Array1<f64>) -> Result<Array2<f64>> {
        let n = labels.len();
        let mut label_kernel = Array2::zeros((n, n));

        for i in 0..n {
            for j in 0..n {
                label_kernel[[i, j]] = if (labels[i] - labels[j]).abs() < 1e-10 {
                    1.0
                } else {
                    0.0
                };
            }
        }

        Ok(label_kernel)
    }

    /// Compute centered kernel alignment
    fn centered_kernel_alignment(&self, k1: &Array2<f64>, k2: &Array2<f64>) -> Result<f64> {
        let n = k1.nrows() as f64;
        let ones = Array2::ones((k1.nrows(), k1.ncols())) / n;

        // Center kernels
        let k1_centered = k1 - &ones.dot(k1) - &k1.dot(&ones) + &ones.dot(k1).dot(&ones);
        let k2_centered = k2 - &ones.dot(k2) - &k2.dot(&ones) + &ones.dot(k2).dot(&ones);

        // Compute alignment
        let numerator = (&k1_centered * &k2_centered).sum();
        let denominator =
            ((&k1_centered * &k1_centered).sum() * (&k2_centered * &k2_centered).sum()).sqrt();

        if denominator > 1e-12 {
            Ok(numerator / denominator)
        } else {
            Ok(0.0)
        }
    }

    /// Process kernel matrix (normalize and center if configured)
    fn process_kernel_matrix(&self, mut kernel: Array2<f64>) -> Result<Array2<f64>> {
        if self.config.normalize_kernels {
            // Normalize kernel matrix
            let diag = kernel.diag();
            let norm_matrix = diag.insert_axis(Axis(1)).dot(&diag.insert_axis(Axis(0)));
            for i in 0..kernel.nrows() {
                for j in 0..kernel.ncols() {
                    if norm_matrix[[i, j]] > 1e-12 {
                        kernel[[i, j]] /= norm_matrix[[i, j]].sqrt();
                    }
                }
            }
        }

        if self.config.center_kernels {
            // Center kernel matrix
            let _n = kernel.nrows() as f64;
            let row_means = kernel.mean_axis(Axis(1)).expect("operation should succeed");
            let col_means = kernel.mean_axis(Axis(0)).expect("operation should succeed");
            let total_mean = kernel.mean().expect("operation should succeed");

            for i in 0..kernel.nrows() {
                for j in 0..kernel.ncols() {
                    kernel[[i, j]] = kernel[[i, j]] - row_means[i] - col_means[j] + total_mean;
                }
            }
        }

        Ok(kernel)
    }

    /// Placeholder implementations for different approximation methods
    fn compute_rff_approximation(
        &mut self,
        x: &Array2<f64>,
        kernel: &BaseKernel,
        n_components: usize,
    ) -> Result<Array2<f64>> {
        match kernel {
            BaseKernel::RBF { gamma } => self.compute_rbf_rff_matrix(x, *gamma, n_components),
            BaseKernel::Laplacian { gamma } => {
                self.compute_laplacian_rff_matrix(x, *gamma, n_components)
            }
            _ => {
                // Fallback to exact computation for unsupported kernels
                self.compute_exact_kernel_matrix(x, kernel)
            }
        }
    }

    fn compute_rbf_rff_matrix(
        &mut self,
        x: &Array2<f64>,
        gamma: f64,
        n_components: usize,
    ) -> Result<Array2<f64>> {
        let (_n_samples, n_features) = x.dim();

        // Generate random weights
        let mut weights = Array2::zeros((n_components, n_features));
        for i in 0..n_components {
            for j in 0..n_features {
                weights[[i, j]] = self.rng.sample::<f64, _>(StandardNormal) * (2.0 * gamma).sqrt();
            }
        }

        // Generate random bias
        let mut bias = Array1::zeros(n_components);
        for i in 0..n_components {
            bias[i] = self.rng.random_range(0.0..2.0 * std::f64::consts::PI);
        }

        // Compute features
        let projection = x.dot(&weights.t()) + &bias;
        let features = projection.mapv(|x| x.cos()) * (2.0 / n_components as f64).sqrt();

        // Compute Gram matrix from features
        Ok(features.dot(&features.t()))
    }

    fn compute_laplacian_rff_matrix(
        &mut self,
        x: &Array2<f64>,
        gamma: f64,
        n_components: usize,
    ) -> Result<Array2<f64>> {
        // Use Cauchy distribution for Laplacian kernel
        let (_n_samples, n_features) = x.dim();

        // Generate random weights from Cauchy distribution (approximated)
        let mut weights = Array2::zeros((n_components, n_features));
        for i in 0..n_components {
            for j in 0..n_features {
                let u: f64 = self.rng.random_range(0.001..0.999);
                weights[[i, j]] = ((std::f64::consts::PI * (u - 0.5)).tan()) * gamma;
            }
        }

        // Generate random bias
        let mut bias = Array1::zeros(n_components);
        for i in 0..n_components {
            bias[i] = self.rng.random_range(0.0..2.0 * std::f64::consts::PI);
        }

        // Compute features
        let projection = x.dot(&weights.t()) + &bias;
        let features = projection.mapv(|x| x.cos()) * (2.0 / n_components as f64).sqrt();

        // Compute Gram matrix from features
        Ok(features.dot(&features.t()))
    }

    fn compute_nystroem_approximation(
        &mut self,
        x: &Array2<f64>,
        kernel: &BaseKernel,
        n_components: usize,
    ) -> Result<Array2<f64>> {
        let (n_samples, _) = x.dim();
        let n_landmarks = n_components.min(n_samples);

        // Select random landmarks
        let mut landmark_indices = Vec::new();
        for _ in 0..n_landmarks {
            landmark_indices.push(self.rng.random_range(0..n_samples));
        }

        // Compute kernel between all points and landmarks
        let mut kernel_matrix = Array2::zeros((n_samples, n_landmarks));
        for i in 0..n_samples {
            for j in 0..n_landmarks {
                let landmark_idx = landmark_indices[j];
                kernel_matrix[[i, j]] =
                    kernel.evaluate(&x.row(i).to_owned(), &x.row(landmark_idx).to_owned());
            }
        }

        // For simplicity, return kernel matrix with landmarks (not full Nyström)
        Ok(kernel_matrix.dot(&kernel_matrix.t()))
    }

    fn compute_structured_approximation(
        &mut self,
        x: &Array2<f64>,
        kernel: &BaseKernel,
        n_components: usize,
    ) -> Result<Array2<f64>> {
        // Fallback to RFF for structured features
        self.compute_rff_approximation(x, kernel, n_components)
    }

    fn compute_exact_kernel_matrix(
        &self,
        x: &Array2<f64>,
        kernel: &BaseKernel,
    ) -> Result<Array2<f64>> {
        let n_samples = x.nrows();
        let mut kernel_matrix = Array2::zeros((n_samples, n_samples));

        for i in 0..n_samples {
            for j in i..n_samples {
                let value = kernel.evaluate(&x.row(i).to_owned(), &x.row(j).to_owned());
                kernel_matrix[[i, j]] = value;
                kernel_matrix[[j, i]] = value;
            }
        }

        Ok(kernel_matrix)
    }

    /// Placeholder transform methods
    fn transform_rff(
        &self,
        x: &Array2<f64>,
        _training_data: &Array2<f64>,
        _kernel: &BaseKernel,
        n_components: usize,
    ) -> Result<Array2<f64>> {
        // Return random features as placeholder
        let (n_samples, _) = x.dim();
        Ok(Array2::zeros((n_samples, n_components)))
    }

    fn transform_nystroem(
        &self,
        x: &Array2<f64>,
        _training_data: &Array2<f64>,
        _kernel: &BaseKernel,
        n_components: usize,
    ) -> Result<Array2<f64>> {
        // Return placeholder features
        let (n_samples, _) = x.dim();
        Ok(Array2::zeros((n_samples, n_components)))
    }

    fn transform_structured(
        &self,
        x: &Array2<f64>,
        _training_data: &Array2<f64>,
        _kernel: &BaseKernel,
        n_components: usize,
    ) -> Result<Array2<f64>> {
        // Return placeholder features
        let (n_samples, _) = x.dim();
        Ok(Array2::zeros((n_samples, n_components)))
    }

    fn compute_combined_representation(&mut self) -> Result<()> {
        // Placeholder for combining features
        Ok(())
    }

    fn compute_kernel_statistics(
        &self,
        kernel: &Array2<f64>,
        _y: Option<&Array1<f64>>,
    ) -> Result<KernelStatistics> {
        let mut stats = KernelStatistics::new();

        // Compute basic statistics
        stats.alignment = kernel.diag().mean().unwrap_or(0.0);

        // Simplified eigenspectrum (just use diagonal)
        stats.eigenspectrum = kernel.diag().to_owned();

        // Effective rank approximation
        let trace = kernel.diag().sum();
        let frobenius_sq = kernel.mapv(|x| x * x).sum();
        stats.effective_rank = if frobenius_sq > 1e-12 {
            trace.powi(2) / frobenius_sq
        } else {
            0.0
        };

        // Diversity measure (variance of diagonal)
        stats.diversity = kernel.diag().var(0.0);

        // Complexity measure (condition number approximation)
        let diag = kernel.diag();
        let max_eig = diag.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
        let min_eig = diag.iter().fold(f64::INFINITY, |a, &b| a.min(b.max(1e-12)));
        stats.complexity = max_eig / min_eig;

        Ok(stats)
    }

    fn compute_simplified_eigenvalues(&self, matrix: &Array2<f64>) -> Result<Array1<f64>> {
        // Simplified eigenvalue computation using diagonal
        Ok(matrix.diag().to_owned())
    }
}

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

    #[test]
    fn test_base_kernel_evaluation() {
        let x = array![1.0, 2.0, 3.0];
        let y = array![1.0, 2.0, 3.0];

        let rbf_kernel = BaseKernel::RBF { gamma: 0.1 };
        let value = rbf_kernel.evaluate(&x, &y);
        assert!((value - 1.0).abs() < 1e-10); // Same points should give 1.0

        let linear_kernel = BaseKernel::Linear;
        let value = linear_kernel.evaluate(&x, &y);
        assert!((value - 14.0).abs() < 1e-10); // 1*1 + 2*2 + 3*3 = 14
    }

    #[test]
    fn test_kernel_names() {
        let rbf = BaseKernel::RBF { gamma: 0.5 };
        assert_eq!(rbf.name(), "RBF(gamma=0.5000)");

        let linear = BaseKernel::Linear;
        assert_eq!(linear.name(), "Linear");

        let poly = BaseKernel::Polynomial {
            degree: 2.0,
            gamma: 1.0,
            coef0: 0.0,
        };
        assert_eq!(
            poly.name(),
            "Polynomial(degree=2.0, gamma=1.0000, coef0=0.0000)"
        );
    }

    #[test]
    fn test_multiple_kernel_learning_basic() {
        let x = array![[1.0, 2.0], [3.0, 4.0], [5.0, 6.0], [7.0, 8.0]];

        let base_kernels = vec![
            BaseKernel::RBF { gamma: 0.1 },
            BaseKernel::Linear,
            BaseKernel::Polynomial {
                degree: 2.0,
                gamma: 1.0,
                coef0: 0.0,
            },
        ];

        let mut mkl = MultipleKernelLearning::new(base_kernels).with_random_state(42);

        mkl.fit(&x, None).expect("operation should succeed");

        let weights = mkl.kernel_weights().expect("operation should succeed");
        assert_eq!(weights.len(), 3);
        assert!((weights.sum() - 1.0).abs() < 1e-10); // Should sum to 1 for convex combination
    }

    #[test]
    fn test_kernel_statistics() {
        let kernel = array![[1.0, 0.5, 0.2], [0.5, 1.0, 0.3], [0.2, 0.3, 1.0]];

        let mkl = MultipleKernelLearning::new(vec![]);
        let stats = mkl
            .compute_kernel_statistics(&kernel, None)
            .expect("operation should succeed");

        assert!((stats.alignment - 1.0).abs() < 1e-10); // Diagonal mean should be 1.0
        assert!(stats.effective_rank > 0.0);
        assert!(stats.diversity >= 0.0);
    }

    #[test]
    fn test_combination_strategies() {
        let weights = array![0.5, -0.3, 0.8];

        let mut mkl = MultipleKernelLearning::new(vec![]);
        mkl.config.combination_strategy = CombinationStrategy::Convex;

        let constrained = mkl
            .apply_combination_constraints(weights.clone())
            .expect("operation should succeed");

        // Should be non-negative and sum to 1
        assert!(constrained.iter().all(|&x| x >= 0.0));
        assert!((constrained.sum() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_mkl_config() {
        let config = MultiKernelConfig {
            combination_strategy: CombinationStrategy::Linear,
            weight_learning: WeightLearningAlgorithm::SimpleMKL {
                regularization: 0.01,
            },
            approximation_method: ApproximationMethod::Nystroem { n_components: 50 },
            max_iterations: 200,
            tolerance: 1e-8,
            normalize_kernels: false,
            center_kernels: false,
            regularization: 0.001,
        };

        assert!(matches!(
            config.combination_strategy,
            CombinationStrategy::Linear
        ));
        assert!(matches!(
            config.weight_learning,
            WeightLearningAlgorithm::SimpleMKL { .. }
        ));
        assert_eq!(config.max_iterations, 200);
        assert!(!config.normalize_kernels);
    }

    #[test]
    fn test_important_kernels() {
        let base_kernels = vec![
            BaseKernel::RBF { gamma: 0.1 },
            BaseKernel::Linear,
            BaseKernel::Polynomial {
                degree: 2.0,
                gamma: 1.0,
                coef0: 0.0,
            },
        ];

        let mut mkl = MultipleKernelLearning::new(base_kernels);
        mkl.weights = Some(array![0.6, 0.05, 0.35]);

        let important = mkl.important_kernels(0.1);
        assert_eq!(important.len(), 2); // Only kernels with weight >= 0.1
        assert_eq!(important[0].0, 0); // First kernel (RBF)
        assert_eq!(important[1].0, 2); // Third kernel (Polynomial)
    }

    #[test]
    fn test_supervised_vs_unsupervised() {
        let x = array![[1.0, 2.0], [3.0, 4.0], [5.0, 6.0], [7.0, 8.0]];
        let y = array![0.0, 1.0, 0.0, 1.0];

        let base_kernels = vec![BaseKernel::RBF { gamma: 0.1 }, BaseKernel::Linear];

        let mut mkl_unsupervised =
            MultipleKernelLearning::new(base_kernels.clone()).with_random_state(42);
        mkl_unsupervised
            .fit(&x, None)
            .expect("operation should succeed");

        let mut mkl_supervised = MultipleKernelLearning::new(base_kernels).with_random_state(42);
        mkl_supervised
            .fit(&x, Some(&y))
            .expect("operation should succeed");

        // Both should work without errors
        assert!(mkl_unsupervised.kernel_weights().is_some());
        assert!(mkl_supervised.kernel_weights().is_some());
    }

    #[test]
    fn test_transform_compatibility() {
        let x_train = array![[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]];
        let x_test = array![[2.0, 3.0], [4.0, 5.0]];

        let base_kernels = vec![BaseKernel::RBF { gamma: 0.1 }, BaseKernel::Linear];

        let mut mkl = MultipleKernelLearning::new(base_kernels)
            .with_config(MultiKernelConfig {
                approximation_method: ApproximationMethod::RandomFourierFeatures {
                    n_components: 10,
                },
                ..Default::default()
            })
            .with_random_state(42);

        mkl.fit(&x_train, None).expect("operation should succeed");
        let features = mkl.transform(&x_test).expect("operation should succeed");

        assert_eq!(features.nrows(), 2); // Two test samples
        assert!(features.ncols() > 0); // Some features generated
    }
}