quantrs2_core/qml/

advanced_algorithms.rs

//! Advanced Quantum Machine Learning Algorithms
//!
//! This module provides sophisticated QML algorithms including:
//! - Quantum Kernel Methods for SVM and kernel-based classifiers
//! - Quantum Transfer Learning for pre-trained circuit reuse
//! - Quantum Ensemble Methods for combining multiple quantum models
//! - Quantum Feature Maps with advanced embedding strategies

use crate::error::{QuantRS2Error, QuantRS2Result};
use crate::qml::{EncodingStrategy, EntanglementPattern, QMLConfig, QMLLayer};
use scirs2_core::ndarray::{Array1, Array2};
use scirs2_core::Complex64;
use std::collections::HashMap;
use std::sync::Arc;

// =============================================================================
// Quantum Kernel Methods
// =============================================================================

/// Quantum kernel configuration
#[derive(Debug, Clone)]
pub struct QuantumKernelConfig {
    /// Number of qubits
    pub num_qubits: usize,
    /// Feature map type
    pub feature_map: FeatureMapType,
    /// Number of repetitions
    pub reps: usize,
    /// Entanglement pattern
    pub entanglement: EntanglementPattern,
    /// Parameter scaling
    pub parameter_scaling: f64,
}

impl Default for QuantumKernelConfig {
    fn default() -> Self {
        Self {
            num_qubits: 4,
            feature_map: FeatureMapType::ZZFeatureMap,
            reps: 2,
            entanglement: EntanglementPattern::Full,
            parameter_scaling: 2.0,
        }
    }
}

/// Feature map types for quantum kernels
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum FeatureMapType {
    /// ZZ feature map with entanglement
    ZZFeatureMap,
    /// Pauli feature map
    PauliFeatureMap,
    /// IQP feature map
    IQPFeatureMap,
    /// Custom trainable feature map
    TrainableFeatureMap,
}

/// Quantum kernel for kernel-based machine learning
pub struct QuantumKernel {
    /// Configuration
    config: QuantumKernelConfig,
    /// Cached kernel matrix
    kernel_cache: Option<Array2<f64>>,
    /// Training data (for caching)
    training_data: Option<Array2<f64>>,
}

impl QuantumKernel {
    /// Create a new quantum kernel
    pub fn new(config: QuantumKernelConfig) -> Self {
        Self {
            config,
            kernel_cache: None,
            training_data: None,
        }
    }

    /// Compute kernel value between two data points
    pub fn kernel(&self, x1: &[f64], x2: &[f64]) -> QuantRS2Result<f64> {
        if x1.len() != self.config.num_qubits || x2.len() != self.config.num_qubits {
            return Err(QuantRS2Error::InvalidInput(format!(
                "Data dimension {} doesn't match num_qubits {}",
                x1.len(),
                self.config.num_qubits
            )));
        }

        // Compute kernel as |<φ(x1)|φ(x2)>|²
        // This is a simplified implementation
        let state1 = self.encode_data(x1)?;
        let state2 = self.encode_data(x2)?;

        // Inner product
        let inner: Complex64 = state1
            .iter()
            .zip(state2.iter())
            .map(|(a, b)| a.conj() * b)
            .sum();

        Ok(inner.norm_sqr())
    }

    /// Compute kernel matrix for dataset
    pub fn kernel_matrix(&mut self, data: &Array2<f64>) -> QuantRS2Result<Array2<f64>> {
        let n_samples = data.nrows();
        let mut kernel_matrix = Array2::zeros((n_samples, n_samples));

        for i in 0..n_samples {
            for j in i..n_samples {
                let x_i = data.row(i).to_vec();
                let x_j = data.row(j).to_vec();

                let k_ij = self.kernel(&x_i, &x_j)?;
                kernel_matrix[[i, j]] = k_ij;
                kernel_matrix[[j, i]] = k_ij; // Symmetric
            }
        }

        // Cache for later use
        self.kernel_cache = Some(kernel_matrix.clone());
        self.training_data = Some(data.clone());

        Ok(kernel_matrix)
    }

    /// Encode data into quantum state using feature map
    fn encode_data(&self, data: &[f64]) -> QuantRS2Result<Array1<Complex64>> {
        let dim = 1 << self.config.num_qubits;
        let mut state = Array1::zeros(dim);
        state[0] = Complex64::new(1.0, 0.0);

        // Apply feature map encoding
        match self.config.feature_map {
            FeatureMapType::ZZFeatureMap => {
                self.apply_zz_feature_map(&mut state, data)?;
            }
            FeatureMapType::PauliFeatureMap => {
                self.apply_pauli_feature_map(&mut state, data)?;
            }
            FeatureMapType::IQPFeatureMap => {
                self.apply_iqp_feature_map(&mut state, data)?;
            }
            FeatureMapType::TrainableFeatureMap => {
                self.apply_trainable_feature_map(&mut state, data)?;
            }
        }

        Ok(state)
    }

    fn apply_zz_feature_map(
        &self,
        state: &mut Array1<Complex64>,
        data: &[f64],
    ) -> QuantRS2Result<()> {
        for _ in 0..self.config.reps {
            // Single-qubit rotations
            for (i, &x) in data.iter().enumerate() {
                let angle = self.config.parameter_scaling * x;
                self.apply_rz(state, i, angle);
                self.apply_ry(state, i, angle);
            }

            // Two-qubit interactions
            for i in 0..self.config.num_qubits - 1 {
                let angle = self.config.parameter_scaling
                    * (std::f64::consts::PI - data[i])
                    * (std::f64::consts::PI - data[i + 1]);
                self.apply_rzz(state, i, i + 1, angle);
            }
        }
        Ok(())
    }

    fn apply_pauli_feature_map(
        &self,
        state: &mut Array1<Complex64>,
        data: &[f64],
    ) -> QuantRS2Result<()> {
        for _ in 0..self.config.reps {
            for (i, &x) in data.iter().enumerate() {
                let angle = self.config.parameter_scaling * x;
                self.apply_rx(state, i, angle);
                self.apply_rz(state, i, angle);
            }
        }
        Ok(())
    }

    fn apply_iqp_feature_map(
        &self,
        state: &mut Array1<Complex64>,
        data: &[f64],
    ) -> QuantRS2Result<()> {
        // Hadamard layer
        for i in 0..self.config.num_qubits {
            self.apply_hadamard(state, i);
        }

        // Diagonal gates
        for (i, &x) in data.iter().enumerate() {
            let angle = self.config.parameter_scaling * x * x;
            self.apply_rz(state, i, angle);
        }

        Ok(())
    }

    fn apply_trainable_feature_map(
        &self,
        state: &mut Array1<Complex64>,
        data: &[f64],
    ) -> QuantRS2Result<()> {
        // Combination of encoding and trainable parameters
        for _ in 0..self.config.reps {
            for (i, &x) in data.iter().enumerate() {
                self.apply_ry(state, i, x);
                self.apply_rz(state, i, x);
            }
        }
        Ok(())
    }

    // Gate implementations
    fn apply_rx(&self, state: &mut Array1<Complex64>, qubit: usize, angle: f64) {
        let cos = (angle / 2.0).cos();
        let sin = (angle / 2.0).sin();
        let dim = state.len();
        let mask = 1 << qubit;

        for i in 0..dim / 2 {
            let idx0 = (i & !(mask >> 1)) | ((i & (mask >> 1)) << 1);
            let idx1 = idx0 | mask;

            if idx1 < dim {
                let a = state[idx0];
                let b = state[idx1];

                state[idx0] = Complex64::new(cos, 0.0) * a + Complex64::new(0.0, -sin) * b;
                state[idx1] = Complex64::new(0.0, -sin) * a + Complex64::new(cos, 0.0) * b;
            }
        }
    }

    fn apply_ry(&self, state: &mut Array1<Complex64>, qubit: usize, angle: f64) {
        let cos = (angle / 2.0).cos();
        let sin = (angle / 2.0).sin();
        let dim = state.len();
        let mask = 1 << qubit;

        for i in 0..dim / 2 {
            let idx0 = (i & !(mask >> 1)) | ((i & (mask >> 1)) << 1);
            let idx1 = idx0 | mask;

            if idx1 < dim {
                let a = state[idx0];
                let b = state[idx1];

                state[idx0] = Complex64::new(cos, 0.0) * a - Complex64::new(sin, 0.0) * b;
                state[idx1] = Complex64::new(sin, 0.0) * a + Complex64::new(cos, 0.0) * b;
            }
        }
    }

    fn apply_rz(&self, state: &mut Array1<Complex64>, qubit: usize, angle: f64) {
        let dim = state.len();
        let mask = 1 << qubit;

        for i in 0..dim {
            if i & mask != 0 {
                state[i] *= Complex64::new(0.0, angle / 2.0).exp();
            } else {
                state[i] *= Complex64::new(0.0, -angle / 2.0).exp();
            }
        }
    }

    fn apply_rzz(&self, state: &mut Array1<Complex64>, q1: usize, q2: usize, angle: f64) {
        let dim = state.len();
        let mask1 = 1 << q1;
        let mask2 = 1 << q2;

        for i in 0..dim {
            let bit1 = (i & mask1) != 0;
            let bit2 = (i & mask2) != 0;
            let parity = if bit1 == bit2 { 1.0 } else { -1.0 };

            state[i] *= Complex64::new(0.0, parity * angle / 2.0).exp();
        }
    }

    fn apply_hadamard(&self, state: &mut Array1<Complex64>, qubit: usize) {
        let inv_sqrt2 = 1.0 / std::f64::consts::SQRT_2;
        let dim = state.len();
        let mask = 1 << qubit;

        for i in 0..dim / 2 {
            let idx0 = (i & !(mask >> 1)) | ((i & (mask >> 1)) << 1);
            let idx1 = idx0 | mask;

            if idx1 < dim {
                let a = state[idx0];
                let b = state[idx1];

                state[idx0] = Complex64::new(inv_sqrt2, 0.0) * (a + b);
                state[idx1] = Complex64::new(inv_sqrt2, 0.0) * (a - b);
            }
        }
    }
}

// =============================================================================
// Quantum Support Vector Machine
// =============================================================================

/// Quantum SVM classifier
pub struct QuantumSVM {
    /// Quantum kernel
    kernel: QuantumKernel,
    /// Support vector indices
    support_vectors: Vec<usize>,
    /// Dual coefficients
    alphas: Vec<f64>,
    /// Bias term
    bias: f64,
    /// Training labels
    labels: Vec<f64>,
    /// Training data
    training_data: Option<Array2<f64>>,
}

impl QuantumSVM {
    /// Create a new Quantum SVM
    pub fn new(kernel_config: QuantumKernelConfig) -> Self {
        Self {
            kernel: QuantumKernel::new(kernel_config),
            support_vectors: Vec::new(),
            alphas: Vec::new(),
            bias: 0.0,
            labels: Vec::new(),
            training_data: None,
        }
    }

    /// Train the QSVM on data
    pub fn fit(&mut self, data: &Array2<f64>, labels: &[f64], c: f64) -> QuantRS2Result<()> {
        let n_samples = data.nrows();

        // Compute kernel matrix
        let kernel_matrix = self.kernel.kernel_matrix(data)?;

        // Simplified SMO-like training
        self.alphas = vec![0.0; n_samples];
        self.labels = labels.to_vec();
        self.training_data = Some(data.clone());

        // Simple gradient descent on dual problem
        let learning_rate = 0.01;
        let max_iter = 100;

        for _ in 0..max_iter {
            for i in 0..n_samples {
                let mut grad = 1.0;
                for j in 0..n_samples {
                    grad -= self.alphas[j] * labels[i] * labels[j] * kernel_matrix[[i, j]];
                }

                self.alphas[i] += learning_rate * grad;
                self.alphas[i] = self.alphas[i].clamp(0.0, c);
            }
        }

        // Find support vectors
        let epsilon = 1e-6;
        self.support_vectors = (0..n_samples)
            .filter(|&i| self.alphas[i] > epsilon)
            .collect();

        // Compute bias
        if !self.support_vectors.is_empty() {
            let sv = self.support_vectors[0];
            let mut b = labels[sv];
            for j in 0..n_samples {
                b -= self.alphas[j] * labels[j] * kernel_matrix[[sv, j]];
            }
            self.bias = b;
        }

        Ok(())
    }

    /// Predict class for new data point
    pub fn predict(&self, x: &[f64]) -> QuantRS2Result<f64> {
        let training_data = self
            .training_data
            .as_ref()
            .ok_or_else(|| QuantRS2Error::RuntimeError("Model not trained".to_string()))?;

        let mut decision = self.bias;

        for &i in &self.support_vectors {
            let x_i = training_data.row(i).to_vec();
            let k = self.kernel.kernel(&x_i, x)?;
            decision += self.alphas[i] * self.labels[i] * k;
        }

        Ok(if decision >= 0.0 { 1.0 } else { -1.0 })
    }

    /// Predict probabilities using Platt scaling approximation
    pub fn predict_proba(&self, x: &[f64]) -> QuantRS2Result<f64> {
        let training_data = self
            .training_data
            .as_ref()
            .ok_or_else(|| QuantRS2Error::RuntimeError("Model not trained".to_string()))?;

        let mut decision = self.bias;
        for &i in &self.support_vectors {
            let x_i = training_data.row(i).to_vec();
            let k = self.kernel.kernel(&x_i, x)?;
            decision += self.alphas[i] * self.labels[i] * k;
        }

        // Sigmoid transformation
        Ok(1.0 / (1.0 + (-decision).exp()))
    }
}

// =============================================================================
// Quantum Transfer Learning
// =============================================================================

/// Transfer learning configuration
#[derive(Debug, Clone)]
pub struct TransferLearningConfig {
    /// Freeze pre-trained layers
    pub freeze_pretrained: bool,
    /// Number of fine-tuning epochs
    pub fine_tune_epochs: usize,
    /// Learning rate for fine-tuning
    pub fine_tune_lr: f64,
    /// Layer to split at (pretrained | new)
    pub split_layer: usize,
}

impl Default for TransferLearningConfig {
    fn default() -> Self {
        Self {
            freeze_pretrained: true,
            fine_tune_epochs: 50,
            fine_tune_lr: 0.01,
            split_layer: 2,
        }
    }
}

/// Quantum transfer learning model
pub struct QuantumTransferLearning {
    /// Pre-trained circuit parameters
    pretrained_params: Vec<f64>,
    /// New trainable parameters
    new_params: Vec<f64>,
    /// Configuration
    config: TransferLearningConfig,
    /// Number of qubits
    num_qubits: usize,
}

impl QuantumTransferLearning {
    /// Create transfer learning model from pre-trained parameters
    pub fn from_pretrained(
        pretrained_params: Vec<f64>,
        num_qubits: usize,
        config: TransferLearningConfig,
    ) -> Self {
        // Initialize new trainable layers
        let new_param_count = num_qubits * 3; // Simple ansatz
        let new_params = vec![0.0; new_param_count];

        Self {
            pretrained_params,
            new_params,
            config,
            num_qubits,
        }
    }

    /// Get all parameters (pretrained + new)
    pub fn parameters(&self) -> Vec<f64> {
        let mut params = self.pretrained_params.clone();
        params.extend(self.new_params.clone());
        params
    }

    /// Get trainable parameters only
    pub fn trainable_parameters(&self) -> &[f64] {
        if self.config.freeze_pretrained {
            &self.new_params
        } else {
            // Would return all, but for simplicity return new only
            &self.new_params
        }
    }

    /// Update trainable parameters
    pub fn update_parameters(&mut self, new_values: &[f64]) -> QuantRS2Result<()> {
        if new_values.len() != self.new_params.len() {
            return Err(QuantRS2Error::InvalidInput(format!(
                "Expected {} parameters, got {}",
                self.new_params.len(),
                new_values.len()
            )));
        }

        self.new_params.copy_from_slice(new_values);
        Ok(())
    }

    /// Get number of trainable parameters
    pub fn num_trainable(&self) -> usize {
        if self.config.freeze_pretrained {
            self.new_params.len()
        } else {
            self.pretrained_params.len() + self.new_params.len()
        }
    }
}

// =============================================================================
// Quantum Ensemble Methods
// =============================================================================

/// Ensemble voting strategy
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum VotingStrategy {
    /// Hard voting (majority vote)
    Hard,
    /// Soft voting (probability averaging)
    Soft,
    /// Weighted voting
    Weighted,
}

/// Quantum ensemble classifier
pub struct QuantumEnsemble {
    /// Individual models (parameters)
    models: Vec<Vec<f64>>,
    /// Model weights
    weights: Vec<f64>,
    /// Voting strategy
    voting: VotingStrategy,
    /// Number of qubits per model
    num_qubits: usize,
}

impl QuantumEnsemble {
    /// Create a new ensemble
    pub fn new(num_qubits: usize, voting: VotingStrategy) -> Self {
        Self {
            models: Vec::new(),
            weights: Vec::new(),
            voting,
            num_qubits,
        }
    }

    /// Add a model to the ensemble
    pub fn add_model(&mut self, params: Vec<f64>, weight: f64) {
        self.models.push(params);
        self.weights.push(weight);
    }

    /// Get number of models in ensemble
    pub fn num_models(&self) -> usize {
        self.models.len()
    }

    /// Combine predictions using voting strategy
    pub fn combine_predictions(&self, predictions: &[f64]) -> QuantRS2Result<f64> {
        if predictions.len() != self.models.len() {
            return Err(QuantRS2Error::InvalidInput(
                "Predictions count doesn't match models".to_string(),
            ));
        }

        match self.voting {
            VotingStrategy::Hard => {
                // Majority vote
                let sum: f64 = predictions.iter().sum();
                Ok(if sum > 0.5 * predictions.len() as f64 {
                    1.0
                } else {
                    0.0
                })
            }
            VotingStrategy::Soft => {
                // Average probabilities
                let avg = predictions.iter().sum::<f64>() / predictions.len() as f64;
                Ok(avg)
            }
            VotingStrategy::Weighted => {
                // Weighted average
                let total_weight: f64 = self.weights.iter().sum();
                let weighted_sum: f64 = predictions
                    .iter()
                    .zip(self.weights.iter())
                    .map(|(p, w)| p * w)
                    .sum();
                Ok(weighted_sum / total_weight)
            }
        }
    }

    /// Bootstrap aggregating (bagging) for ensemble diversity
    pub fn bagging_sample(data: &Array2<f64>, sample_size: usize, seed: u64) -> Array2<f64> {
        use scirs2_core::random::prelude::*;
        let mut rng = seeded_rng(seed);

        let n_samples = data.nrows();
        let n_features = data.ncols();

        let mut sampled = Array2::zeros((sample_size, n_features));
        for i in 0..sample_size {
            let idx = rng.gen_range(0..n_samples);
            sampled.row_mut(i).assign(&data.row(idx));
        }

        sampled
    }
}

// =============================================================================
// Metrics and Evaluation
// =============================================================================

/// QML metrics for model evaluation
pub struct QMLMetrics;

impl QMLMetrics {
    /// Compute accuracy
    pub fn accuracy(predictions: &[f64], labels: &[f64]) -> f64 {
        if predictions.len() != labels.len() {
            return 0.0;
        }

        let correct: usize = predictions
            .iter()
            .zip(labels.iter())
            .filter(|(&p, &l)| (p - l).abs() < 0.5)
            .count();

        correct as f64 / predictions.len() as f64
    }

    /// Compute precision
    pub fn precision(predictions: &[f64], labels: &[f64]) -> f64 {
        let (tp, fp, _, _) = Self::confusion_counts(predictions, labels);
        if tp + fp == 0 {
            0.0
        } else {
            tp as f64 / (tp + fp) as f64
        }
    }

    /// Compute recall
    pub fn recall(predictions: &[f64], labels: &[f64]) -> f64 {
        let (tp, _, _, fn_) = Self::confusion_counts(predictions, labels);
        if tp + fn_ == 0 {
            0.0
        } else {
            tp as f64 / (tp + fn_) as f64
        }
    }

    /// Compute F1 score
    pub fn f1_score(predictions: &[f64], labels: &[f64]) -> f64 {
        let precision = Self::precision(predictions, labels);
        let recall = Self::recall(predictions, labels);

        if precision + recall == 0.0 {
            0.0
        } else {
            2.0 * precision * recall / (precision + recall)
        }
    }

    fn confusion_counts(predictions: &[f64], labels: &[f64]) -> (usize, usize, usize, usize) {
        let mut tp = 0;
        let mut fp = 0;
        let mut tn = 0;
        let mut fn_ = 0;

        for (&p, &l) in predictions.iter().zip(labels.iter()) {
            let pred_pos = p >= 0.5;
            let label_pos = l >= 0.5;

            match (pred_pos, label_pos) {
                (true, true) => tp += 1,
                (true, false) => fp += 1,
                (false, true) => fn_ += 1,
                (false, false) => tn += 1,
            }
        }

        (tp, fp, tn, fn_)
    }
}

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

    #[test]
    fn test_quantum_kernel_config_default() {
        let config = QuantumKernelConfig::default();
        assert_eq!(config.num_qubits, 4);
        assert_eq!(config.reps, 2);
    }

    #[test]
    fn test_quantum_kernel_creation() {
        let config = QuantumKernelConfig {
            num_qubits: 2,
            ..Default::default()
        };
        let kernel = QuantumKernel::new(config);
        assert!(kernel.kernel_cache.is_none());
    }

    #[test]
    fn test_quantum_kernel_value() {
        let config = QuantumKernelConfig {
            num_qubits: 2,
            reps: 1,
            ..Default::default()
        };
        let kernel = QuantumKernel::new(config);

        let x1 = vec![0.5, 0.3];
        let x2 = vec![0.5, 0.3];

        let k = kernel.kernel(&x1, &x2).unwrap();
        assert!(k >= 0.0, "Kernel value should be non-negative");
        // Note: Quantum kernel values are |<φ(x1)|φ(x2)>|², which is bounded by 1 for normalized states
        // but during intermediate computations, the state may not be fully normalized
    }

    #[test]
    fn test_quantum_svm_creation() {
        let config = QuantumKernelConfig {
            num_qubits: 2,
            ..Default::default()
        };
        let qsvm = QuantumSVM::new(config);
        assert!(qsvm.support_vectors.is_empty());
    }

    #[test]
    fn test_transfer_learning_creation() {
        let pretrained = vec![0.1, 0.2, 0.3, 0.4];
        let config = TransferLearningConfig::default();
        let model = QuantumTransferLearning::from_pretrained(pretrained.clone(), 2, config);

        assert_eq!(model.pretrained_params.len(), 4);
        assert!(model.new_params.len() > 0);
    }

    #[test]
    fn test_ensemble_creation() {
        let mut ensemble = QuantumEnsemble::new(2, VotingStrategy::Soft);
        ensemble.add_model(vec![0.1, 0.2], 1.0);
        ensemble.add_model(vec![0.3, 0.4], 1.0);

        assert_eq!(ensemble.num_models(), 2);
    }

    #[test]
    fn test_ensemble_voting() {
        let mut ensemble = QuantumEnsemble::new(2, VotingStrategy::Hard);
        ensemble.add_model(vec![0.1], 1.0);
        ensemble.add_model(vec![0.2], 1.0);
        ensemble.add_model(vec![0.3], 1.0);

        // Majority says 0
        let predictions = vec![0.2, 0.3, 0.4];
        let result = ensemble.combine_predictions(&predictions).unwrap();
        assert_eq!(result, 0.0);

        // Majority says 1
        let predictions = vec![0.6, 0.7, 0.8];
        let result = ensemble.combine_predictions(&predictions).unwrap();
        assert_eq!(result, 1.0);
    }

    #[test]
    fn test_metrics_accuracy() {
        let predictions = vec![1.0, 0.0, 1.0, 1.0];
        let labels = vec![1.0, 0.0, 0.0, 1.0];

        let acc = QMLMetrics::accuracy(&predictions, &labels);
        assert_eq!(acc, 0.75);
    }

    #[test]
    fn test_metrics_precision_recall() {
        let predictions = vec![1.0, 1.0, 0.0, 0.0];
        let labels = vec![1.0, 0.0, 0.0, 1.0];

        let precision = QMLMetrics::precision(&predictions, &labels);
        let recall = QMLMetrics::recall(&predictions, &labels);

        assert_eq!(precision, 0.5); // TP=1, FP=1
        assert_eq!(recall, 0.5); // TP=1, FN=1
    }

    #[test]
    fn test_bagging_sample() {
        let data = Array2::from_shape_vec((10, 3), (0..30).map(|x| x as f64).collect()).unwrap();
        let sample = QuantumEnsemble::bagging_sample(&data, 5, 42);
        assert_eq!(sample.nrows(), 5);
        assert_eq!(sample.ncols(), 3);
    }
}