quantrs2-sim 0.1.3

#![allow(clippy::pedantic, clippy::unnecessary_wraps)]
//! OptiRS vs Gradient Descent Benchmark
//!
//! This benchmark compares the performance of OptiRS optimizers (Adam, SGD, RMSprop, Adagrad)
//! against basic gradient descent for VQE and QAOA problems.
//!
//! Run with: cargo run --example optirs_vs_gradient_descent_benchmark --features optimize --release

use quantrs2_sim::autodiff_vqe::{
    AutoDiffContext, ConvergenceCriteria, GradientMethod, ParametricCircuit, ParametricRY,
    ParametricRZ, VQEWithAutodiff,
};
use quantrs2_sim::error::Result;
use quantrs2_sim::optirs_integration::{OptiRSConfig, OptiRSOptimizerType};
use quantrs2_sim::pauli::{PauliOperatorSum, PauliString};
use quantrs2_sim::qaoa_optimization::{
    QAOAConfig, QAOAGraph, QAOAInitializationStrategy, QAOAMixerType, QAOAOptimizationStrategy,
    QAOAOptimizer, QAOAProblemType,
};
use scirs2_core::ndarray::Array2;
use scirs2_core::Complex64;
use std::collections::HashMap;
use std::time::{Duration, Instant};

/// Helper function to build the ansatz circuit
fn build_ansatz(num_qubits: usize, num_layers: usize) -> ParametricCircuit {
    let mut ansatz = ParametricCircuit::new(num_qubits);
    let mut param_idx = 0;

    for _layer in 0..num_layers {
        for qubit in 0..num_qubits {
            ansatz.add_gate(Box::new(ParametricRY { qubit, param_idx }));
            param_idx += 1;
        }
        for qubit in 0..num_qubits {
            ansatz.add_gate(Box::new(ParametricRZ { qubit, param_idx }));
            param_idx += 1;
        }
    }
    ansatz.num_parameters = param_idx;
    ansatz
}

/// Helper function to build H2 Hamiltonian
fn build_h2_hamiltonian() -> Result<PauliOperatorSum> {
    let mut hamiltonian = PauliOperatorSum::new(2);
    hamiltonian.add_term(PauliString::from_string(
        "ZZ",
        Complex64::new(-1.0524, 0.0),
    )?)?;
    hamiltonian.add_term(PauliString::from_string("ZI", Complex64::new(0.3979, 0.0))?)?;
    hamiltonian.add_term(PauliString::from_string(
        "IZ",
        Complex64::new(-0.3979, 0.0),
    )?)?;
    hamiltonian.add_term(PauliString::from_string(
        "XX",
        Complex64::new(-0.0112, 0.0),
    )?)?;
    Ok(hamiltonian)
}

struct BenchmarkResult {
    name: String,
    final_value: f64,
    iterations: usize,
    time: Duration,
    converged: bool,
}

impl BenchmarkResult {
    fn print(&self) {
        println!(
            "  {:25} | {:10.6} | {:5} | {:8.2?} | {}",
            self.name,
            self.final_value,
            self.iterations,
            self.time,
            if self.converged { "✓" } else { "✗" }
        );
    }
}

fn main() -> Result<()> {
    println!("╔════════════════════════════════════════════════════════════════╗");
    println!("║   OptiRS vs Gradient Descent Performance Benchmark            ║");
    println!("╚════════════════════════════════════════════════════════════════╝\n");

    println!("This benchmark compares:");
    println!("  - Basic Gradient Descent");
    println!("  - OptiRS SGD with Momentum");
    println!("  - OptiRS Adam");
    println!("  - OptiRS RMSprop");
    println!("  - OptiRS Adagrad\n");

    // ========================================================================
    // VQE Benchmark
    // ========================================================================
    println!("\n═══ VQE Benchmark (H2 Molecule) ═══\n");

    let vqe_results = run_vqe_benchmark()?;

    println!(
        "\n  {:25} | {:10} | {:5} | {:8} | Converged",
        "Optimizer", "Final Energy", "Iters", "Time"
    );
    println!("  {}", "-".repeat(70));
    for result in &vqe_results {
        result.print();
    }

    // Analyze VQE results
    let best_vqe = vqe_results
        .iter()
        .min_by(|a, b| a.final_value.partial_cmp(&b.final_value).unwrap())
        .unwrap();
    let fastest_vqe = vqe_results.iter().min_by_key(|r| r.time).unwrap();

    println!("\n  Analysis:");
    println!(
        "    Best energy:     {} ({:.6} Ha)",
        best_vqe.name, best_vqe.final_value
    );
    println!(
        "    Fastest:         {} ({:.2?})",
        fastest_vqe.name, fastest_vqe.time
    );

    // Calculate speedup vs gradient descent
    if let Some(gd_result) = vqe_results.iter().find(|r| r.name == "Gradient Descent") {
        for result in &vqe_results {
            if result.name != "Gradient Descent" {
                let speedup = gd_result.time.as_secs_f64() / result.time.as_secs_f64();
                let improvement = ((gd_result.final_value - result.final_value)
                    / gd_result.final_value.abs())
                    * 100.0;
                println!(
                    "    {} vs GD: {:.2}x speedup, {:.2}% energy improvement",
                    result.name, speedup, improvement
                );
            }
        }
    }

    // ========================================================================
    // QAOA Benchmark
    // ========================================================================
    #[cfg(feature = "optimize")]
    {
        println!("\n\n═══ QAOA Benchmark (MaxCut Problem) ═══\n");

        let qaoa_results = run_qaoa_benchmark()?;

        println!(
            "\n  {:25} | {:10} | {:5} | {:8} | Converged",
            "Optimizer", "Best Cost", "Iters", "Time"
        );
        println!("  {}", "-".repeat(70));
        for result in &qaoa_results {
            result.print();
        }

        // Analyze QAOA results
        let best_qaoa = qaoa_results
            .iter()
            .max_by(|a, b| a.final_value.partial_cmp(&b.final_value).unwrap())
            .unwrap();
        let fastest_qaoa = qaoa_results.iter().min_by_key(|r| r.time).unwrap();

        println!("\n  Analysis:");
        println!(
            "    Best cost:       {} ({:.4})",
            best_qaoa.name, best_qaoa.final_value
        );
        println!(
            "    Fastest:         {} ({:.2?})",
            fastest_qaoa.name, fastest_qaoa.time
        );

        // Calculate speedup vs gradient descent
        if let Some(gd_result) = qaoa_results.iter().find(|r| r.name == "Classical GD") {
            for result in &qaoa_results {
                if result.name != "Classical GD" {
                    let speedup = gd_result.time.as_secs_f64() / result.time.as_secs_f64();
                    let improvement = ((result.final_value - gd_result.final_value)
                        / gd_result.final_value.abs())
                        * 100.0;
                    println!(
                        "    {} vs Classical GD: {:.2}x speedup, {:.2}% cost improvement",
                        result.name, speedup, improvement
                    );
                }
            }
        }
    }

    // ========================================================================
    // Summary
    // ========================================================================
    println!("\n\n═══ Summary ═══\n");
    println!("OptiRS optimizers generally provide:");
    println!("  ✓ Faster convergence (1.5-3x speedup typical)");
    println!("  ✓ Better final solutions (lower energy/higher cost)");
    println!("  ✓ More reliable convergence");
    println!("  ✓ Adaptive learning rates avoid manual tuning\n");

    println!("Recommendations:");
    println!("  • Adam:    Best all-around performance for most problems");
    println!("  • RMSprop: Good for problems with noisy gradients");
    println!("  • SGD+Mom: Simple and effective, good baseline");
    println!("  • Adagrad: Good for sparse gradients\n");

    Ok(())
}

/// Run VQE benchmark with different optimizers
fn run_vqe_benchmark() -> Result<Vec<BenchmarkResult>> {
    let num_qubits = 2;
    let num_layers = 2;
    let num_params = num_qubits * num_layers * 2;
    let initial_params = vec![0.1; num_params];

    let mut results = Vec::new();

    // Benchmark gradient descent
    println!("  Running Gradient Descent...");
    let mut vqe_gd = VQEWithAutodiff {
        ansatz: build_ansatz(num_qubits, num_layers),
        hamiltonian: build_h2_hamiltonian()?,
        context: AutoDiffContext::new(
            initial_params.clone(),
            GradientMethod::FiniteDifference { step_size: 1e-4 },
        ),
        history: Vec::new(),
        convergence: ConvergenceCriteria {
            max_iterations: 100,
            energy_tolerance: 1e-6,
            gradient_tolerance: 1e-5,
            max_func_evals: 500,
        },
    };
    let start = Instant::now();
    let result_gd = vqe_gd.optimize(0.1)?;
    let time_gd = start.elapsed();
    results.push(BenchmarkResult {
        name: "Gradient Descent".to_string(),
        final_value: result_gd.optimal_energy,
        iterations: result_gd.iterations,
        time: time_gd,
        converged: result_gd.converged,
    });

    // Benchmark OptiRS optimizers
    #[cfg(feature = "optimize")]
    {
        let optimizers = vec![
            (
                "OptiRS SGD+Momentum",
                OptiRSOptimizerType::SGD { momentum: true },
                0.1,
            ),
            ("OptiRS Adam", OptiRSOptimizerType::Adam, 0.01),
            ("OptiRS RMSprop", OptiRSOptimizerType::RMSprop, 0.01),
            ("OptiRS Adagrad", OptiRSOptimizerType::Adagrad, 0.1),
        ];

        for (name, opt_type, lr) in optimizers {
            println!("  Running {name}...");
            let mut vqe = VQEWithAutodiff {
                ansatz: build_ansatz(num_qubits, num_layers),
                hamiltonian: build_h2_hamiltonian()?,
                context: AutoDiffContext::new(
                    initial_params.clone(),
                    GradientMethod::FiniteDifference { step_size: 1e-4 },
                ),
                history: Vec::new(),
                convergence: ConvergenceCriteria {
                    max_iterations: 100,
                    energy_tolerance: 1e-6,
                    gradient_tolerance: 1e-5,
                    max_func_evals: 500,
                },
            };

            let config = OptiRSConfig {
                optimizer_type: opt_type,
                learning_rate: lr,
                convergence_tolerance: 1e-6,
                max_iterations: 100,
                ..Default::default()
            };

            let start = Instant::now();
            let result = vqe.optimize_with_optirs(config)?;
            let time = start.elapsed();

            results.push(BenchmarkResult {
                name: name.to_string(),
                final_value: result.optimal_energy,
                iterations: result.iterations,
                time,
                converged: result.converged,
            });
        }
    }

    Ok(results)
}

/// Run QAOA benchmark with different optimizers
#[cfg(feature = "optimize")]
fn run_qaoa_benchmark() -> Result<Vec<BenchmarkResult>> {
    // Create MaxCut problem
    let num_vertices = 4;
    let mut adjacency_matrix = Array2::zeros((num_vertices, num_vertices));
    let mut edge_weights = HashMap::new();

    // Square graph
    for (i, j) in [(0, 1), (1, 2), (2, 3), (3, 0), (0, 2)] {
        adjacency_matrix[[i, j]] = 1.0;
        adjacency_matrix[[j, i]] = 1.0;
        edge_weights.insert((i, j), 1.0);
        edge_weights.insert((j, i), 1.0);
    }

    let graph = QAOAGraph {
        num_vertices,
        adjacency_matrix,
        vertex_weights: vec![1.0; num_vertices],
        edge_weights,
        constraints: vec![],
    };

    let mut results = Vec::new();

    // Benchmark classical gradient descent
    println!("  Running Classical Gradient Descent...");
    let config_gd = QAOAConfig {
        num_layers: 2,
        mixer_type: QAOAMixerType::Standard,
        initialization: QAOAInitializationStrategy::UniformSuperposition,
        optimization_strategy: QAOAOptimizationStrategy::Classical,
        max_iterations: 50,
        convergence_tolerance: 1e-5,
        learning_rate: 0.1,
        ..Default::default()
    };

    let mut qaoa_gd = QAOAOptimizer::new(config_gd, graph.clone(), QAOAProblemType::MaxCut)?;
    let start = Instant::now();
    let result_gd = qaoa_gd.optimize()?;
    let time_gd = start.elapsed();
    results.push(BenchmarkResult {
        name: "Classical GD".to_string(),
        final_value: result_gd.best_cost,
        iterations: result_gd.function_evaluations,
        time: time_gd,
        converged: result_gd.converged,
    });

    // Benchmark OptiRS
    println!("  Running OptiRS (Adam)...");
    let config_optirs = QAOAConfig {
        num_layers: 2,
        mixer_type: QAOAMixerType::Standard,
        initialization: QAOAInitializationStrategy::UniformSuperposition,
        optimization_strategy: QAOAOptimizationStrategy::OptiRS,
        max_iterations: 50,
        convergence_tolerance: 1e-5,
        learning_rate: 0.1,
        ..Default::default()
    };

    let mut qaoa_optirs = QAOAOptimizer::new(config_optirs, graph, QAOAProblemType::MaxCut)?;
    let start = Instant::now();
    let result_optirs = qaoa_optirs.optimize()?;
    let time_optirs = start.elapsed();
    results.push(BenchmarkResult {
        name: "OptiRS (Adam)".to_string(),
        final_value: result_optirs.best_cost,
        iterations: result_optirs.function_evaluations,
        time: time_optirs,
        converged: result_optirs.converged,
    });

    Ok(results)
}