quantum/algorithms.py

"""
FireEcho Quantum Gold - Standard Quantum Algorithms

Implements common quantum algorithms and state preparations:
  - Bell states (maximally entangled 2-qubit states)
  - GHZ states (n-qubit entangled states)
  - Quantum Fourier Transform (QFT)
  - Grover's search algorithm primitives

These serve as both utilities and benchmarks for the simulator.
"""

import torch
import math
from typing import Optional, List

from .circuit import QuantumCircuit
from .simulator import QuantumSimulator, StateVector


def bell_state(variant: int = 0, device: str = 'cuda:0') -> StateVector:
    """
    Create one of the four Bell states.
    
    Bell states are maximally entangled 2-qubit states:
      - Φ⁺ = (|00⟩ + |11⟩)/√2  (variant=0)
      - Φ⁻ = (|00⟩ - |11⟩)/√2  (variant=1)
      - Ψ⁺ = (|01⟩ + |10⟩)/√2  (variant=2)
      - Ψ⁻ = (|01⟩ - |10⟩)/√2  (variant=3)
    
    Args:
        variant: Which Bell state (0-3)
        device: CUDA device
    
    Returns:
        Bell state vector
    
    Example:
        state = bell_state(0)  # (|00⟩ + |11⟩)/√2
        
        # Verify entanglement
        from .measurement import entanglement_entropy
        S = entanglement_entropy(state, [0])  # Should be 1.0
    """
    qc = QuantumCircuit(2, f"bell_{variant}")
    
    # Start with |00⟩
    # Apply H to qubit 0: (|0⟩ + |1⟩)/√2 ⊗ |0⟩
    qc.h(0)
    
    # CNOT creates entanglement: (|00⟩ + |11⟩)/√2
    qc.cx(0, 1)
    
    # Variants modify the state
    if variant == 1:
        # Φ⁻: Apply Z to add relative phase
        qc.z(0)
    elif variant == 2:
        # Ψ⁺: Apply X to flip second qubit
        qc.x(1)
    elif variant == 3:
        # Ψ⁻: Apply both
        qc.z(0)
        qc.x(1)
    
    sim = QuantumSimulator(device)
    return sim.run(qc)


def ghz_state(num_qubits: int, device: str = 'cuda:0') -> StateVector:
    """
    Create a GHZ (Greenberger-Horne-Zeilinger) state.
    
    GHZ state: (|00...0⟩ + |11...1⟩)/√2
    
    This is the maximally entangled n-qubit state, generalizing
    the Bell state to n qubits.
    
    Args:
        num_qubits: Number of qubits (≥2)
        device: CUDA device
    
    Returns:
        GHZ state vector
    
    Example:
        state = ghz_state(3)  # (|000⟩ + |111⟩)/√2
        
        # Sample measurements - only "000" or "111"
        counts = sample(state, shots=1000)
    """
    if num_qubits < 2:
        raise ValueError("GHZ state requires at least 2 qubits")
    
    qc = QuantumCircuit(num_qubits, f"ghz_{num_qubits}")
    
    # Hadamard on first qubit
    qc.h(0)
    
    # CNOT cascade
    for i in range(1, num_qubits):
        qc.cx(0, i)
    
    sim = QuantumSimulator(device)
    return sim.run(qc)


def w_state(num_qubits: int, device: str = 'cuda:0') -> StateVector:
    """
    Create a W state.
    
    W state: (|100...0⟩ + |010...0⟩ + ... + |000...1⟩)/√n
    
    W states are entangled but more robust to qubit loss than GHZ.
    
    Args:
        num_qubits: Number of qubits (≥2)
        device: CUDA device
    
    Returns:
        W state vector
    """
    if num_qubits < 2:
        raise ValueError("W state requires at least 2 qubits")
    
    # Direct construction
    state = StateVector.zeros(num_qubits, device)
    
    norm = 1.0 / math.sqrt(num_qubits)
    for i in range(num_qubits):
        idx = 1 << i  # Single 1 in position i
        state.amplitudes[idx] = norm
    
    state.amplitudes[0] = 0  # Clear |000...0⟩
    
    return state


def qft(num_qubits: int) -> QuantumCircuit:
    """
    Create Quantum Fourier Transform circuit.
    
    QFT transforms computational basis states to Fourier basis:
    |j⟩ → (1/√N) Σₖ e^(2πijk/N) |k⟩
    
    QFT is a key subroutine in Shor's algorithm and quantum
    phase estimation.
    
    Args:
        num_qubits: Number of qubits
    
    Returns:
        QFT circuit
    
    Example:
        qc = qft(4)
        sim = QuantumSimulator()
        state = sim.run(qc)
    """
    qc = QuantumCircuit(num_qubits, f"qft_{num_qubits}")
    
    for i in range(num_qubits):
        # Hadamard on qubit i
        qc.h(i)
        
        # Controlled rotations
        for j in range(i + 1, num_qubits):
            angle = math.pi / (2 ** (j - i))
            qc.cp(angle, j, i)
    
    # Swap qubits to reverse order (standard QFT convention)
    for i in range(num_qubits // 2):
        qc.swap(i, num_qubits - 1 - i)
    
    return qc


def inverse_qft(num_qubits: int) -> QuantumCircuit:
    """
    Create inverse Quantum Fourier Transform circuit.
    
    QFT† is the adjoint (inverse) of QFT:
    QFT† · QFT = I
    
    Args:
        num_qubits: Number of qubits
    
    Returns:
        Inverse QFT circuit
    """
    return qft(num_qubits).inverse()


def grover_diffusion(num_qubits: int) -> QuantumCircuit:
    """
    Create Grover diffusion operator circuit.
    
    D = 2|s⟩⟨s| - I where |s⟩ is uniform superposition.
    Also known as the "inversion about the mean" operator.
    
    Args:
        num_qubits: Number of qubits
    
    Returns:
        Diffusion operator circuit
    """
    qc = QuantumCircuit(num_qubits, "grover_diffusion")
    
    # H⊗n
    for i in range(num_qubits):
        qc.h(i)
    
    # X⊗n
    for i in range(num_qubits):
        qc.x(i)
    
    # Multi-controlled Z (via decomposition)
    if num_qubits == 2:
        qc.cz(0, 1)
    elif num_qubits == 3:
        # CCZ = H-CCX-H on target
        qc.h(2)
        qc.ccx(0, 1, 2)
        qc.h(2)
    else:
        # General multi-controlled Z
        # Use H on last qubit, multi-controlled X, H again
        qc.h(num_qubits - 1)
        # Decompose multi-controlled X (simplified)
        for i in range(num_qubits - 2):
            qc.ccx(i, i + 1, num_qubits - 1)
        qc.h(num_qubits - 1)
    
    # X⊗n
    for i in range(num_qubits):
        qc.x(i)
    
    # H⊗n
    for i in range(num_qubits):
        qc.h(i)
    
    return qc


def quantum_phase_estimation(num_counting_qubits: int, unitary_circuit: QuantumCircuit) -> QuantumCircuit:
    """
    Create Quantum Phase Estimation circuit.
    
    QPE estimates the phase φ in U|ψ⟩ = e^(2πiφ)|ψ⟩.
    
    Args:
        num_counting_qubits: Precision qubits for phase estimate
        unitary_circuit: Circuit implementing unitary U
    
    Returns:
        QPE circuit
    
    Note: The unitary eigenstate should be prepared separately.
    """
    total_qubits = num_counting_qubits + unitary_circuit.num_qubits
    qc = QuantumCircuit(total_qubits, "qpe")
    
    # Hadamard on counting qubits
    for i in range(num_counting_qubits):
        qc.h(i)
    
    # Controlled-U^(2^k) operations
    for k in range(num_counting_qubits):
        # Apply U^(2^k) controlled by qubit k
        repetitions = 2 ** k
        for _ in range(repetitions):
            # Add controlled version of unitary
            # (simplified - actual implementation needs controlled gates)
            for gate in unitary_circuit.gates:
                if gate.name == "RZ":
                    qc.crz(gate.params[0], k, num_counting_qubits + gate.targets[0])
    
    # Inverse QFT on counting qubits
    inv_qft = inverse_qft(num_counting_qubits)
    qc.compose(inv_qft, list(range(num_counting_qubits)))
    
    return qc


def random_circuit(num_qubits: int, depth: int, seed: Optional[int] = None) -> QuantumCircuit:
    """
    Create a random quantum circuit.
    
    Useful for benchmarking and testing.
    
    Args:
        num_qubits: Number of qubits
        depth: Circuit depth
        seed: Random seed
    
    Returns:
        Random circuit
    """
    import random
    if seed is not None:
        random.seed(seed)
    
    qc = QuantumCircuit(num_qubits, f"random_{num_qubits}x{depth}")
    
    single_gates = ['h', 'x', 'y', 'z', 's', 't']
    rotation_gates = ['rx', 'ry', 'rz']
    
    for _ in range(depth):
        # Single-qubit layer
        for q in range(num_qubits):
            gate_type = random.choice(single_gates + rotation_gates)
            
            if gate_type in single_gates:
                getattr(qc, gate_type)(q)
            else:
                angle = random.uniform(0, 2 * math.pi)
                getattr(qc, gate_type)(angle, q)
        
        # Two-qubit layer (CNOTs on adjacent pairs)
        for q in range(0, num_qubits - 1, 2):
            if random.random() > 0.5:
                qc.cx(q, q + 1)
    
    return qc


# =============================================================================
# Variational Circuits (for VQE/QAOA)
# =============================================================================

def variational_ansatz(num_qubits: int, num_layers: int, params: List[float]) -> QuantumCircuit:
    """
    Create a variational ansatz circuit for VQE.
    
    Hardware-efficient ansatz with Ry-CNOT structure.
    
    Args:
        num_qubits: Number of qubits
        num_layers: Number of variational layers
        params: Rotation parameters (length = num_qubits * num_layers * 2)
    
    Returns:
        Parameterized circuit
    """
    expected_params = num_qubits * num_layers * 2
    if len(params) != expected_params:
        raise ValueError(f"Expected {expected_params} parameters, got {len(params)}")
    
    qc = QuantumCircuit(num_qubits, f"vqe_ansatz_{num_layers}L")
    
    param_idx = 0
    for layer in range(num_layers):
        # Rotation layer
        for q in range(num_qubits):
            qc.ry(params[param_idx], q)
            param_idx += 1
            qc.rz(params[param_idx], q)
            param_idx += 1
        
        # Entangling layer (linear connectivity)
        for q in range(num_qubits - 1):
            qc.cx(q, q + 1)
    
    return qc


def qaoa_circuit(num_qubits: int, p: int, gamma: List[float], beta: List[float]) -> QuantumCircuit:
    """
    Create QAOA (Quantum Approximate Optimization Algorithm) circuit.
    
    Standard QAOA ansatz for combinatorial optimization.
    
    Args:
        num_qubits: Number of qubits
        p: Number of QAOA layers
        gamma: Cost unitary parameters
        beta: Mixer unitary parameters
    
    Returns:
        QAOA circuit
    """
    if len(gamma) != p or len(beta) != p:
        raise ValueError(f"Expected {p} gamma and beta values each")
    
    qc = QuantumCircuit(num_qubits, f"qaoa_p{p}")
    
    # Initial superposition
    for q in range(num_qubits):
        qc.h(q)
    
    for layer in range(p):
        # Cost unitary (problem-dependent ZZ interactions)
        for i in range(num_qubits - 1):
            qc.cx(i, i + 1)
            qc.rz(gamma[layer], i + 1)
            qc.cx(i, i + 1)
        
        # Mixer unitary (X rotations)
        for q in range(num_qubits):
            qc.rx(2 * beta[layer], q)
    
    return qc


# =============================================================================
# VQE (Variational Quantum Eigensolver)
# =============================================================================

class VQE:
    """
    Variational Quantum Eigensolver for finding ground state energies.
    
    VQE is a hybrid quantum-classical algorithm that uses:
    1. Quantum circuit to prepare trial wavefunctions
    2. Classical optimizer to minimize energy expectation value
    
    Example:
        from quantum.algorithms import VQE
        
        # Define Hamiltonian (e.g., H2 molecule)
        hamiltonian = [
            (0.5, 'ZZ', [0, 1]),
            (-0.5, 'X', [0]),
            (-0.5, 'X', [1]),
        ]
        
        vqe = VQE(num_qubits=2, num_layers=2)
        energy, params = vqe.run(hamiltonian)
    """
    
    def __init__(
        self,
        num_qubits: int,
        num_layers: int = 2,
        ansatz_type: str = 'hardware_efficient',
        device: str = 'cuda:0'
    ):
        """
        Args:
            num_qubits: Number of qubits
            num_layers: Depth of variational ansatz
            ansatz_type: 'hardware_efficient', 'uccsd', or 'hea'
            device: CUDA device
        """
        self.num_qubits = num_qubits
        self.num_layers = num_layers
        self.ansatz_type = ansatz_type
        self.device = device
        self.sim = QuantumSimulator(device)
        
        # Parameter count depends on ansatz
        if ansatz_type == 'hardware_efficient':
            self.num_params = num_qubits * num_layers * 2
        elif ansatz_type == 'uccsd':
            self.num_params = num_qubits * (num_qubits - 1)
        else:
            self.num_params = num_qubits * num_layers * 3
    
    def build_circuit(self, params: List[float]) -> QuantumCircuit:
        """Build variational circuit with given parameters."""
        if self.ansatz_type == 'hardware_efficient':
            return variational_ansatz(self.num_qubits, self.num_layers, params)
        elif self.ansatz_type == 'uccsd':
            return self._build_uccsd_ansatz(params)
        else:
            return self._build_hea_ansatz(params)
    
    def _build_uccsd_ansatz(self, params: List[float]) -> QuantumCircuit:
        """Build UCCSD (Unitary Coupled Cluster) ansatz."""
        qc = QuantumCircuit(self.num_qubits, "uccsd")
        
        # Hartree-Fock initial state (alternating |1⟩|0⟩)
        for i in range(0, self.num_qubits, 2):
            qc.x(i)
        
        # Single excitations
        param_idx = 0
        for p in range(0, self.num_qubits, 2):
            for q in range(1, self.num_qubits, 2):
                if param_idx < len(params):
                    # Givens rotation for single excitation
                    theta = params[param_idx]
                    qc.cx(p, q)
                    qc.ry(theta, p)
                    qc.cx(p, q)
                    param_idx += 1
        
        return qc
    
    def _build_hea_ansatz(self, params: List[float]) -> QuantumCircuit:
        """Build Hardware-Efficient Ansatz with Ry-Rz-CNOT."""
        qc = QuantumCircuit(self.num_qubits, "hea")
        
        param_idx = 0
        for layer in range(self.num_layers):
            # Ry layer
            for q in range(self.num_qubits):
                qc.ry(params[param_idx], q)
                param_idx += 1
            
            # Rz layer
            for q in range(self.num_qubits):
                qc.rz(params[param_idx], q)
                param_idx += 1
            
            # Rx layer
            for q in range(self.num_qubits):
                qc.rx(params[param_idx], q)
                param_idx += 1
            
            # Entangling layer
            for q in range(self.num_qubits - 1):
                qc.cx(q, q + 1)
            if self.num_qubits > 2:
                qc.cx(self.num_qubits - 1, 0)  # Ring topology
        
        return qc
    
    def compute_expectation(
        self,
        params: List[float],
        hamiltonian: List[tuple]
    ) -> float:
        """
        Compute expectation value ⟨ψ|H|ψ⟩.
        
        Args:
            params: Variational parameters
            hamiltonian: List of (coeff, pauli_string, qubits) tuples
                        e.g., [(0.5, 'ZZ', [0,1]), (-0.3, 'X', [0])]
        
        Returns:
            Energy expectation value
        """
        qc = self.build_circuit(params)
        state = self.sim.run(qc)
        
        total_energy = 0.0
        
        for coeff, pauli_string, qubits in hamiltonian:
            # Measure in appropriate basis
            exp_val = self._measure_pauli_string(state, pauli_string, qubits)
            total_energy += coeff * exp_val
        
        return total_energy
    
    def _measure_pauli_string(
        self,
        state: StateVector,
        pauli_string: str,
        qubits: List[int]
    ) -> float:
        """
        Measure expectation of Pauli string on specified qubits.
        
        For multi-qubit states, we need to properly handle the tensor structure.
        """
        import torch
        
        if len(pauli_string) != len(qubits):
            raise ValueError("Pauli string length must match qubit count")
        
        # Define Pauli matrices
        I = torch.eye(2, dtype=torch.complex64, device=state.amplitudes.device)
        X = torch.tensor([[0, 1], [1, 0]], dtype=torch.complex64, device=state.amplitudes.device)
        Y = torch.tensor([[0, -1j], [1j, 0]], dtype=torch.complex64, device=state.amplitudes.device)
        Z = torch.tensor([[1, 0], [0, -1]], dtype=torch.complex64, device=state.amplitudes.device)
        
        paulis = {'I': I, 'X': X, 'Y': Y, 'Z': Z}
        
        # Build full observable by tensor product
        num_qubits = state.num_qubits
        
        # Start with identity on all qubits
        ops = [I.clone() for _ in range(num_qubits)]
        
        # Place Paulis on specified qubits
        for pauli, qubit in zip(pauli_string, qubits):
            ops[qubit] = paulis[pauli]
        
        # Compute tensor product
        full_obs = ops[0]
        for op in ops[1:]:
            full_obs = torch.kron(full_obs, op)
        
        # Compute expectation: ⟨ψ|O|ψ⟩
        psi = state.amplitudes
        o_psi = torch.mv(full_obs, psi)
        expectation = torch.vdot(psi, o_psi).real
        
        return expectation.item()
    
    def run(
        self,
        hamiltonian: List[tuple],
        max_iters: int = 100,
        learning_rate: float = 0.1,
        callback: Optional[callable] = None
    ) -> tuple:
        """
        Run VQE optimization.
        
        Args:
            hamiltonian: List of (coeff, pauli_string, qubits)
            max_iters: Maximum optimization iterations
            learning_rate: Gradient descent step size
            callback: Optional callback(iter, energy, params)
        
        Returns:
            (final_energy, optimal_params)
        """
        import random
        
        # Initialize random parameters
        params = [random.uniform(-math.pi, math.pi) for _ in range(self.num_params)]
        
        best_energy = float('inf')
        best_params = params.copy()
        
        for iteration in range(max_iters):
            # Compute energy
            energy = self.compute_expectation(params, hamiltonian)
            
            if energy < best_energy:
                best_energy = energy
                best_params = params.copy()
            
            if callback:
                callback(iteration, energy, params)
            
            # Parameter shift gradient
            gradients = []
            for i in range(len(params)):
                params_plus = params.copy()
                params_minus = params.copy()
                params_plus[i] += math.pi / 2
                params_minus[i] -= math.pi / 2
                
                e_plus = self.compute_expectation(params_plus, hamiltonian)
                e_minus = self.compute_expectation(params_minus, hamiltonian)
                
                grad = (e_plus - e_minus) / 2
                gradients.append(grad)
            
            # Update parameters
            for i in range(len(params)):
                params[i] -= learning_rate * gradients[i]
        
        return best_energy, best_params


# =============================================================================
# QSVM (Quantum Support Vector Machine)
# =============================================================================

class QSVM:
    """
    Quantum Support Vector Machine for classification.
    
    Uses quantum feature maps to encode classical data into quantum states,
    enabling kernel-based classification in exponentially large Hilbert space.
    
    Example:
        from quantum.algorithms import QSVM
        
        # Binary classification
        qsvm = QSVM(num_features=4, num_qubits=4)
        
        # Train
        X_train = [[0.1, 0.2, 0.3, 0.4], ...]
        y_train = [0, 1, 0, 1, ...]
        qsvm.fit(X_train, y_train)
        
        # Predict
        predictions = qsvm.predict(X_test)
    """
    
    def __init__(
        self,
        num_features: int,
        num_qubits: Optional[int] = None,
        feature_map: str = 'zz',
        num_layers: int = 2,
        device: str = 'cuda:0'
    ):
        """
        Args:
            num_features: Dimension of input data
            num_qubits: Number of qubits (defaults to num_features)
            feature_map: 'zz', 'pauli', or 'iqp'
            num_layers: Feature map depth
            device: CUDA device
        """
        self.num_features = num_features
        self.num_qubits = num_qubits or num_features
        self.feature_map_type = feature_map
        self.num_layers = num_layers
        self.device = device
        self.sim = QuantumSimulator(device)
        
        # Storage for training data (for kernel computation)
        self.X_train = None
        self.y_train = None
        self.alpha = None  # SVM dual coefficients
    
    def _build_feature_map(self, x: List[float]) -> QuantumCircuit:
        """Build quantum feature map circuit for data point x."""
        if len(x) < self.num_qubits:
            x = list(x) + [0.0] * (self.num_qubits - len(x))
        
        qc = QuantumCircuit(self.num_qubits, "feature_map")
        
        if self.feature_map_type == 'zz':
            return self._zz_feature_map(qc, x)
        elif self.feature_map_type == 'pauli':
            return self._pauli_feature_map(qc, x)
        else:
            return self._iqp_feature_map(qc, x)
    
    def _zz_feature_map(self, qc: QuantumCircuit, x: List[float]) -> QuantumCircuit:
        """ZZ Feature Map - encodes data through ZZ interactions."""
        for layer in range(self.num_layers):
            # Hadamard layer
            for q in range(self.num_qubits):
                qc.h(q)
            
            # Feature encoding layer
            for q in range(self.num_qubits):
                qc.rz(2 * x[q], q)
            
            # Entangling layer with product features
            for i in range(self.num_qubits - 1):
                qc.cx(i, i + 1)
                qc.rz(2 * (math.pi - x[i]) * (math.pi - x[i + 1]), i + 1)
                qc.cx(i, i + 1)
        
        return qc
    
    def _pauli_feature_map(self, qc: QuantumCircuit, x: List[float]) -> QuantumCircuit:
        """Pauli Feature Map - encodes through Pauli rotations."""
        for layer in range(self.num_layers):
            # Hadamard layer
            for q in range(self.num_qubits):
                qc.h(q)
            
            # Z rotations
            for q in range(self.num_qubits):
                qc.rz(x[q], q)
            
            # ZZ interactions
            for i in range(self.num_qubits - 1):
                qc.cx(i, i + 1)
                qc.rz(x[i] * x[i + 1], i + 1)
                qc.cx(i, i + 1)
        
        return qc
    
    def _iqp_feature_map(self, qc: QuantumCircuit, x: List[float]) -> QuantumCircuit:
        """IQP (Instantaneous Quantum Polynomial) Feature Map."""
        for layer in range(self.num_layers):
            # Hadamard layer
            for q in range(self.num_qubits):
                qc.h(q)
            
            # Diagonal gates encoding
            for q in range(self.num_qubits):
                qc.rz(x[q], q)
            
            # Cross terms
            for i in range(self.num_qubits):
                for j in range(i + 1, self.num_qubits):
                    qc.cx(i, j)
                    qc.rz(x[i] * x[j], j)
                    qc.cx(i, j)
        
        return qc
    
    def compute_kernel(self, x1: List[float], x2: List[float]) -> float:
        """
        Compute quantum kernel K(x1, x2) = |⟨φ(x1)|φ(x2)⟩|².
        
        This is the probability of measuring |0...0⟩ after preparing
        the state U†(x2)U(x1)|0⟩.
        """
        # Build circuits
        qc1 = self._build_feature_map(x1)
        qc2 = self._build_feature_map(x2)
        
        # Combined circuit: U(x1) followed by U†(x2)
        combined = QuantumCircuit(self.num_qubits, "kernel")
        combined.compose(qc1, list(range(self.num_qubits)))
        combined.compose(qc2.inverse(), list(range(self.num_qubits)))
        
        # Run and get probability of |0...0⟩
        state = self.sim.run(combined)
        p_zero = (state.amplitudes[0].abs() ** 2).item()
        
        return p_zero
    
    def compute_kernel_matrix(self, X: List[List[float]]) -> torch.Tensor:
        """Compute full kernel matrix for dataset."""
        n = len(X)
        K = torch.zeros(n, n, device=self.device)
        
        for i in range(n):
            for j in range(i, n):
                k_ij = self.compute_kernel(X[i], X[j])
                K[i, j] = k_ij
                K[j, i] = k_ij
        
        return K
    
    def fit(self, X: List[List[float]], y: List[int], C: float = 1.0):
        """
        Fit QSVM to training data.
        
        Args:
            X: Training features, shape [n_samples, n_features]
            y: Training labels, {0, 1} or {-1, 1}
            C: Regularization parameter
        """
        self.X_train = X
        self.y_train = [1 if label > 0 else -1 for label in y]
        n = len(X)
        
        # Compute kernel matrix
        K = self.compute_kernel_matrix(X)
        
        # Convert to numpy for SVM solver
        K_np = K.cpu().numpy()
        y_np = torch.tensor(self.y_train, dtype=torch.float32).numpy()
        
        # Simple gradient descent for dual SVM
        self.alpha = torch.zeros(n, device=self.device)
        
        for iteration in range(100):
            for i in range(n):
                # Compute gradient for alpha[i]
                grad = 1.0
                for j in range(n):
                    grad -= self.alpha[j].item() * self.y_train[j] * self.y_train[i] * K[i, j].item()
                
                # Update
                self.alpha[i] = max(0, min(C, self.alpha[i] + 0.01 * grad))
    
    def predict(self, X: List[List[float]]) -> List[int]:
        """Predict labels for new data."""
        if self.X_train is None:
            raise ValueError("Model not fitted. Call fit() first.")
        
        predictions = []
        
        for x in X:
            # Compute kernel with all training points
            decision = 0.0
            for i, x_train in enumerate(self.X_train):
                k = self.compute_kernel(x, x_train)
                decision += self.alpha[i].item() * self.y_train[i] * k
            
            predictions.append(1 if decision > 0 else 0)
        
        return predictions
    
    def score(self, X: List[List[float]], y: List[int]) -> float:
        """Compute classification accuracy."""
        predictions = self.predict(X)
        correct = sum(1 for p, t in zip(predictions, y) if p == t)
        return correct / len(y)


# =============================================================================
# Quantum Autoencoder
# =============================================================================

def quantum_autoencoder_circuit(
    num_qubits: int,
    latent_qubits: int,
    params: List[float]
) -> QuantumCircuit:
    """
    Create a quantum autoencoder circuit.
    
    Compresses num_qubits down to latent_qubits through a trash-latent separation.
    
    Args:
        num_qubits: Input dimension
        latent_qubits: Compressed dimension
        params: Variational parameters
    
    Returns:
        Autoencoder circuit
    """
    if latent_qubits >= num_qubits:
        raise ValueError("latent_qubits must be < num_qubits")
    
    trash_qubits = num_qubits - latent_qubits
    
    qc = QuantumCircuit(num_qubits, f"qae_{num_qubits}to{latent_qubits}")
    
    # Encoder - variational layers
    param_idx = 0
    num_layers = min(2, len(params) // num_qubits)
    
    for layer in range(num_layers):
        for q in range(num_qubits):
            if param_idx < len(params):
                qc.ry(params[param_idx], q)
                param_idx += 1
        
        for q in range(num_qubits - 1):
            qc.cx(q, q + 1)
    
    # Compression: SWAP latent qubits to the beginning
    # This moves qubits [0..latent-1] to the front
    for i in range(trash_qubits):
        for j in range(trash_qubits - i):
            if j < num_qubits - 1:
                qc.swap(j, j + 1)
    
    return qc


# =============================================================================
# Quantum Principal Component Analysis
# =============================================================================

def quantum_pca_circuit(num_qubits: int, num_components: int) -> QuantumCircuit:
    """
    Create a quantum PCA circuit using quantum phase estimation.
    
    Args:
        num_qubits: Data dimension
        num_components: Number of principal components to extract
    
    Returns:
        QPCA circuit
    """
    total_qubits = num_qubits + num_components
    qc = QuantumCircuit(total_qubits, f"qpca_{num_components}")
    
    # Prepare superposition on counting qubits
    for i in range(num_components):
        qc.h(i)
    
    # Controlled rotations encoding covariance structure
    for k in range(num_components):
        for q in range(num_qubits):
            angle = math.pi / (2 ** (k + 1))
            qc.crz(angle, k, num_components + q)
    
    # Inverse QFT on counting qubits
    for i in range(num_components // 2):
        qc.swap(i, num_components - 1 - i)
    
    for i in range(num_components):
        qc.h(i)
        for j in range(i + 1, num_components):
            angle = -math.pi / (2 ** (j - i))
            qc.cp(angle, j, i)
    
    return qc