quantrs2-anneal 0.2.1

//! Infinite-depth QAOA (∞-QAOA) optimizer implementation
//!
//! This module implements the infinite-depth Quantum Approximate Optimization Algorithm
//! with adaptive depth selection, parameter optimization, and convergence detection.

use scirs2_core::random::prelude::*;
use scirs2_core::random::ChaCha8Rng;
use scirs2_core::random::{Rng, SeedableRng};
use scirs2_core::Complex64;
use scirs2_core::RngExt;
use std::collections::HashMap;
use std::f64::consts::PI;
use std::time::{Duration, Instant};

use super::error::{AdvancedQuantumError, AdvancedQuantumResult};
use super::utils::{calculate_relative_improvement, normalize_parameters, validate_parameters};
use crate::ising::{IsingModel, QuboModel};
use crate::simulator::{AnnealingParams, AnnealingResult, AnnealingSolution};

/// Infinite-depth QAOA (∞-QAOA) optimizer
#[derive(Debug, Clone)]
pub struct InfiniteDepthQAOA {
    /// Configuration for ∞-QAOA
    pub config: InfiniteQAOAConfig,
    /// Current parameter history
    pub parameter_history: Vec<Vec<f64>>,
    /// Energy history
    pub energy_history: Vec<f64>,
    /// Depth progression
    pub depth_progression: Vec<usize>,
    /// Convergence metrics
    pub convergence_metrics: ConvergenceMetrics,
    /// Adaptive depth control
    pub depth_controller: AdaptiveDepthController,
    /// Performance statistics
    pub performance_stats: InfiniteQAOAStats,
}

/// Configuration for infinite-depth QAOA
#[derive(Debug, Clone)]
pub struct InfiniteQAOAConfig {
    /// Initial depth
    pub initial_depth: usize,
    /// Maximum depth (for practical limits)
    pub max_depth: usize,
    /// Depth increment strategy
    pub depth_strategy: DepthIncrementStrategy,
    /// Parameter initialization method
    pub initialization_method: ParameterInitializationMethod,
    /// Optimization tolerance
    pub optimization_tolerance: f64,
    /// Maximum optimization iterations per depth
    pub max_iterations_per_depth: usize,
    /// Convergence criteria
    pub convergence_criteria: ConvergenceCriteria,
    /// Classical optimizer configuration
    pub classical_optimizer: ClassicalOptimizerConfig,
    /// Measurement strategy
    pub measurement_strategy: MeasurementStrategy,
    /// Noise mitigation settings
    pub noise_mitigation: NoiseMitigationConfig,
}

impl Default for InfiniteQAOAConfig {
    fn default() -> Self {
        Self {
            initial_depth: 1,
            max_depth: 100,
            depth_strategy: DepthIncrementStrategy::Adaptive,
            initialization_method: ParameterInitializationMethod::Heuristic,
            optimization_tolerance: 1e-6,
            max_iterations_per_depth: 1000,
            convergence_criteria: ConvergenceCriteria::default(),
            classical_optimizer: ClassicalOptimizerConfig::default(),
            measurement_strategy: MeasurementStrategy::default(),
            noise_mitigation: NoiseMitigationConfig::default(),
        }
    }
}

/// Depth increment strategies
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum DepthIncrementStrategy {
    /// Linear increment
    Linear,
    /// Exponential increment
    Exponential,
    /// Adaptive based on convergence
    Adaptive,
    /// Golden ratio increment
    GoldenRatio,
    /// Fibonacci sequence
    Fibonacci,
}

/// Parameter initialization methods
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum ParameterInitializationMethod {
    /// Random initialization
    Random,
    /// Heuristic initialization
    Heuristic,
    /// Transfer from previous depth
    Transfer,
    /// Interpolation-based
    Interpolation,
    /// Machine learning guided
    MLGuided,
}

/// Convergence criteria
#[derive(Debug, Clone)]
pub struct ConvergenceCriteria {
    /// Energy improvement threshold
    pub energy_threshold: f64,
    /// Parameter change threshold
    pub parameter_threshold: f64,
    /// Gradient norm threshold
    pub gradient_threshold: f64,
    /// Maximum stagnation iterations
    pub max_stagnation: usize,
    /// Relative improvement threshold
    pub relative_improvement: f64,
}

impl Default for ConvergenceCriteria {
    fn default() -> Self {
        Self {
            energy_threshold: 1e-8,
            parameter_threshold: 1e-6,
            gradient_threshold: 1e-6,
            max_stagnation: 50,
            relative_improvement: 1e-6,
        }
    }
}

/// Classical optimizer configuration
#[derive(Debug, Clone)]
pub struct ClassicalOptimizerConfig {
    /// Optimizer type
    pub optimizer_type: ClassicalOptimizerType,
    /// Learning rate
    pub learning_rate: f64,
    /// Momentum coefficient
    pub momentum: f64,
    /// L-BFGS memory size
    pub lbfgs_memory: usize,
    /// Maximum function evaluations
    pub max_evaluations: usize,
}

impl Default for ClassicalOptimizerConfig {
    fn default() -> Self {
        Self {
            optimizer_type: ClassicalOptimizerType::LBFGS,
            learning_rate: 0.01,
            momentum: 0.9,
            lbfgs_memory: 10,
            max_evaluations: 1000,
        }
    }
}

/// Classical optimizer types
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum ClassicalOptimizerType {
    /// Gradient Descent
    GradientDescent,
    /// Adam optimizer
    Adam,
    /// L-BFGS
    LBFGS,
    /// Nelder-Mead
    NelderMead,
    /// Powell's method
    Powell,
    /// Differential Evolution
    DifferentialEvolution,
}

/// Measurement strategy
#[derive(Debug, Clone)]
pub struct MeasurementStrategy {
    /// Number of shots per measurement
    pub shots: usize,
    /// Observable decomposition method
    pub observable_decomposition: ObservableDecomposition,
    /// Error mitigation for measurements
    pub error_mitigation: MeasurementErrorMitigation,
    /// Grouping strategy for observables
    pub grouping_strategy: ObservableGrouping,
}

impl Default for MeasurementStrategy {
    fn default() -> Self {
        Self {
            shots: 8192,
            observable_decomposition: ObservableDecomposition::PauliStrings,
            error_mitigation: MeasurementErrorMitigation::ZeroNoiseExtrapolation,
            grouping_strategy: ObservableGrouping::QubitWise,
        }
    }
}

/// Observable decomposition methods
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum ObservableDecomposition {
    /// Pauli string decomposition
    PauliStrings,
    /// Tensor network decomposition
    TensorNetwork,
    /// Clifford decomposition
    Clifford,
    /// Fermionic decomposition
    Fermionic,
}

/// Measurement error mitigation techniques
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum MeasurementErrorMitigation {
    /// No mitigation
    None,
    /// Zero-noise extrapolation
    ZeroNoiseExtrapolation,
    /// Readout error correction
    ReadoutCorrection,
    /// Symmetry verification
    SymmetryVerification,
    /// Virtual distillation
    VirtualDistillation,
}

/// Observable grouping strategies
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum ObservableGrouping {
    /// Group by qubit
    QubitWise,
    /// Group by commutation
    Commuting,
    /// Graph coloring based
    GraphColoring,
    /// Tensor factorization
    TensorFactorization,
}

/// Noise mitigation configuration
#[derive(Debug, Clone)]
pub struct NoiseMitigationConfig {
    /// Enable noise mitigation
    pub enabled: bool,
    /// Mitigation techniques
    pub techniques: Vec<NoiseMitigationTechnique>,
    /// Noise characterization
    pub noise_characterization: NoiseCharacterization,
    /// Error threshold
    pub error_threshold: f64,
}

impl Default for NoiseMitigationConfig {
    fn default() -> Self {
        Self {
            enabled: true,
            techniques: vec![
                NoiseMitigationTechnique::ZeroNoiseExtrapolation,
                NoiseMitigationTechnique::SymmetryVerification,
            ],
            noise_characterization: NoiseCharacterization::default(),
            error_threshold: 0.01,
        }
    }
}

/// Noise mitigation techniques
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum NoiseMitigationTechnique {
    /// Zero-noise extrapolation
    ZeroNoiseExtrapolation,
    /// Probabilistic error cancellation
    ProbabilisticErrorCancellation,
    /// Symmetry verification
    SymmetryVerification,
    /// Virtual distillation
    VirtualDistillation,
    /// Dynamical decoupling
    DynamicalDecoupling,
    /// Composite pulses
    CompositePulses,
}

/// Noise characterization
#[derive(Debug, Clone)]
pub struct NoiseCharacterization {
    /// Gate error rates
    pub gate_errors: HashMap<String, f64>,
    /// Readout error rates
    pub readout_errors: Vec<f64>,
    /// Coherence times
    pub coherence_times: Vec<f64>,
    /// Cross-talk matrix
    pub crosstalk_matrix: Vec<Vec<f64>>,
}

impl Default for NoiseCharacterization {
    fn default() -> Self {
        Self {
            gate_errors: HashMap::new(),
            readout_errors: Vec::new(),
            coherence_times: Vec::new(),
            crosstalk_matrix: Vec::new(),
        }
    }
}

/// Adaptive depth controller
#[derive(Debug, Clone)]
pub struct AdaptiveDepthController {
    /// Current depth
    pub current_depth: usize,
    /// Depth increment factor
    pub increment_factor: f64,
    /// Performance history
    pub performance_history: Vec<DepthPerformance>,
    /// Depth selection strategy
    pub selection_strategy: DepthSelectionStrategy,
    /// Convergence detector
    pub convergence_detector: DepthConvergenceDetector,
}

/// Performance at specific depth
#[derive(Debug, Clone)]
pub struct DepthPerformance {
    /// Depth level
    pub depth: usize,
    /// Best energy achieved
    pub best_energy: f64,
    /// Optimization iterations
    pub iterations: usize,
    /// Convergence time
    pub convergence_time: Duration,
    /// Parameter count
    pub parameter_count: usize,
    /// Improvement over previous depth
    pub improvement: f64,
}

/// Depth selection strategies
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum DepthSelectionStrategy {
    /// Conservative increment
    Conservative,
    /// Aggressive increment
    Aggressive,
    /// Performance-based
    PerformanceBased,
    /// Resource-aware
    ResourceAware,
    /// Theoretical limit guided
    TheoreticalGuided,
}

/// Depth convergence detector
#[derive(Debug, Clone)]
pub struct DepthConvergenceDetector {
    /// Energy improvement history
    pub improvement_history: Vec<f64>,
    /// Convergence threshold
    pub convergence_threshold: f64,
    /// Minimum depths to check
    pub min_depths: usize,
    /// Detection confidence
    pub confidence_level: f64,
}

/// Convergence metrics
#[derive(Debug, Clone)]
pub struct ConvergenceMetrics {
    /// Current energy
    pub current_energy: f64,
    /// Best energy found
    pub best_energy: f64,
    /// Energy gradient norm
    pub gradient_norm: f64,
    /// Parameter change magnitude
    pub parameter_change: f64,
    /// Convergence score
    pub convergence_score: f64,
    /// Iterations without improvement
    pub stagnation_count: usize,
}

/// Performance statistics for ∞-QAOA
#[derive(Debug, Clone)]
pub struct InfiniteQAOAStats {
    /// Total depths explored
    pub depths_explored: usize,
    /// Total optimization time
    pub total_time: Duration,
    /// Best energy achieved
    pub best_energy: f64,
    /// Optimal depth found
    pub optimal_depth: usize,
    /// Average convergence time per depth
    pub avg_convergence_time: Duration,
    /// Parameter optimization efficiency
    pub optimization_efficiency: f64,
}

impl InfiniteDepthQAOA {
    /// Create a new infinite-depth QAOA optimizer
    #[must_use]
    pub fn new(config: InfiniteQAOAConfig) -> Self {
        Self {
            config: config.clone(),
            parameter_history: Vec::new(),
            energy_history: Vec::new(),
            depth_progression: Vec::new(),
            convergence_metrics: ConvergenceMetrics {
                current_energy: f64::INFINITY,
                best_energy: f64::INFINITY,
                gradient_norm: f64::INFINITY,
                parameter_change: f64::INFINITY,
                convergence_score: 0.0,
                stagnation_count: 0,
            },
            depth_controller: AdaptiveDepthController {
                current_depth: config.initial_depth,
                increment_factor: 1.5,
                performance_history: Vec::new(),
                selection_strategy: DepthSelectionStrategy::PerformanceBased,
                convergence_detector: DepthConvergenceDetector {
                    improvement_history: Vec::new(),
                    convergence_threshold: 1e-6,
                    min_depths: 3,
                    confidence_level: 0.95,
                },
            },
            performance_stats: InfiniteQAOAStats {
                depths_explored: 0,
                total_time: Duration::from_secs(0),
                best_energy: f64::INFINITY,
                optimal_depth: 0,
                avg_convergence_time: Duration::from_secs(0),
                optimization_efficiency: 0.0,
            },
        }
    }

    /// Solve problem using infinite-depth QAOA
    pub fn solve<P>(&mut self, problem: &P) -> AdvancedQuantumResult<AnnealingResult<Vec<i32>>>
    where
        P: Clone + 'static,
    {
        // Convert the caller's problem to an Ising model, then optimize it.
        let ising_problem = self.convert_to_ising(problem)?;
        let solution = self.optimize(&ising_problem)?;
        match solution {
            Ok(annealing_solution) => {
                let spins: Vec<i32> = annealing_solution
                    .best_spins
                    .iter()
                    .map(|&s| i32::from(s))
                    .collect();
                Ok(Ok(spins))
            }
            Err(err) => Ok(Err(err)),
        }
    }

    /// Convert a supported problem type into an [`IsingModel`].
    ///
    /// The optimizer operates natively on Ising models. Inputs that already are an
    /// `IsingModel` (owned or by reference) are used directly, and `QuboModel`
    /// inputs are converted exactly via the standard `x = (1 + s)/2` substitution
    /// (`QuboModel::to_ising`). The constant energy offset produced by that
    /// conversion is dropped because it shifts every configuration's energy by the
    /// same amount and therefore does not change the optimal spin assignment.
    ///
    /// Any other type is rejected with an honest error rather than being silently
    /// replaced by fabricated problem data.
    fn convert_to_ising<P: 'static>(
        &self,
        problem: &P,
    ) -> Result<IsingModel, AdvancedQuantumError> {
        use std::any::Any;
        let any_problem = problem as &dyn Any;

        // Already an Ising model (owned or reference).
        if let Some(ising) = any_problem.downcast_ref::<IsingModel>() {
            return Ok(ising.clone());
        }
        if let Some(ising_ref) = any_problem.downcast_ref::<&IsingModel>() {
            return Ok((*ising_ref).clone());
        }

        // QUBO models convert exactly to the Ising picture.
        if let Some(qubo) = any_problem.downcast_ref::<QuboModel>() {
            return Ok(qubo.to_ising().0);
        }
        if let Some(qubo_ref) = any_problem.downcast_ref::<&QuboModel>() {
            return Ok((*qubo_ref).to_ising().0);
        }

        Err(AdvancedQuantumError::ParameterError(
            "Unsupported problem type for infinite-depth QAOA: expected IsingModel or QuboModel; \
             convert the problem to an IsingModel before solving"
                .to_string(),
        ))
    }

    /// Optimize using infinite-depth QAOA
    pub fn optimize(
        &mut self,
        problem: &IsingModel,
    ) -> AdvancedQuantumResult<AnnealingResult<AnnealingSolution>> {
        println!("Starting ∞-QAOA optimization");
        let start_time = Instant::now();

        let mut current_depth = self.config.initial_depth;
        let mut best_result = None;
        let mut converged = false;

        while current_depth <= self.config.max_depth && !converged {
            println!("Optimizing at depth {current_depth}");
            let depth_start_time = Instant::now();

            // Initialize parameters for current depth
            let initial_params = self.initialize_parameters(current_depth)?;

            // Optimize parameters at current depth
            let (optimized_params, energy) =
                self.optimize_at_depth(problem, current_depth, initial_params)?;

            // Record performance
            let depth_performance = DepthPerformance {
                depth: current_depth,
                best_energy: energy,
                iterations: self.config.max_iterations_per_depth, // Simplified
                convergence_time: depth_start_time.elapsed(),
                parameter_count: optimized_params.len(),
                improvement: if let Some(last_perf) =
                    self.depth_controller.performance_history.last()
                {
                    last_perf.best_energy - energy
                } else {
                    0.0
                },
            };

            self.depth_controller
                .performance_history
                .push(depth_performance);
            self.parameter_history.push(optimized_params.clone());
            self.energy_history.push(energy);
            self.depth_progression.push(current_depth);

            // Update best result
            if energy < self.convergence_metrics.best_energy {
                self.convergence_metrics.best_energy = energy;
                self.performance_stats.best_energy = energy;
                self.performance_stats.optimal_depth = current_depth;

                // Create result (simplified)
                best_result = Some(Ok(AnnealingSolution {
                    best_energy: energy,
                    best_spins: self.extract_solution_from_params(&optimized_params, problem)?,
                    repetitions: 1,
                    total_sweeps: current_depth * self.config.max_iterations_per_depth,
                    runtime: start_time.elapsed(),
                    info: format!("Infinite-depth QAOA with depth {current_depth}"),
                }));
            }

            // Check convergence
            converged = self.check_depth_convergence()?;

            // Determine next depth
            if !converged {
                current_depth = self.determine_next_depth(current_depth)?;
            }

            self.performance_stats.depths_explored += 1;
        }

        self.performance_stats.total_time = start_time.elapsed();
        self.performance_stats.avg_convergence_time = Duration::from_nanos(
            self.performance_stats.total_time.as_nanos() as u64
                / self.performance_stats.depths_explored.max(1) as u64,
        );

        println!(
            "∞-QAOA completed. Best energy: {:.6} at depth {}",
            self.performance_stats.best_energy, self.performance_stats.optimal_depth
        );

        best_result.ok_or_else(|| {
            AdvancedQuantumError::ConvergenceError("No valid result found".to_string())
        })
    }

    /// Initialize parameters for given depth
    fn initialize_parameters(&self, depth: usize) -> AdvancedQuantumResult<Vec<f64>> {
        let num_params = 2 * depth; // gamma and beta parameters

        match self.config.initialization_method {
            ParameterInitializationMethod::Random => {
                let mut rng = ChaCha8Rng::seed_from_u64(thread_rng().random());
                Ok((0..num_params)
                    .map(|_| rng.random_range(0.0..2.0 * PI))
                    .collect())
            }
            ParameterInitializationMethod::Heuristic => {
                let mut params = Vec::new();
                for i in 0..depth {
                    // Gamma parameters (problem Hamiltonian)
                    params.push(0.5 * PI / depth as f64);
                    // Beta parameters (mixer Hamiltonian)
                    params.push(0.25 * PI / depth as f64);
                }
                Ok(params)
            }
            ParameterInitializationMethod::Transfer => {
                if let Some(prev_params) = self.parameter_history.last() {
                    // Interpolate from previous depth
                    Ok(self.interpolate_parameters(prev_params, depth)?)
                } else {
                    // Fallback to heuristic
                    self.initialize_parameters_heuristic(depth)
                }
            }
            _ => {
                // Default to heuristic
                self.initialize_parameters_heuristic(depth)
            }
        }
    }

    /// Initialize parameters using heuristic method
    fn initialize_parameters_heuristic(&self, depth: usize) -> AdvancedQuantumResult<Vec<f64>> {
        let mut params = Vec::new();
        for _i in 0..depth {
            // Gamma parameters (problem Hamiltonian)
            params.push(0.5 * PI / depth as f64);
            // Beta parameters (mixer Hamiltonian)
            params.push(0.25 * PI / depth as f64);
        }
        Ok(params)
    }

    /// Interpolate parameters from previous depth
    fn interpolate_parameters(
        &self,
        prev_params: &[f64],
        new_depth: usize,
    ) -> AdvancedQuantumResult<Vec<f64>> {
        let prev_depth = prev_params.len() / 2;
        let new_param_count = 2 * new_depth;

        if new_depth <= prev_depth {
            // Truncate if new depth is smaller
            Ok(prev_params[..new_param_count].to_vec())
        } else {
            // Extend with interpolated values
            let mut new_params = prev_params.to_vec();

            for i in prev_depth..new_depth {
                // Simple interpolation strategy
                let gamma = if prev_depth > 0 {
                    prev_params[2 * (prev_depth - 1)]
                } else {
                    0.5 * PI / new_depth as f64
                };
                let beta = if prev_depth > 0 {
                    prev_params[2 * (prev_depth - 1) + 1]
                } else {
                    0.25 * PI / new_depth as f64
                };

                new_params.push(gamma * 0.8); // Scale down for stability
                new_params.push(beta * 0.8);
            }

            Ok(new_params)
        }
    }

    /// Optimize parameters at specific depth
    fn optimize_at_depth(
        &self,
        problem: &IsingModel,
        depth: usize,
        initial_params: Vec<f64>,
    ) -> AdvancedQuantumResult<(Vec<f64>, f64)> {
        let mut current_params = initial_params;
        let mut best_energy = f64::INFINITY;
        let mut best_params = current_params.clone();

        // Simple gradient-free optimization (in practice would use sophisticated methods)
        let mut rng = ChaCha8Rng::seed_from_u64(thread_rng().random());

        for iteration in 0..self.config.max_iterations_per_depth {
            // Evaluate current parameters
            let energy = self.evaluate_qaoa_energy(problem, depth, &current_params)?;

            if energy < best_energy {
                best_energy = energy;
                best_params.clone_from(&current_params);
            }

            // Simple parameter update (placeholder for actual optimization)
            for param in &mut current_params {
                *param += rng.random_range(-0.1..0.1);
                *param = param.clamp(0.0, 2.0 * PI); // Keep in valid range
            }

            // Check convergence
            if iteration > 10 && (best_energy - energy).abs() < self.config.optimization_tolerance {
                break;
            }
        }

        Ok((best_params, best_energy))
    }

    /// Evaluate QAOA energy expectation using improved quantum simulation
    fn evaluate_qaoa_energy(
        &self,
        problem: &IsingModel,
        depth: usize,
        params: &[f64],
    ) -> AdvancedQuantumResult<f64> {
        // Improved energy evaluation using quantum state evolution principles
        if params.len() != 2 * depth {
            return Err(AdvancedQuantumError::ParameterError(format!(
                "Expected {} parameters for depth {}, got {}",
                2 * depth,
                depth,
                params.len()
            )));
        }

        // For large systems, use approximation to avoid exponential memory
        if problem.num_qubits > 12 {
            return self.evaluate_qaoa_energy_approximation(problem, depth, params);
        }

        // Initialize state to |+⟩⊗n superposition
        let num_qubits = problem.num_qubits;
        let mut state_amplitudes = self.initialize_plus_state(num_qubits);

        // Apply alternating QAOA layers
        for layer in 0..depth {
            let gamma = params[2 * layer]; // Problem Hamiltonian angle
            let beta = params[2 * layer + 1]; // Mixer Hamiltonian angle

            // Apply problem Hamiltonian evolution: exp(-i * gamma * H_C)
            self.apply_problem_hamiltonian(&mut state_amplitudes, problem, gamma);

            // Apply mixer Hamiltonian evolution: exp(-i * beta * H_B)
            self.apply_mixer_hamiltonian(&mut state_amplitudes, num_qubits, beta);
        }

        // Calculate energy expectation value
        let energy = self.calculate_energy_expectation(&state_amplitudes, problem);
        Ok(energy)
    }

    /// Evaluate QAOA energy using approximation for large systems
    fn evaluate_qaoa_energy_approximation(
        &self,
        problem: &IsingModel,
        depth: usize,
        params: &[f64],
    ) -> AdvancedQuantumResult<f64> {
        // Use improved approximation for large systems
        let mut energy = 0.0;

        // Calculate approximate expectation values using quantum-inspired methods
        for i in 0..problem.num_qubits {
            if let Ok(bias) = problem.get_bias(i) {
                energy +=
                    bias * self.estimate_qubit_expectation_improved(i, params, depth, problem);
            }
        }

        for i in 0..problem.num_qubits {
            for j in (i + 1)..problem.num_qubits {
                if let Ok(coupling) = problem.get_coupling(i, j) {
                    if coupling.abs() > 1e-10 {
                        energy += coupling
                            * self.estimate_coupling_expectation_improved(
                                i, j, params, depth, problem,
                            );
                    }
                }
            }
        }

        Ok(energy)
    }

    /// Initialize state to |+⟩⊗n superposition state
    fn initialize_plus_state(&self, num_qubits: usize) -> Vec<Complex64> {
        let state_size = 1 << num_qubits; // 2^n
        let amplitude = Complex64::new(1.0 / (state_size as f64).sqrt(), 0.0);
        vec![amplitude; state_size]
    }

    /// Apply problem Hamiltonian evolution
    fn apply_problem_hamiltonian(&self, state: &mut [Complex64], problem: &IsingModel, gamma: f64) {
        let num_qubits = problem.num_qubits;
        let state_size = 1 << num_qubits;

        // For each computational basis state
        for basis_state in 0..state_size {
            if state[basis_state].norm() < 1e-12 {
                continue;
            }

            // Calculate energy of this basis state
            let mut energy = 0.0;

            // Add bias terms
            for i in 0..num_qubits {
                let spin = if (basis_state >> i) & 1 == 0 {
                    -1.0
                } else {
                    1.0
                };
                if let Ok(bias) = problem.get_bias(i) {
                    energy += bias * spin;
                }
            }

            // Add coupling terms
            for i in 0..num_qubits {
                for j in (i + 1)..num_qubits {
                    if let Ok(coupling) = problem.get_coupling(i, j) {
                        if coupling.abs() > 1e-10 {
                            let spin_i = if (basis_state >> i) & 1 == 0 {
                                -1.0
                            } else {
                                1.0
                            };
                            let spin_j = if (basis_state >> j) & 1 == 0 {
                                -1.0
                            } else {
                                1.0
                            };
                            energy += coupling * spin_i * spin_j;
                        }
                    }
                }
            }

            // Apply phase evolution: exp(-i * gamma * energy)
            let phase = Complex64::new(0.0, -gamma * energy).exp();
            state[basis_state] *= phase;
        }
    }

    /// Apply mixer Hamiltonian evolution (X rotations)
    fn apply_mixer_hamiltonian(&self, state: &mut [Complex64], num_qubits: usize, beta: f64) {
        let state_size = 1 << num_qubits;
        let mut new_state = vec![Complex64::new(0.0, 0.0); state_size];

        let cos_half_beta = (beta / 2.0).cos();
        let sin_half_beta = (beta / 2.0).sin();

        // Apply product of X rotations
        for basis_state in 0..state_size {
            if state[basis_state].norm() < 1e-12 {
                continue;
            }

            // For each qubit, apply X rotation
            let mut current_amplitude = state[basis_state];
            let mut current_state = basis_state;

            // Simplified: apply average effect of X rotations
            new_state[current_state] += current_amplitude * cos_half_beta.powi(num_qubits as i32);

            // Add contributions from flipped states (simplified)
            for qubit in 0..num_qubits {
                let flipped_state = current_state ^ (1 << qubit);
                new_state[flipped_state] += current_amplitude
                    * cos_half_beta.powi((num_qubits - 1) as i32)
                    * Complex64::new(0.0, -sin_half_beta);
            }
        }

        // Normalize
        let norm = new_state
            .iter()
            .map(scirs2_core::Complex::norm_sqr)
            .sum::<f64>()
            .sqrt();
        if norm > 1e-12 {
            for amplitude in &mut new_state {
                *amplitude /= norm;
            }
        }

        state.copy_from_slice(&new_state);
    }

    /// Calculate energy expectation value from quantum state
    fn calculate_energy_expectation(&self, state: &[Complex64], problem: &IsingModel) -> f64 {
        let num_qubits = problem.num_qubits;
        let state_size = 1 << num_qubits;
        let mut expectation = 0.0;

        for basis_state in 0..state_size {
            let probability = state[basis_state].norm_sqr();
            if probability < 1e-12 {
                continue;
            }

            let mut energy = 0.0;

            // Add bias terms
            for i in 0..num_qubits {
                let spin = if (basis_state >> i) & 1 == 0 {
                    -1.0
                } else {
                    1.0
                };
                if let Ok(bias) = problem.get_bias(i) {
                    energy += bias * spin;
                }
            }

            // Add coupling terms
            for i in 0..num_qubits {
                for j in (i + 1)..num_qubits {
                    if let Ok(coupling) = problem.get_coupling(i, j) {
                        if coupling.abs() > 1e-10 {
                            let spin_i = if (basis_state >> i) & 1 == 0 {
                                -1.0
                            } else {
                                1.0
                            };
                            let spin_j = if (basis_state >> j) & 1 == 0 {
                                -1.0
                            } else {
                                1.0
                            };
                            energy += coupling * spin_i * spin_j;
                        }
                    }
                }
            }

            expectation += probability * energy;
        }

        expectation
    }

    /// Effective local (mean-field) field acting on `qubit`.
    ///
    /// Combines the qubit's own bias with the sum of its coupling strengths to all
    /// other qubits. This is the standard mean-field local field `h_i + Σ_j J_ij`
    /// and makes the large-system estimators genuinely problem- and qubit-aware.
    fn local_field(&self, qubit: usize, problem: &IsingModel) -> f64 {
        let mut field = problem.get_bias(qubit).unwrap_or(0.0);
        for other in 0..problem.num_qubits {
            if other == qubit {
                continue;
            }
            field += problem.get_coupling(qubit, other).unwrap_or(0.0);
        }
        field
    }

    /// Mean-field single-qubit ⟨Z⟩ estimator for large systems.
    ///
    /// Propagates a single-qubit magnetization through the alternating QAOA layers
    /// using the qubit's actual effective local field, so the estimate depends on
    /// both the qubit index and the problem data.
    fn estimate_qubit_expectation_improved(
        &self,
        qubit: usize,
        params: &[f64],
        depth: usize,
        problem: &IsingModel,
    ) -> f64 {
        // Enhanced expectation value estimation using QAOA theory
        let mut state_prob_up = 0.5; // Start in equal superposition
        let local_field = self.local_field(qubit, problem);

        for layer in 0..depth {
            let gamma = params[2 * layer];
            let beta = params[2 * layer + 1];

            // Apply problem Hamiltonian effect (single-qubit mean-field approximation)
            // driven by the qubit's actual effective local field.
            state_prob_up = (0.5 * 2.0f64.mul_add(state_prob_up, -1.0))
                .mul_add((gamma * local_field).cos(), 0.5);

            // Apply mixer Hamiltonian effect (X rotation)
            let x_expectation = 2.0f64.mul_add(state_prob_up, -1.0); // Convert to [-1, 1]
            let z_expectation = (beta * x_expectation).cos();
            state_prob_up = 0.5f64.mul_add(z_expectation, 0.5);
        }

        let expectation = 2.0f64.mul_add(state_prob_up, -1.0); // Convert to Z expectation in [-1, 1]
                                                               // Bias the sign toward the energetically favorable orientation: a positive
                                                               // local field lowers energy for a -1 spin, so nudge the magnetization with
                                                               // the field direction to break the degeneracy of the symmetric estimator.
        (expectation - local_field.tanh()).tanh()
    }

    /// Improved two-qubit coupling expectation value estimation
    fn estimate_coupling_expectation_improved(
        &self,
        qubit1: usize,
        qubit2: usize,
        params: &[f64],
        depth: usize,
        problem: &IsingModel,
    ) -> f64 {
        // Enhanced two-qubit expectation using correlation functions
        let exp1 = self.estimate_qubit_expectation_improved(qubit1, params, depth, problem);
        let exp2 = self.estimate_qubit_expectation_improved(qubit2, params, depth, problem);

        // Calculate correlation based on QAOA dynamics
        let mut correlation_factor = 1.0;

        for layer in 0..depth {
            let gamma = params[2 * layer];
            let beta = params[2 * layer + 1];

            // Reduce correlation due to mixing
            correlation_factor *= (beta / 2.0).cos().powi(2);

            // Problem Hamiltonian can increase or decrease correlation
            correlation_factor *= 0.1f64.mul_add(-gamma.abs(), 1.0); // Simple approximation
        }

        // Return correlated expectation
        let independent_correlation = exp1 * exp2;
        let qaoa_correlation = correlation_factor.clamp(0.1, 1.0);

        independent_correlation * qaoa_correlation
    }

    /// Extract a discrete spin configuration from the optimized QAOA parameters.
    ///
    /// For simulable problem sizes the full QAOA state is reconstructed with the
    /// optimized angles (`|+⟩⊗n` followed by the alternating problem/mixer layers),
    /// and the returned bit-string is the most-probable measurement outcome of that
    /// state — a genuine, problem- and qubit-dependent readout. For larger systems
    /// the mean-field, qubit-aware estimator (which uses each qubit's actual local
    /// field) is used instead of an exponential-memory simulation.
    fn extract_solution_from_params(
        &self,
        params: &[f64],
        problem: &IsingModel,
    ) -> AdvancedQuantumResult<Vec<i8>> {
        let num_qubits = problem.num_qubits;
        if num_qubits == 0 {
            return Ok(Vec::new());
        }

        let depth = params.len() / 2;

        // Simulable regime: rebuild the actual QAOA state and read it out.
        if num_qubits <= 12 {
            let mut state = self.initialize_plus_state(num_qubits);
            for layer in 0..depth {
                let gamma = params[2 * layer];
                let beta = params[2 * layer + 1];
                self.apply_problem_hamiltonian(&mut state, problem, gamma);
                self.apply_mixer_hamiltonian(&mut state, num_qubits, beta);
            }
            return Ok(self.extract_spins_from_state(&state, num_qubits));
        }

        // Large-system regime: mean-field per-qubit magnetization sign, using the
        // qubit-aware estimator (which incorporates each qubit's local field).
        let mut solution = Vec::with_capacity(num_qubits);
        for qubit in 0..num_qubits {
            let expectation =
                self.estimate_qubit_expectation_improved(qubit, params, depth, problem);
            solution.push(if expectation >= 0.0 { 1 } else { -1 });
        }
        Ok(solution)
    }

    /// Read a spin configuration out of a QAOA amplitude vector.
    ///
    /// Returns the spins of the most-probable computational basis state (the mode
    /// of `|amplitude|²`). Bit `i` set corresponds to spin `+1`, cleared to `-1`,
    /// matching the convention used throughout the energy evaluation. Using the
    /// mode (rather than per-qubit marginals) correctly resolves degenerate,
    /// spin-flip-symmetric ground states instead of collapsing them to a symmetric
    /// average.
    fn extract_spins_from_state(&self, state: &[Complex64], num_qubits: usize) -> Vec<i8> {
        let mut best_index = 0usize;
        let mut best_probability = -1.0;
        for (index, amplitude) in state.iter().enumerate() {
            let probability = amplitude.norm_sqr();
            if probability > best_probability {
                best_probability = probability;
                best_index = index;
            }
        }

        (0..num_qubits)
            .map(|qubit| {
                if (best_index >> qubit) & 1 == 1 {
                    1i8
                } else {
                    -1i8
                }
            })
            .collect()
    }

    /// Check convergence across depths
    fn check_depth_convergence(&self) -> AdvancedQuantumResult<bool> {
        let history_len = self.depth_controller.performance_history.len();

        if history_len < self.depth_controller.convergence_detector.min_depths {
            return Ok(false);
        }

        // Check energy improvement trend
        let recent_improvements: Vec<f64> = self
            .depth_controller
            .performance_history
            .windows(2)
            .map(|window| window[0].best_energy - window[1].best_energy)
            .collect();

        let avg_improvement =
            recent_improvements.iter().sum::<f64>() / recent_improvements.len() as f64;

        // Converged if improvements are below threshold
        Ok(avg_improvement
            < self
                .depth_controller
                .convergence_detector
                .convergence_threshold)
    }

    /// Determine next depth to explore
    fn determine_next_depth(&self, current_depth: usize) -> AdvancedQuantumResult<usize> {
        match self.config.depth_strategy {
            DepthIncrementStrategy::Linear => Ok(current_depth + 1),
            DepthIncrementStrategy::Exponential => Ok((current_depth as f64 * 1.5) as usize),
            DepthIncrementStrategy::Adaptive => {
                // Base decision on recent performance
                if let Some(last_perf) = self.depth_controller.performance_history.last() {
                    if last_perf.improvement > 0.01 {
                        Ok(current_depth + 1) // Small increment for good improvement
                    } else {
                        Ok(current_depth + 2) // Larger increment for poor improvement
                    }
                } else {
                    Ok(current_depth + 1)
                }
            }
            DepthIncrementStrategy::GoldenRatio => Ok((current_depth as f64 * 1.618) as usize),
            DepthIncrementStrategy::Fibonacci => {
                // Simple Fibonacci-like increment
                Ok(current_depth + (current_depth / 2).max(1))
            }
        }
    }
}

/// Create default infinite-depth QAOA optimizer
#[must_use]
pub fn create_infinite_qaoa_optimizer() -> InfiniteDepthQAOA {
    InfiniteDepthQAOA::new(InfiniteQAOAConfig::default())
}

/// Create infinite-depth QAOA with custom configuration
#[must_use]
pub fn create_custom_infinite_qaoa(
    max_depth: usize,
    depth_strategy: DepthIncrementStrategy,
    initialization_method: ParameterInitializationMethod,
) -> InfiniteDepthQAOA {
    let mut config = InfiniteQAOAConfig::default();
    config.max_depth = max_depth;
    config.depth_strategy = depth_strategy;
    config.initialization_method = initialization_method;

    InfiniteDepthQAOA::new(config)
}

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

    #[test]
    fn test_infinite_qaoa_creation() {
        let optimizer = create_infinite_qaoa_optimizer();
        assert_eq!(optimizer.config.initial_depth, 1);
        assert_eq!(optimizer.config.max_depth, 100);
        assert_eq!(optimizer.depth_controller.current_depth, 1);
    }

    #[test]
    fn test_parameter_initialization() {
        let optimizer = create_infinite_qaoa_optimizer();
        let params = optimizer
            .initialize_parameters(3)
            .expect("should initialize parameters for depth 3");
        assert_eq!(params.len(), 6); // 2 * depth

        for &param in &params {
            assert!(param >= 0.0 && param <= 2.0 * PI);
        }
    }

    #[test]
    fn test_convert_to_ising_uses_actual_ising() {
        let optimizer = create_infinite_qaoa_optimizer();
        let mut ising = IsingModel::new(3);
        ising.set_bias(0, 0.75).expect("set bias");
        ising.set_coupling(1, 2, -0.5).expect("set coupling");

        let converted = optimizer
            .convert_to_ising(&ising)
            .expect("Ising input should convert");
        assert_eq!(converted.num_qubits, 3);
        assert!((converted.get_bias(0).expect("bias") - 0.75).abs() < 1e-12);
        assert!((converted.get_coupling(1, 2).expect("coupling") + 0.5).abs() < 1e-12);
    }

    #[test]
    fn test_convert_to_ising_from_qubo() {
        let optimizer = create_infinite_qaoa_optimizer();
        // Build a small QUBO and confirm the derived Ising matches to_ising().
        let mut qubo = QuboModel::new(2);
        qubo.set_linear(0, 1.0).expect("set linear");
        qubo.set_linear(1, -2.0).expect("set linear");
        qubo.set_quadratic(0, 1, 0.5).expect("set quadratic");

        let converted = optimizer
            .convert_to_ising(&qubo)
            .expect("QUBO input should convert");
        let (expected, _offset) = qubo.to_ising();

        assert_eq!(converted.num_qubits, expected.num_qubits);
        for i in 0..expected.num_qubits {
            assert!(
                (converted.get_bias(i).expect("bias") - expected.get_bias(i).expect("bias")).abs()
                    < 1e-12
            );
        }
        assert!(
            (converted.get_coupling(0, 1).expect("coupling")
                - expected.get_coupling(0, 1).expect("coupling"))
            .abs()
                < 1e-12
        );
    }

    #[test]
    fn test_convert_to_ising_rejects_unsupported_type() {
        let optimizer = create_infinite_qaoa_optimizer();
        // A plain integer is neither an IsingModel nor a QuboModel.
        let bogus: usize = 42;
        let result = optimizer.convert_to_ising(&bogus);
        assert!(
            result.is_err(),
            "unsupported problem types must return an honest error, not fabricated data"
        );
    }

    #[test]
    fn test_extract_spins_from_state_reads_amplitudes() {
        let optimizer = create_infinite_qaoa_optimizer();
        // A 2-qubit state concentrated on basis index 2 = binary 10:
        // bit 0 = 0 -> spin -1, bit 1 = 1 -> spin +1.
        let mut state = vec![Complex64::new(0.0, 0.0); 4];
        state[2] = Complex64::new(1.0, 0.0);

        let spins = optimizer.extract_spins_from_state(&state, 2);
        assert_eq!(spins, vec![-1i8, 1i8]);

        // A different mode yields different, qubit-specific spins (not uniform).
        let mut state_b = vec![Complex64::new(0.0, 0.0); 4];
        state_b[1] = Complex64::new(1.0, 0.0); // index 1 = binary 01 -> spins (+1, -1)
        let spins_b = optimizer.extract_spins_from_state(&state_b, 2);
        assert_eq!(spins_b, vec![1i8, -1i8]);
        assert_ne!(spins, spins_b);
    }

    #[test]
    fn test_extract_solution_derives_from_qaoa_state() {
        let optimizer = create_infinite_qaoa_optimizer();
        let mut problem = IsingModel::new(2);
        problem.set_bias(0, 1.5).expect("set bias");
        problem.set_bias(1, -2.0).expect("set bias");
        problem.set_coupling(0, 1, 0.5).expect("set coupling");

        // depth-2 parameters
        let params = vec![0.7, 0.3, 0.4, 0.9];

        // Independently reconstruct the QAOA state with the same primitives.
        let mut state = optimizer.initialize_plus_state(2);
        for layer in 0..2 {
            optimizer.apply_problem_hamiltonian(&mut state, &problem, params[2 * layer]);
            optimizer.apply_mixer_hamiltonian(&mut state, 2, params[2 * layer + 1]);
        }
        let expected = optimizer.extract_spins_from_state(&state, 2);

        let solution = optimizer
            .extract_solution_from_params(&params, &problem)
            .expect("solution extraction should succeed");

        // The extraction must come from the actual amplitudes, not a qubit-index
        // independent formula.
        assert_eq!(solution, expected);
        assert_eq!(solution.len(), 2);
        assert!(solution.iter().all(|&s| s == 1 || s == -1));
    }

    #[test]
    fn test_solve_reports_consistent_spins_for_ising() {
        // End-to-end: solving a 2-qubit Ising must return one spin per qubit
        // (each ±1), derived from the actual state rather than a uniform vector.
        let mut optimizer = create_custom_infinite_qaoa(
            2,
            DepthIncrementStrategy::Linear,
            ParameterInitializationMethod::Heuristic,
        );
        let mut problem = IsingModel::new(2);
        problem.set_bias(0, 1.0).expect("set bias");
        problem.set_bias(1, -1.0).expect("set bias");

        let result = optimizer.solve(&problem).expect("solve should succeed");
        let spins = result.expect("annealing should produce a solution");
        assert_eq!(spins.len(), 2);
        assert!(spins.iter().all(|&s| s == 1 || s == -1));
    }

    #[test]
    fn test_parameter_interpolation() {
        let optimizer = create_infinite_qaoa_optimizer();
        let prev_params = vec![1.0, 2.0, 3.0, 4.0]; // depth 2
        let interpolated = optimizer
            .interpolate_parameters(&prev_params, 3)
            .expect("should interpolate parameters from depth 2 to 3");

        assert_eq!(interpolated.len(), 6); // 2 * 3
        assert_eq!(interpolated[0], 1.0);
        assert_eq!(interpolated[1], 2.0);
        assert_eq!(interpolated[2], 3.0);
        assert_eq!(interpolated[3], 4.0);
        // Last two should be scaled versions
        assert!(interpolated[4] < 3.0);
        assert!(interpolated[5] < 4.0);
    }

    #[test]
    fn test_depth_increment_strategies() {
        let mut optimizer = create_infinite_qaoa_optimizer();

        // Test linear increment
        optimizer.config.depth_strategy = DepthIncrementStrategy::Linear;
        assert_eq!(
            optimizer
                .determine_next_depth(5)
                .expect("should determine next depth for linear strategy"),
            6
        );

        // Test exponential increment
        optimizer.config.depth_strategy = DepthIncrementStrategy::Exponential;
        assert_eq!(
            optimizer
                .determine_next_depth(4)
                .expect("should determine next depth for exponential strategy"),
            6
        ); // 4 * 1.5 = 6

        // Test golden ratio increment
        optimizer.config.depth_strategy = DepthIncrementStrategy::GoldenRatio;
        assert_eq!(
            optimizer
                .determine_next_depth(3)
                .expect("should determine next depth for golden ratio strategy"),
            4
        ); // 3 * 1.618 ≈ 4.85 -> 4
    }
}