optirs_bench/

cross_framework.rs

// Cross-framework benchmarking against PyTorch and TensorFlow optimizers
//
// This module provides comprehensive benchmarking capabilities to compare
// SciRS2 optimizers against their PyTorch and TensorFlow counterparts.

use crate::error::{OptimError, Result};
use crate::TestFunction;
use scirs2_core::ndarray::Array1;
use scirs2_core::numeric::Float;
// use serde::{Deserialize, Serialize}; // Commented out for now
use std::collections::HashMap;
use std::fmt::Debug;
use std::process::Command;
use std::time::{Duration, Instant};

/// Type alias for optimizer functions taking parameters and gradients
type OptimizerFn<A> = Box<dyn Fn(&Array1<A>, &Array1<A>) -> Array1<A>>;

/// Type alias for optimizer function references (trait objects)
type OptimizerFnRef<'a, A> = &'a dyn Fn(&Array1<A>, &Array1<A>) -> Array1<A>;

/// Cross-framework benchmark configuration
#[derive(Debug, Clone)]
pub struct CrossFrameworkConfig {
    /// Enable PyTorch comparison
    pub enable_pytorch: bool,
    /// Enable TensorFlow comparison
    pub enable_tensorflow: bool,
    /// Python executable path
    pub python_path: String,
    /// Temporary directory for Python scripts
    pub temp_dir: String,
    /// Benchmark precision (f32 or f64)
    pub precision: Precision,
    /// Maximum iterations per test
    pub max_iterations: usize,
    /// Convergence tolerance
    pub tolerance: f64,
    /// Random seed for reproducibility
    pub random_seed: u64,
    /// Batch sizes to test
    pub batch_sizes: Vec<usize>,
    /// Problem dimensions to test
    pub problem_dimensions: Vec<usize>,
    /// Number of runs per test for statistical significance
    pub num_runs: usize,
}

impl Default for CrossFrameworkConfig {
    fn default() -> Self {
        Self {
            enable_pytorch: true,
            enable_tensorflow: true,
            python_path: "python3".to_string(),
            temp_dir: "/tmp/scirs2_benchmark".to_string(),
            precision: Precision::F64,
            max_iterations: 1000,
            tolerance: 1e-6,
            random_seed: 42,
            batch_sizes: vec![1, 32, 128, 512],
            problem_dimensions: vec![10, 100, 1000],
            num_runs: 5,
        }
    }
}

/// Precision options for benchmarking
#[derive(Debug, Clone, Copy)]
pub enum Precision {
    F32,
    F64,
}

/// Framework identifier
#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub enum Framework {
    SciRS2,
    PyTorch,
    TensorFlow,
}

impl std::fmt::Display for Framework {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Framework::SciRS2 => write!(f, "SciRS2"),
            Framework::PyTorch => write!(f, "PyTorch"),
            Framework::TensorFlow => write!(f, "TensorFlow"),
        }
    }
}

/// Optimizer identifier for cross-framework comparison
#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub struct OptimizerIdentifier {
    pub framework: Framework,
    pub name: String,
    pub version: Option<String>,
}

impl std::fmt::Display for OptimizerIdentifier {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        if let Some(ref version) = self.version {
            write!(f, "{}-{}-v{}", self.framework, self.name, version)
        } else {
            write!(f, "{}-{}", self.framework, self.name)
        }
    }
}

/// Comprehensive benchmark result with framework comparison
#[derive(Debug, Clone)]
pub struct CrossFrameworkBenchmarkResult<A: Float> {
    /// Test configuration
    pub config: CrossFrameworkConfig,
    /// Test function name
    pub function_name: String,
    /// Problem dimension
    pub problem_dim: usize,
    /// Batch size
    pub batch_size: usize,
    /// Results per optimizer
    pub optimizer_results: HashMap<OptimizerIdentifier, OptimizerBenchmarkSummary<A>>,
    /// Statistical comparison
    pub statistical_comparison: StatisticalComparison<A>,
    /// Performance ranking
    pub performance_ranking: Vec<(OptimizerIdentifier, f64)>,
    /// Resource usage comparison
    pub resource_usage: ResourceUsageComparison,
    /// Timestamp
    pub timestamp: std::time::Instant,
}

/// Summary statistics for an optimizer across multiple runs
#[derive(Debug, Clone)]
pub struct OptimizerBenchmarkSummary<A: Float> {
    /// Optimizer identifier
    pub optimizer: OptimizerIdentifier,
    /// Number of successful runs
    pub successful_runs: usize,
    /// Total runs attempted
    pub total_runs: usize,
    /// Success rate (0.0 to 1.0)
    pub success_rate: f64,
    /// Mean convergence time
    pub mean_convergence_time: Duration,
    /// Standard deviation of convergence time
    pub std_convergence_time: Duration,
    /// Mean final function value
    pub mean_final_value: A,
    /// Standard deviation of final function value
    pub std_final_value: A,
    /// Mean iterations to convergence
    pub mean_iterations: f64,
    /// Standard deviation of iterations
    pub std_iterations: f64,
    /// Mean final gradient norm
    pub mean_gradient_norm: A,
    /// Standard deviation of gradient norm
    pub std_gradient_norm: A,
    /// Convergence curves (one per run)
    pub convergence_curves: Vec<Vec<A>>,
    /// Memory usage statistics
    pub memory_stats: MemoryStats,
    /// GPU utilization (if applicable)
    pub gpu_utilization: Option<f64>,
}

/// Statistical comparison between optimizers
#[derive(Debug, Clone)]
pub struct StatisticalComparison<A: Float> {
    /// Pairwise t-test results for convergence time
    pub convergence_time_tests: HashMap<(OptimizerIdentifier, OptimizerIdentifier), TTestResult>,
    /// Pairwise t-test results for final function value
    pub final_value_tests: HashMap<(OptimizerIdentifier, OptimizerIdentifier), TTestResult>,
    /// ANOVA results
    pub anova_results: AnovaResult<A>,
    /// Effect sizes (Cohen's d)
    pub effect_sizes: HashMap<(OptimizerIdentifier, OptimizerIdentifier), f64>,
    /// Confidence intervals
    pub confidence_intervals: HashMap<OptimizerIdentifier, ConfidenceInterval<A>>,
}

/// T-test result for pairwise comparison
#[derive(Debug, Clone)]
pub struct TTestResult {
    /// T-statistic
    pub t_statistic: f64,
    /// P-value
    pub p_value: f64,
    /// Degrees of freedom
    pub degrees_of_freedom: f64,
    /// Is statistically significant (p < 0.05)
    pub is_significant: bool,
}

/// ANOVA result for multiple group comparison
#[derive(Debug, Clone)]
pub struct AnovaResult<A: Float> {
    /// F-statistic
    pub f_statistic: f64,
    /// P-value
    pub p_value: f64,
    /// Between-group sum of squares
    pub between_ss: A,
    /// Within-group sum of squares
    pub within_ss: A,
    /// Total sum of squares
    pub total_ss: A,
    /// Degrees of freedom between groups
    pub df_between: usize,
    /// Degrees of freedom within groups
    pub df_within: usize,
}

/// Confidence interval
#[derive(Debug, Clone)]
pub struct ConfidenceInterval<A: Float> {
    /// Lower bound
    pub lower: A,
    /// Upper bound
    pub upper: A,
    /// Confidence level (e.g., 0.95 for 95%)
    pub confidence_level: f64,
}

/// Resource usage comparison
#[derive(Debug, Clone)]
pub struct ResourceUsageComparison {
    /// Memory usage per optimizer
    pub memory_usage: HashMap<OptimizerIdentifier, MemoryStats>,
    /// CPU usage per optimizer
    pub cpu_usage: HashMap<OptimizerIdentifier, CpuStats>,
    /// GPU usage per optimizer (if applicable)
    pub gpu_usage: HashMap<OptimizerIdentifier, Option<GpuStats>>,
}

/// Memory usage statistics
#[derive(Debug, Clone)]
pub struct MemoryStats {
    /// Peak memory usage (bytes)
    pub peak_memory_bytes: usize,
    /// Average memory usage (bytes)
    pub avg_memory_bytes: usize,
    /// Memory allocations count
    pub allocation_count: usize,
    /// Memory fragmentation ratio
    pub fragmentation_ratio: f64,
}

/// CPU usage statistics
#[derive(Debug, Clone)]
pub struct CpuStats {
    /// CPU utilization percentage
    pub cpu_percent: f64,
    /// Number of CPU cores used
    pub cores_used: usize,
    /// Cache misses
    pub cache_misses: usize,
    /// Context switches
    pub context_switches: usize,
}

/// GPU usage statistics
#[derive(Debug, Clone)]
pub struct GpuStats {
    /// GPU utilization percentage
    pub gpu_percent: f64,
    /// GPU memory usage (bytes)
    pub memory_usage_bytes: usize,
    /// Kernel launches
    pub kernel_launches: usize,
    /// Average kernel execution time (microseconds)
    pub avg_kernel_time_us: f64,
}

/// Cross-framework benchmark suite
pub struct CrossFrameworkBenchmark<A: Float> {
    /// Configuration
    config: CrossFrameworkConfig,
    /// Test functions
    test_functions: Vec<TestFunction<A>>,
    /// Python script templates
    python_scripts: PythonScriptTemplates,
    /// Results storage
    results: Vec<CrossFrameworkBenchmarkResult<A>>,
}

/// Python script templates for external framework benchmarking
struct PythonScriptTemplates {
    /// PyTorch optimizer script template
    pytorch_template: String,
    /// TensorFlow optimizer script template
    tensorflow_template: String,
}

impl<A: Float + Debug + Send + Sync> CrossFrameworkBenchmark<A> {
    /// Create a new cross-framework benchmark suite
    pub fn new(config: CrossFrameworkConfig) -> Result<Self> {
        let python_scripts = PythonScriptTemplates::new();

        // Create temporary directory
        std::fs::create_dir_all(&config.temp_dir).map_err(|e| {
            OptimError::InvalidConfig(format!("Failed to create temp directory: {}", e))
        })?;

        Ok(Self {
            config,
            test_functions: Vec::new(),
            python_scripts,
            results: Vec::new(),
        })
    }

    /// Add a test function to the benchmark suite
    pub fn add_test_function(&mut self, test_function: TestFunction<A>) {
        self.test_functions.push(test_function);
    }

    /// Add standard optimization test functions
    pub fn add_standard_test_functions(&mut self) {
        // Quadratic function
        self.add_test_function(TestFunction {
            name: "Quadratic".to_string(),
            dimension: 10,
            function: Box::new(|x: &Array1<A>| x.mapv(|val| val * val).sum()),
            gradient: Box::new(|x: &Array1<A>| {
                x.mapv(|val| A::from(2.0).expect("unwrap failed") * val)
            }),
            optimal_value: Some(A::zero()),
            optimal_point: Some(Array1::zeros(10)),
        });

        // Rosenbrock function
        self.add_test_function(TestFunction {
            name: "Rosenbrock".to_string(),
            dimension: 2,
            function: Box::new(|x: &Array1<A>| {
                let a = A::one();
                let b = A::from(100.0).expect("unwrap failed");
                let term1 = (a - x[0]) * (a - x[0]);
                let term2 = b * (x[1] - x[0] * x[0]) * (x[1] - x[0] * x[0]);
                term1 + term2
            }),
            gradient: Box::new(|x: &Array1<A>| {
                let a = A::one();
                let b = A::from(100.0).expect("unwrap failed");
                let grad_x = A::from(-2.0).expect("unwrap failed") * (a - x[0])
                    - A::from(4.0).expect("unwrap failed") * b * x[0] * (x[1] - x[0] * x[0]);
                let grad_y = A::from(2.0).expect("unwrap failed") * b * (x[1] - x[0] * x[0]);
                Array1::from_vec(vec![grad_x, grad_y])
            }),
            optimal_value: Some(A::zero()),
            optimal_point: Some(Array1::from_vec(vec![A::one(), A::one()])),
        });

        // Beale function
        self.add_test_function(TestFunction {
            name: "Beale".to_string(),
            dimension: 2,
            function: Box::new(|x: &Array1<A>| {
                let x1 = x[0];
                let x2 = x[1];
                let term1 = (A::from(1.5).expect("unwrap failed") - x1 + x1 * x2)
                    * (A::from(1.5).expect("unwrap failed") - x1 + x1 * x2);
                let term2 = (A::from(2.25).expect("unwrap failed") - x1 + x1 * x2 * x2)
                    * (A::from(2.25).expect("unwrap failed") - x1 + x1 * x2 * x2);
                let term3 = (A::from(2.625).expect("unwrap failed") - x1 + x1 * x2 * x2 * x2)
                    * (A::from(2.625).expect("unwrap failed") - x1 + x1 * x2 * x2 * x2);
                term1 + term2 + term3
            }),
            gradient: Box::new(|x: &Array1<A>| {
                let x1 = x[0];
                let x2 = x[1];
                let dx1 = A::from(2.0).expect("unwrap failed")
                    * (A::from(1.5).expect("unwrap failed") - x1 + x1 * x2)
                    * (x2 - A::one())
                    + A::from(2.0).expect("unwrap failed")
                        * (A::from(2.25).expect("unwrap failed") - x1 + x1 * x2 * x2)
                        * (x2 * x2 - A::one())
                    + A::from(2.0).expect("unwrap failed")
                        * (A::from(2.625).expect("unwrap failed") - x1 + x1 * x2 * x2 * x2)
                        * (x2 * x2 * x2 - A::one());
                let dx2 = A::from(2.0).expect("unwrap failed")
                    * (A::from(1.5).expect("unwrap failed") - x1 + x1 * x2)
                    * x1
                    + A::from(2.0).expect("unwrap failed")
                        * (A::from(2.25).expect("unwrap failed") - x1 + x1 * x2 * x2)
                        * (A::from(2.0).expect("unwrap failed") * x1 * x2)
                    + A::from(2.0).expect("unwrap failed")
                        * (A::from(2.625).expect("unwrap failed") - x1 + x1 * x2 * x2 * x2)
                        * (A::from(3.0).expect("unwrap failed") * x1 * x2 * x2);
                Array1::from_vec(vec![dx1, dx2])
            }),
            optimal_value: Some(A::zero()),
            optimal_point: Some(Array1::from_vec(vec![
                A::from(3.0).expect("unwrap failed"),
                A::from(0.5).expect("unwrap failed"),
            ])),
        });
    }

    /// Run comprehensive cross-framework benchmark
    pub fn run_comprehensive_benchmark(
        &mut self,
        scirs2_optimizers: Vec<(String, OptimizerFn<A>)>,
    ) -> Result<Vec<CrossFrameworkBenchmarkResult<A>>> {
        let mut all_results = Vec::new();

        for test_function in &self.test_functions {
            for &problem_dim in &self.config.problem_dimensions {
                for &batch_size in &self.config.batch_sizes {
                    let result = self.run_single_benchmark(
                        test_function,
                        problem_dim,
                        batch_size,
                        &scirs2_optimizers,
                    )?;
                    all_results.push(result);
                }
            }
        }

        self.results.extend(all_results.clone());
        Ok(all_results)
    }

    /// Run benchmark for a single configuration
    fn run_single_benchmark(
        &self,
        test_function: &TestFunction<A>,
        problem_dim: usize,
        batch_size: usize,
        scirs2_optimizers: &[(String, OptimizerFn<A>)],
    ) -> Result<CrossFrameworkBenchmarkResult<A>> {
        let mut optimizer_results = HashMap::new();

        // Run SciRS2 _optimizers
        for (name, optimizer) in scirs2_optimizers {
            let identifier = OptimizerIdentifier {
                framework: Framework::SciRS2,
                name: name.clone(),
                version: Some("0.1.0".to_string()),
            };

            let summary =
                self.benchmark_scirs2_optimizer(test_function, problem_dim, batch_size, optimizer)?;
            optimizer_results.insert(identifier, summary);
        }

        // Run PyTorch _optimizers
        if self.config.enable_pytorch {
            let pytorch_results =
                self.benchmark_pytorch_optimizers(test_function, problem_dim, batch_size)?;
            optimizer_results.extend(pytorch_results);
        }

        // Run TensorFlow _optimizers
        if self.config.enable_tensorflow {
            let tensorflow_results =
                self.benchmark_tensorflow_optimizers(test_function, problem_dim, batch_size)?;
            optimizer_results.extend(tensorflow_results);
        }

        // Perform statistical analysis
        let statistical_comparison = self.perform_statistical_analysis(&optimizer_results)?;

        // Rank _optimizers by performance
        let performance_ranking = self.rank_optimizers(&optimizer_results);

        // Analyze resource usage
        let resource_usage = self.analyze_resource_usage(&optimizer_results);

        Ok(CrossFrameworkBenchmarkResult {
            config: self.config.clone(),
            function_name: test_function.name.clone(),
            problem_dim,
            batch_size,
            optimizer_results,
            statistical_comparison,
            performance_ranking,
            resource_usage,
            timestamp: std::time::Instant::now(),
        })
    }

    /// Benchmark SciRS2 optimizer
    fn benchmark_scirs2_optimizer(
        &self,
        test_function: &TestFunction<A>,
        problem_dim: usize,
        _batch_size: usize,
        optimizer: &OptimizerFn<A>,
    ) -> Result<OptimizerBenchmarkSummary<A>> {
        let mut convergence_times = Vec::new();
        let mut final_values = Vec::new();
        let mut iterations_counts = Vec::new();
        let mut gradient_norms = Vec::new();
        let mut convergence_curves = Vec::new();
        let mut successful_runs = 0;

        for run in 0..self.config.num_runs {
            // Set random seed for reproducibility
            let mut rng_seed = self.config.random_seed + run as u64;

            // Initialize parameters
            let mut x = Array1::from_vec(
                (0..problem_dim)
                    .map(|_| {
                        rng_seed = rng_seed.wrapping_mul(1103515245).wrapping_add(12345);
                        A::from((rng_seed % 1000) as f64 / 1000.0 - 0.5).expect("unwrap failed")
                    })
                    .collect(),
            );

            let start_time = Instant::now();
            let mut convergence_curve = Vec::new();
            let mut converged = false;

            for iteration in 0..self.config.max_iterations {
                let f_val = (test_function.function)(&x);
                let grad = (test_function.gradient)(&x);
                let grad_norm = grad.mapv(|g| g * g).sum().sqrt();

                convergence_curve.push(f_val);

                // Check convergence
                if grad_norm.to_f64().unwrap_or(f64::INFINITY) < self.config.tolerance {
                    let elapsed = start_time.elapsed();
                    convergence_times.push(elapsed);
                    final_values.push(f_val);
                    iterations_counts.push(iteration as f64);
                    gradient_norms.push(grad_norm);
                    convergence_curves.push(convergence_curve.clone());
                    successful_runs += 1;
                    converged = true;
                    break;
                }

                // Perform optimization step
                x = optimizer(&x, &grad);
            }

            // If didn't converge, record final state
            if !converged {
                let elapsed = start_time.elapsed();
                let f_val = (test_function.function)(&x);
                let grad = (test_function.gradient)(&x);
                let grad_norm = grad.mapv(|g| g * g).sum().sqrt();

                convergence_times.push(elapsed);
                final_values.push(f_val);
                iterations_counts.push(self.config.max_iterations as f64);
                gradient_norms.push(grad_norm);
                convergence_curves.push(convergence_curve);
            }
        }

        // Calculate statistics
        let success_rate = successful_runs as f64 / self.config.num_runs as f64;

        let mean_convergence_time = if !convergence_times.is_empty() {
            convergence_times.iter().sum::<Duration>() / convergence_times.len() as u32
        } else {
            Duration::from_secs(0)
        };

        let mean_final_value = if !final_values.is_empty() {
            final_values.iter().fold(A::zero(), |acc, &x| acc + x)
                / A::from(final_values.len()).expect("unwrap failed")
        } else {
            A::zero()
        };

        let mean_iterations = if !iterations_counts.is_empty() {
            iterations_counts.iter().sum::<f64>() / iterations_counts.len() as f64
        } else {
            0.0
        };

        let mean_gradient_norm = if !gradient_norms.is_empty() {
            gradient_norms.iter().fold(A::zero(), |acc, &x| acc + x)
                / A::from(gradient_norms.len()).expect("unwrap failed")
        } else {
            A::zero()
        };

        // Calculate standard deviations
        let std_convergence_time =
            self.calculate_duration_std(&convergence_times, mean_convergence_time);
        let std_final_value = self.calculate_std(&final_values, mean_final_value);
        let std_iterations = self.calculate_f64std(&iterations_counts, mean_iterations);
        let std_gradient_norm = self.calculate_std(&gradient_norms, mean_gradient_norm);

        Ok(OptimizerBenchmarkSummary {
            optimizer: OptimizerIdentifier {
                framework: Framework::SciRS2,
                name: "SciRS2".to_string(),
                version: Some("0.1.0".to_string()),
            },
            successful_runs,
            total_runs: self.config.num_runs,
            success_rate,
            mean_convergence_time,
            std_convergence_time,
            mean_final_value,
            std_final_value,
            mean_iterations,
            std_iterations,
            mean_gradient_norm,
            std_gradient_norm,
            convergence_curves,
            memory_stats: MemoryStats {
                peak_memory_bytes: 0,
                avg_memory_bytes: 0,
                allocation_count: 0,
                fragmentation_ratio: 0.0,
            },
            gpu_utilization: None,
        })
    }

    /// Benchmark PyTorch optimizers
    fn benchmark_pytorch_optimizers(
        &self,
        test_function: &TestFunction<A>,
        problem_dim: usize,
        batch_size: usize,
    ) -> Result<HashMap<OptimizerIdentifier, OptimizerBenchmarkSummary<A>>> {
        let script_path = format!("{}/pytorch_benchmark.py", self.config.temp_dir);

        // Write PyTorch benchmark script
        let script_content = self.python_scripts.generate_pytorch_script(
            &test_function.name,
            problem_dim,
            batch_size,
            &self.config,
        );

        std::fs::write(&script_path, script_content).map_err(|e| {
            OptimError::InvalidConfig(format!("Failed to write PyTorch script: {}", e))
        })?;

        // Execute PyTorch benchmark
        let output = Command::new(&self.config.python_path)
            .arg(&script_path)
            .output()
            .map_err(|e| {
                OptimError::InvalidConfig(format!("Failed to execute PyTorch benchmark: {}", e))
            })?;

        if !output.status.success() {
            return Err(OptimError::InvalidConfig(format!(
                "PyTorch benchmark failed: {}",
                String::from_utf8_lossy(&output.stderr)
            )));
        }

        // Parse results - simplified for now without serde_json
        let _results_json = String::from_utf8_lossy(&output.stdout);
        // let results: HashMap<String, serde_json::Value> = serde_json:::from_str(&results_json)
        //     .map_err(|e| OptimError::InvalidConfig(format!("Failed to parse PyTorch results: {}", e)))?;
        let results: HashMap<String, HashMap<String, f64>> = HashMap::new(); // Placeholder

        let mut optimizer_results = HashMap::new();

        for (optimizer_name, result_data) in results {
            let identifier = OptimizerIdentifier {
                framework: Framework::PyTorch,
                name: optimizer_name,
                version: Some("2.0".to_string()),
            };

            let summary = self.parse_python_results(identifier.clone(), &result_data)?;
            optimizer_results.insert(identifier, summary);
        }

        Ok(optimizer_results)
    }

    /// Benchmark TensorFlow optimizers
    fn benchmark_tensorflow_optimizers(
        &self,
        test_function: &TestFunction<A>,
        problem_dim: usize,
        batch_size: usize,
    ) -> Result<HashMap<OptimizerIdentifier, OptimizerBenchmarkSummary<A>>> {
        let script_path = format!("{}/tensorflow_benchmark.py", self.config.temp_dir);

        // Write TensorFlow benchmark script
        let script_content = self.python_scripts.generate_tensorflow_script(
            &test_function.name,
            problem_dim,
            batch_size,
            &self.config,
        );

        std::fs::write(&script_path, script_content).map_err(|e| {
            OptimError::InvalidConfig(format!("Failed to write TensorFlow script: {}", e))
        })?;

        // Execute TensorFlow benchmark
        let output = Command::new(&self.config.python_path)
            .arg(&script_path)
            .output()
            .map_err(|e| {
                OptimError::InvalidConfig(format!("Failed to execute TensorFlow benchmark: {}", e))
            })?;

        if !output.status.success() {
            return Err(OptimError::InvalidConfig(format!(
                "TensorFlow benchmark failed: {}",
                String::from_utf8_lossy(&output.stderr)
            )));
        }

        // Parse results - simplified for now without serde_json
        let _results_json = String::from_utf8_lossy(&output.stdout);
        // let results: HashMap<String, serde_json::Value> = serde_json:::from_str(&results_json)
        //     .map_err(|e| OptimError::InvalidConfig(format!("Failed to parse TensorFlow results: {}", e)))?;
        let results: HashMap<String, HashMap<String, f64>> = HashMap::new(); // Placeholder

        let mut optimizer_results = HashMap::new();

        for (optimizer_name, result_data) in results {
            let identifier = OptimizerIdentifier {
                framework: Framework::TensorFlow,
                name: optimizer_name,
                version: Some("2.12".to_string()),
            };

            let summary = self.parse_python_results(identifier.clone(), &result_data)?;
            optimizer_results.insert(identifier, summary);
        }

        Ok(optimizer_results)
    }

    /// Parse Python benchmark results
    fn parse_python_results(
        &self,
        identifier: OptimizerIdentifier,
        result_data: &HashMap<String, f64>,
    ) -> Result<OptimizerBenchmarkSummary<A>> {
        // Simplified parsing - using default values for now
        let successful_runs = *result_data.get("successful_runs").unwrap_or(&0.0) as usize;
        let total_runs = *result_data.get("total_runs").unwrap_or(&5.0) as usize;
        let success_rate = *result_data.get("success_rate").unwrap_or(&1.0);

        let mean_convergence_time_ms = *result_data
            .get("mean_convergence_time_ms")
            .unwrap_or(&100.0);
        let mean_convergence_time = Duration::from_millis(mean_convergence_time_ms as u64);

        let std_convergence_time_ms = *result_data.get("std_convergence_time_ms").unwrap_or(&10.0);
        let std_convergence_time = Duration::from_millis(std_convergence_time_ms as u64);

        let mean_final_value =
            A::from(*result_data.get("mean_final_value").unwrap_or(&0.1)).expect("unwrap failed");
        let std_final_value =
            A::from(*result_data.get("std_final_value").unwrap_or(&0.01)).expect("unwrap failed");

        let mean_iterations = *result_data.get("mean_iterations").unwrap_or(&100.0);
        let std_iterations = *result_data.get("std_iterations").unwrap_or(&10.0);

        let mean_gradient_norm = A::from(*result_data.get("mean_gradient_norm").unwrap_or(&0.01))
            .expect("unwrap failed");
        let std_gradient_norm = A::from(*result_data.get("std_gradient_norm").unwrap_or(&0.001))
            .expect("unwrap failed");

        // Simplified convergence curves
        let convergence_curves: Vec<Vec<A>> = vec![vec![mean_final_value; 100]; total_runs];

        Ok(OptimizerBenchmarkSummary {
            optimizer: identifier,
            successful_runs,
            total_runs,
            success_rate,
            mean_convergence_time,
            std_convergence_time,
            mean_final_value,
            std_final_value,
            mean_iterations,
            std_iterations,
            mean_gradient_norm,
            std_gradient_norm,
            convergence_curves,
            memory_stats: MemoryStats {
                peak_memory_bytes: *result_data.get("peak_memory_bytes").unwrap_or(&1000000.0)
                    as usize,
                avg_memory_bytes: *result_data.get("avg_memory_bytes").unwrap_or(&500000.0)
                    as usize,
                allocation_count: *result_data.get("allocation_count").unwrap_or(&100.0) as usize,
                fragmentation_ratio: *result_data.get("fragmentation_ratio").unwrap_or(&0.1),
            },
            gpu_utilization: result_data.get("gpu_utilization").copied(),
        })
    }

    /// Perform statistical analysis
    fn perform_statistical_analysis(
        &self,
        results: &HashMap<OptimizerIdentifier, OptimizerBenchmarkSummary<A>>,
    ) -> Result<StatisticalComparison<A>> {
        let mut convergence_time_tests = HashMap::new();
        let mut final_value_tests = HashMap::new();
        let mut effect_sizes = HashMap::new();
        let mut confidence_intervals = HashMap::new();

        // Pairwise comparisons
        let optimizers: Vec<_> = results.keys().collect();
        for i in 0..optimizers.len() {
            for j in (i + 1)..optimizers.len() {
                let opt1 = optimizers[i];
                let opt2 = optimizers[j];

                let result1 = &results[opt1];
                let result2 = &results[opt2];

                // T-test for convergence time
                let time_test = self.perform_t_test(
                    &self.extract_convergence_times(&result1.convergence_curves),
                    &self.extract_convergence_times(&result2.convergence_curves),
                );
                convergence_time_tests.insert((opt1.clone(), opt2.clone()), time_test);

                // T-test for final values
                let final_values1 = self.extract_final_values(&result1.convergence_curves);
                let final_values2 = self.extract_final_values(&result2.convergence_curves);
                let value_test = self.perform_t_test(&final_values1, &final_values2);
                final_value_tests.insert((opt1.clone(), opt2.clone()), value_test);

                // Effect size (Cohen's d)
                let effect_size = self.calculate_cohens_d(&final_values1, &final_values2);
                effect_sizes.insert((opt1.clone(), opt2.clone()), effect_size);
            }

            // Confidence intervals
            let result = &results[optimizers[i]];
            let final_values = self.extract_final_values(&result.convergence_curves);
            let ci = self.calculate_confidence_interval(&final_values, 0.95);
            confidence_intervals.insert(optimizers[i].clone(), ci);
        }

        // ANOVA
        let all_final_values: Vec<Vec<f64>> = results
            .values()
            .map(|result| self.extract_final_values(&result.convergence_curves))
            .collect();
        let anova_results = self.perform_anova(&all_final_values);

        Ok(StatisticalComparison {
            convergence_time_tests,
            final_value_tests,
            anova_results,
            effect_sizes,
            confidence_intervals,
        })
    }

    /// Rank optimizers by performance
    fn rank_optimizers(
        &self,
        results: &HashMap<OptimizerIdentifier, OptimizerBenchmarkSummary<A>>,
    ) -> Vec<(OptimizerIdentifier, f64)> {
        let mut rankings: Vec<_> = results
            .iter()
            .map(|(identifier, summary)| {
                // Composite score: success_rate * (1 / mean_final_value) * (1 / mean_convergence_time)
                let time_factor = 1.0 / (summary.mean_convergence_time.as_millis() as f64 + 1.0);
                let value_factor = 1.0 / (summary.mean_final_value.to_f64().unwrap_or(1.0) + 1e-10);
                let score = summary.success_rate * time_factor * value_factor;
                (identifier.clone(), score)
            })
            .collect();

        rankings.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
        rankings
    }

    /// Analyze resource usage
    fn analyze_resource_usage(
        &self,
        results: &HashMap<OptimizerIdentifier, OptimizerBenchmarkSummary<A>>,
    ) -> ResourceUsageComparison {
        let memory_usage = results
            .iter()
            .map(|(id, summary)| (id.clone(), summary.memory_stats.clone()))
            .collect();

        let cpu_usage = results
            .iter()
            .map(|(id_, summary)| {
                (
                    id_.clone(),
                    CpuStats {
                        cpu_percent: 0.0, // Would be measured during actual benchmarking
                        cores_used: 1,
                        cache_misses: 0,
                        context_switches: 0,
                    },
                )
            })
            .collect();

        let gpu_usage = results
            .iter()
            .map(|(id, summary)| {
                let gpu_stats = if summary.gpu_utilization.is_some() {
                    Some(GpuStats {
                        gpu_percent: summary.gpu_utilization.unwrap_or(0.0),
                        memory_usage_bytes: 0,
                        kernel_launches: 0,
                        avg_kernel_time_us: 0.0,
                    })
                } else {
                    None
                };
                (id.clone(), gpu_stats)
            })
            .collect();

        ResourceUsageComparison {
            memory_usage,
            cpu_usage,
            gpu_usage,
        }
    }

    /// Generate comprehensive benchmark report
    pub fn generate_comprehensive_report(&self) -> String {
        let mut report = String::new();

        report.push_str("# Cross-Framework Optimizer Benchmark Report\n\n");
        report.push_str(&format!(
            "Generated: {:?}\n\n",
            std::time::SystemTime::now()
        ));

        if self.results.is_empty() {
            report.push_str("No benchmark results available.\n");
            return report;
        }

        // Executive summary
        report.push_str("## Executive Summary\n\n");
        report.push_str(&format!(
            "Total test configurations: {}\n",
            self.results.len()
        ));

        // Framework coverage
        let frameworks: std::collections::HashSet<_> = self
            .results
            .iter()
            .flat_map(|result| result.optimizer_results.keys())
            .map(|id| &id.framework)
            .collect();
        report.push_str(&format!("Frameworks tested: {:?}\n\n", frameworks));

        // Performance rankings
        report.push_str("## Overall Performance Rankings\n\n");
        for result in &self.results {
            report.push_str(&format!(
                "### {} ({}D, batch={})\n\n",
                result.function_name, result.problem_dim, result.batch_size
            ));

            for (rank, (optimizer, score)) in result.performance_ranking.iter().enumerate() {
                report.push_str(&format!(
                    "{}. {} - Score: {:.6}\n",
                    rank + 1,
                    optimizer,
                    score
                ));
            }
            report.push('\n');
        }

        // Statistical significance
        report.push_str("## Statistical Analysis\n\n");
        for result in &self.results {
            report.push_str(&format!("### {} Results\n\n", result.function_name));

            // ANOVA results
            let anova = &result.statistical_comparison.anova_results;
            report.push_str(&format!(
                "ANOVA F-statistic: {:.4}, p-value: {:.6}\n",
                anova.f_statistic, anova.p_value
            ));

            if anova.p_value < 0.05 {
                report.push_str(
                    "**Statistically significant differences found between optimizers.**\n\n",
                );
            } else {
                report.push_str("No statistically significant differences found.\n\n");
            }
        }

        report
    }

    // Utility functions for statistical calculations

    /// Calculate standard deviation for Duration values
    fn calculate_duration_std(&self, values: &[Duration], mean: Duration) -> Duration {
        if values.len() <= 1 {
            return Duration::from_millis(0);
        }

        let variance = values
            .iter()
            .map(|&v| {
                let diff = v.as_millis() as i64 - mean.as_millis() as i64;
                (diff * diff) as f64
            })
            .sum::<f64>()
            / (values.len() - 1) as f64;

        Duration::from_millis(variance.sqrt() as u64)
    }

    /// Calculate standard deviation for Float values
    fn calculate_std(&self, values: &[A], mean: A) -> A {
        if values.len() <= 1 {
            return A::zero();
        }

        let variance = values
            .iter()
            .map(|&v| (v - mean) * (v - mean))
            .fold(A::zero(), |acc, x| acc + x)
            / A::from(values.len() - 1).expect("unwrap failed");

        variance.sqrt()
    }

    /// Calculate standard deviation for f64 values
    fn calculate_f64std(&self, values: &[f64], mean: f64) -> f64 {
        if values.len() <= 1 {
            return 0.0;
        }

        let variance = values.iter().map(|&v| (v - mean) * (v - mean)).sum::<f64>()
            / (values.len() - 1) as f64;

        variance.sqrt()
    }

    /// Extract convergence times from convergence curves
    fn extract_convergence_times(&self, curves: &[Vec<A>]) -> Vec<f64> {
        curves
            .iter()
            .map(|curve| curve.len() as f64) // Proxy for convergence time
            .collect()
    }

    /// Extract final values from convergence curves
    fn extract_final_values(&self, curves: &[Vec<A>]) -> Vec<f64> {
        curves
            .iter()
            .filter_map(|curve| curve.last())
            .map(|&val| val.to_f64().unwrap_or(0.0))
            .collect()
    }

    /// Perform t-test between two samples
    fn perform_t_test(&self, sample1: &[f64], sample2: &[f64]) -> TTestResult {
        if sample1.is_empty() || sample2.is_empty() {
            return TTestResult {
                t_statistic: 0.0,
                p_value: 1.0,
                degrees_of_freedom: 0.0,
                is_significant: false,
            };
        }

        let n1 = sample1.len() as f64;
        let n2 = sample2.len() as f64;

        let mean1 = sample1.iter().sum::<f64>() / n1;
        let mean2 = sample2.iter().sum::<f64>() / n2;

        let var1 = sample1.iter().map(|&x| (x - mean1).powi(2)).sum::<f64>() / (n1 - 1.0);
        let var2 = sample2.iter().map(|&x| (x - mean2).powi(2)).sum::<f64>() / (n2 - 1.0);

        let pooled_se = ((var1 / n1) + (var2 / n2)).sqrt();
        let t_statistic = (mean1 - mean2) / pooled_se;
        let degrees_of_freedom = n1 + n2 - 2.0;

        // Simplified p-value calculation (two-tailed)
        let p_value = 2.0 * (1.0 - self.t_distribution_cdf(t_statistic.abs(), degrees_of_freedom));

        TTestResult {
            t_statistic,
            p_value,
            degrees_of_freedom,
            is_significant: p_value < 0.05,
        }
    }

    /// Simplified t-distribution CDF approximation
    fn t_distribution_cdf(&self, t: f64, df: f64) -> f64 {
        // Simple approximation - in practice would use proper statistical library
        let x = t / (t * t + df).sqrt();
        0.5 + 0.5 * x / (1.0 + 0.33 * x * x)
    }

    /// Calculate Cohen's d effect size
    fn calculate_cohens_d(&self, sample1: &[f64], sample2: &[f64]) -> f64 {
        if sample1.is_empty() || sample2.is_empty() {
            return 0.0;
        }

        let mean1 = sample1.iter().sum::<f64>() / sample1.len() as f64;
        let mean2 = sample2.iter().sum::<f64>() / sample2.len() as f64;

        let var1 =
            sample1.iter().map(|&x| (x - mean1).powi(2)).sum::<f64>() / (sample1.len() - 1) as f64;
        let var2 =
            sample2.iter().map(|&x| (x - mean2).powi(2)).sum::<f64>() / (sample2.len() - 1) as f64;

        let pooled_std = ((var1 + var2) / 2.0).sqrt();
        (mean1 - mean2) / pooled_std
    }

    /// Calculate confidence interval
    fn calculate_confidence_interval(
        &self,
        values: &[f64],
        confidence_level: f64,
    ) -> ConfidenceInterval<A> {
        if values.is_empty() {
            return ConfidenceInterval {
                lower: A::zero(),
                upper: A::zero(),
                confidence_level,
            };
        }

        let mean = values.iter().sum::<f64>() / values.len() as f64;
        let std_err = self.calculate_f64std(values, mean) / (values.len() as f64).sqrt();

        // Simplified critical value (should use proper t-distribution)
        let _alpha = 1.0 - confidence_level;
        let critical_value = 1.96; // Approximate for 95% confidence

        let margin_of_error = critical_value * std_err;

        ConfidenceInterval {
            lower: A::from(mean - margin_of_error).expect("unwrap failed"),
            upper: A::from(mean + margin_of_error).expect("unwrap failed"),
            confidence_level,
        }
    }

    /// Perform ANOVA
    fn perform_anova(&self, groups: &[Vec<f64>]) -> AnovaResult<A> {
        if groups.len() < 2 {
            return AnovaResult {
                f_statistic: 0.0,
                p_value: 1.0,
                between_ss: A::zero(),
                within_ss: A::zero(),
                total_ss: A::zero(),
                df_between: 0,
                df_within: 0,
            };
        }

        let total_n: usize = groups.iter().map(|g| g.len()).sum();
        let grand_mean = groups.iter().flat_map(|g| g.iter()).sum::<f64>() / total_n as f64;

        // Between-group sum of squares
        let between_ss = groups
            .iter()
            .map(|group| {
                let group_mean = group.iter().sum::<f64>() / group.len() as f64;
                group.len() as f64 * (group_mean - grand_mean).powi(2)
            })
            .sum::<f64>();

        // Within-group sum of squares
        let within_ss = groups
            .iter()
            .flat_map(|group| {
                let group_mean = group.iter().sum::<f64>() / group.len() as f64;
                group.iter().map(move |&x| (x - group_mean).powi(2))
            })
            .sum::<f64>();

        let total_ss = between_ss + within_ss;
        let df_between = groups.len() - 1;
        let df_within = total_n - groups.len();

        let ms_between = between_ss / df_between as f64;
        let ms_within = within_ss / df_within as f64;

        let f_statistic = ms_between / ms_within;

        // Simplified p-value calculation
        let p_value = if f_statistic > 3.0 { 0.01 } else { 0.1 }; // Very rough approximation

        AnovaResult {
            f_statistic,
            p_value,
            between_ss: A::from(between_ss).expect("unwrap failed"),
            within_ss: A::from(within_ss).expect("unwrap failed"),
            total_ss: A::from(total_ss).expect("unwrap failed"),
            df_between,
            df_within,
        }
    }
}

impl PythonScriptTemplates {
    fn new() -> Self {
        Self {
            pytorch_template: r#"# PyTorch benchmark template
import torch
import torch.optim as optim

# Configuration
FUNCTION_NAME = "{{FUNCTION_NAME}}"
PROBLEM_DIM = {{PROBLEM_DIM}}
BATCH_SIZE = {{BATCH_SIZE}}
MAX_ITERATIONS = {{MAX_ITERATIONS}}
TOLERANCE = {{TOLERANCE}}
NUM_RUNS = {{NUM_RUNS}}
RANDOM_SEED = {{RANDOM_SEED}}

print(f"Running PyTorch benchmark for {FUNCTION_NAME}")
print(f"Problem dimension: {PROBLEM_DIM}")
print(f"Batch size: {BATCH_SIZE}")
"#
            .to_string(),
            tensorflow_template: r#"# TensorFlow benchmark template
import tensorflow as tf

# Configuration
FUNCTION_NAME = "{{FUNCTION_NAME}}"
PROBLEM_DIM = {{PROBLEM_DIM}}
BATCH_SIZE = {{BATCH_SIZE}}
MAX_ITERATIONS = {{MAX_ITERATIONS}}
TOLERANCE = {{TOLERANCE}}
NUM_RUNS = {{NUM_RUNS}}
RANDOM_SEED = {{RANDOM_SEED}}

print(f"Running TensorFlow benchmark for {FUNCTION_NAME}")
print(f"Problem dimension: {PROBLEM_DIM}")
print(f"Batch size: {BATCH_SIZE}")
"#
            .to_string(),
        }
    }

    fn generate_pytorch_script(
        &self,
        function_name: &str,
        problem_dim: usize,
        batch_size: usize,
        config: &CrossFrameworkConfig,
    ) -> String {
        self.pytorch_template
            .replace("{{FUNCTION_NAME}}", function_name)
            .replace("{{PROBLEM_DIM}}", &problem_dim.to_string())
            .replace("{{BATCH_SIZE}}", &batch_size.to_string())
            .replace("{{MAX_ITERATIONS}}", &config.max_iterations.to_string())
            .replace("{{TOLERANCE}}", &config.tolerance.to_string())
            .replace("{{NUM_RUNS}}", &config.num_runs.to_string())
            .replace("{{RANDOM_SEED}}", &config.random_seed.to_string())
    }

    fn generate_tensorflow_script(
        &self,
        function_name: &str,
        problem_dim: usize,
        batch_size: usize,
        config: &CrossFrameworkConfig,
    ) -> String {
        self.tensorflow_template
            .replace("{{FUNCTION_NAME}}", function_name)
            .replace("{{PROBLEM_DIM}}", &problem_dim.to_string())
            .replace("{{BATCH_SIZE}}", &batch_size.to_string())
            .replace("{{MAX_ITERATIONS}}", &config.max_iterations.to_string())
            .replace("{{TOLERANCE}}", &config.tolerance.to_string())
            .replace("{{NUM_RUNS}}", &config.num_runs.to_string())
            .replace("{{RANDOM_SEED}}", &config.random_seed.to_string())
    }
}

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

    #[test]
    fn test_cross_framework_config() {
        let config = CrossFrameworkConfig::default();
        assert!(config.enable_pytorch);
        assert!(config.enable_tensorflow);
        assert_eq!(config.max_iterations, 1000);
        assert_eq!(config.tolerance, 1e-6);
    }

    #[test]
    fn test_optimizer_identifier() {
        let id = OptimizerIdentifier {
            framework: Framework::SciRS2,
            name: "Adam".to_string(),
            version: Some("0.1.0".to_string()),
        };
        assert_eq!(id.to_string(), "SciRS2-Adam-v0.1.0");
    }

    #[test]
    fn test_precision_enum() {
        let precision = Precision::F64;
        assert!(matches!(precision, Precision::F64));
    }

    #[test]
    fn test_framework_display() {
        assert_eq!(Framework::SciRS2.to_string(), "SciRS2");
        assert_eq!(Framework::PyTorch.to_string(), "PyTorch");
        assert_eq!(Framework::TensorFlow.to_string(), "TensorFlow");
    }

    #[test]
    fn test_python_script_generation() {
        let templates = PythonScriptTemplates::new();
        let config = CrossFrameworkConfig::default();

        let script = templates.generate_pytorch_script("Quadratic", 10, 32, &config);
        assert!(script.contains("10")); // problem_dim
        assert!(script.contains("32")); // batch_size
        assert!(script.contains("1000")); // max_iterations
    }
}