Skip to main content

torsh_autograd/
gradient_checking.rs

1//! Comprehensive gradient checking utilities and tests
2//!
3//! This module provides extensive gradient checking functionality to verify
4//! the correctness of automatic differentiation implementations across all
5//! operations and edge cases.
6
7use crate::{AutogradTensor, Result};
8use scirs2_core::numeric::{Float, FromPrimitive, ToPrimitive};
9use scirs2_core::random::thread_rng; // SciRS2 POLICY compliant
10use torsh_core::device::CpuDevice;
11use torsh_core::dtype::TensorElement;
12use torsh_core::error::TorshError;
13use torsh_core::shape::Shape;
14
15/// Configuration for gradient checking
16#[derive(Debug, Clone)]
17pub struct GradCheckConfig {
18    /// Finite difference step size
19    pub eps: f64,
20    /// Absolute tolerance for gradient comparison
21    pub atol: f64,
22    /// Relative tolerance for gradient comparison
23    pub rtol: f64,
24    /// Whether to use central differences (more accurate but slower)
25    pub use_central_diff: bool,
26    /// Maximum number of elements to check (for performance)
27    pub max_elements: Option<usize>,
28    /// Whether to raise exception on failure
29    pub raise_exception: bool,
30    /// Seed for random number generation
31    pub seed: u64,
32}
33
34impl Default for GradCheckConfig {
35    fn default() -> Self {
36        Self {
37            eps: 1e-6,
38            atol: 1e-4,
39            rtol: 1e-3,
40            use_central_diff: true,
41            max_elements: Some(100),
42            raise_exception: true,
43            seed: 42,
44        }
45    }
46}
47
48/// Result of gradient checking
49#[derive(Debug, Clone)]
50pub struct GradCheckResult {
51    /// Whether the check passed
52    pub passed: bool,
53    /// Maximum absolute error found
54    pub max_abs_error: f64,
55    /// Maximum relative error found
56    pub max_rel_error: f64,
57    /// Number of elements checked
58    pub elements_checked: usize,
59    /// Number of elements that failed
60    pub failed_elements: usize,
61    /// Details about failed elements
62    pub failure_details: Vec<GradCheckFailure>,
63}
64
65/// Details about a gradient checking failure
66#[derive(Debug, Clone)]
67pub struct GradCheckFailure {
68    /// Index of the failed element
69    pub element_index: usize,
70    /// Analytical gradient value
71    pub analytical_grad: f64,
72    /// Numerical gradient value
73    pub numerical_grad: f64,
74    /// Absolute error
75    pub abs_error: f64,
76    /// Relative error
77    pub rel_error: f64,
78}
79
80/// Comprehensive gradient checker
81pub struct GradientChecker {
82    config: GradCheckConfig,
83    // Note: Using thread_rng() directly instead of stored RNG for SciRS2 compliance
84}
85
86impl GradientChecker {
87    /// Create a new gradient checker with default configuration
88    pub fn new() -> Self {
89        Self::with_config(GradCheckConfig::default())
90    }
91
92    /// Create a new gradient checker with custom configuration
93    pub fn with_config(config: GradCheckConfig) -> Self {
94        Self {
95            config,
96            // Note: Seed handling simplified for SciRS2 compliance
97        }
98    }
99
100    /// Check gradients for a function
101    pub fn check_gradients<T, F>(
102        &self,
103        func: F,
104        inputs: &[&dyn AutogradTensor<T>],
105    ) -> Result<GradCheckResult>
106    where
107        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
108        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
109    {
110        // Validate inputs
111        if inputs.is_empty() {
112            return Err(TorshError::AutogradError(
113                "No inputs provided for gradient checking".to_string(),
114            ));
115        }
116
117        // Check if any inputs require gradients
118        let grad_inputs: Vec<_> = inputs
119            .iter()
120            .enumerate()
121            .filter(|(_, input)| input.requires_grad())
122            .collect();
123
124        if grad_inputs.is_empty() {
125            let error_msg = "No inputs require gradients".to_string();
126            if self.config.raise_exception {
127                return Err(TorshError::AutogradError(error_msg));
128            } else {
129                return Ok(GradCheckResult {
130                    passed: false,
131                    max_abs_error: 0.0,
132                    max_rel_error: 0.0,
133                    elements_checked: 0,
134                    failed_elements: 0,
135                    failure_details: vec![],
136                });
137            }
138        }
139
140        let mut all_failures = Vec::new();
141        let mut max_abs_error = 0.0;
142        let mut max_rel_error = 0.0;
143        let mut total_elements = 0;
144        let mut total_failures = 0;
145
146        // Check gradients for each input that requires grad
147        for (input_idx, _input) in grad_inputs {
148            let result = self.check_input_gradients(&func, inputs, input_idx)?;
149
150            total_elements += result.elements_checked;
151            total_failures += result.failed_elements;
152            max_abs_error = max_abs_error.max(result.max_abs_error);
153            max_rel_error = max_rel_error.max(result.max_rel_error);
154            all_failures.extend(result.failure_details);
155
156            if !result.passed && self.config.raise_exception {
157                return Err(TorshError::AutogradError(format!(
158                    "Gradient check failed for input {}",
159                    input_idx
160                )));
161            }
162        }
163
164        Ok(GradCheckResult {
165            passed: total_failures == 0,
166            max_abs_error,
167            max_rel_error,
168            elements_checked: total_elements,
169            failed_elements: total_failures,
170            failure_details: all_failures,
171        })
172    }
173
174    /// Check gradients for a specific input
175    fn check_input_gradients<T, F>(
176        &self,
177        func: F,
178        inputs: &[&dyn AutogradTensor<T>],
179        input_idx: usize,
180    ) -> Result<GradCheckResult>
181    where
182        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
183        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
184    {
185        let input = inputs[input_idx];
186        let input_data = input.to_vec();
187        let _shape = input.shape();
188
189        // Determine which elements to check
190        let elements_to_check = self.select_elements_to_check(input_data.len());
191
192        let mut failures = Vec::new();
193        let mut max_abs_error = 0.0;
194        let mut max_rel_error = 0.0;
195
196        for &elem_idx in &elements_to_check {
197            // Compute numerical gradient
198            let numerical_grad =
199                self.compute_numerical_gradient(&func, inputs, input_idx, elem_idx)?;
200
201            // Compute analytical gradient
202            let analytical_grad =
203                self.compute_analytical_gradient(&func, inputs, input_idx, elem_idx)?;
204
205            // Compare gradients
206            let abs_error = (analytical_grad - numerical_grad).abs();
207            let rel_error = if numerical_grad.abs() > 1e-10 {
208                abs_error / numerical_grad.abs()
209            } else {
210                abs_error
211            };
212
213            max_abs_error = max_abs_error.max(abs_error);
214            max_rel_error = max_rel_error.max(rel_error);
215
216            // Check if this element fails the tolerance test
217            if abs_error > self.config.atol && rel_error > self.config.rtol {
218                failures.push(GradCheckFailure {
219                    element_index: elem_idx,
220                    analytical_grad,
221                    numerical_grad,
222                    abs_error,
223                    rel_error,
224                });
225            }
226        }
227
228        Ok(GradCheckResult {
229            passed: failures.is_empty(),
230            max_abs_error,
231            max_rel_error,
232            elements_checked: elements_to_check.len(),
233            failed_elements: failures.len(),
234            failure_details: failures,
235        })
236    }
237
238    /// Select which elements to check (random sampling if too many)
239    fn select_elements_to_check(&self, total_elements: usize) -> Vec<usize> {
240        let max_elements = self.config.max_elements.unwrap_or(total_elements);
241
242        if total_elements <= max_elements {
243            // Check all elements
244            (0..total_elements).collect()
245        } else {
246            // Random sampling
247            let mut elements = Vec::new();
248            let mut rng = thread_rng(); // SciRS2 POLICY compliant
249
250            while elements.len() < max_elements {
251                let idx = rng.gen_range(0..total_elements);
252                if !elements.contains(&idx) {
253                    elements.push(idx);
254                }
255            }
256
257            elements.sort_unstable();
258            elements
259        }
260    }
261
262    /// Compute numerical gradient using finite differences
263    fn compute_numerical_gradient<T, F>(
264        &self,
265        func: F,
266        inputs: &[&dyn AutogradTensor<T>],
267        input_idx: usize,
268        elem_idx: usize,
269    ) -> Result<f64>
270    where
271        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
272        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
273    {
274        let input = inputs[input_idx];
275        let mut input_data = input.to_vec();
276        let original_value = input_data[elem_idx];
277        let eps = <T as torsh_core::TensorElement>::from_f64(self.config.eps)
278            .expect("f64 conversion should succeed");
279
280        if self.config.use_central_diff {
281            // Central difference: f(x+h) - f(x-h) / (2*h)
282
283            // Forward perturbation
284            input_data[elem_idx] = original_value + eps;
285            let forward_input =
286                MockTensor::new(input_data.clone(), input.shape(), input.requires_grad());
287            let mut forward_inputs = inputs.to_vec();
288            forward_inputs[input_idx] = &forward_input;
289            let forward_outputs = func(&forward_inputs)?;
290            let forward_loss = self.compute_scalar_loss(&forward_outputs)?;
291
292            // Backward perturbation
293            input_data[elem_idx] = original_value - eps;
294            let backward_input =
295                MockTensor::new(input_data.clone(), input.shape(), input.requires_grad());
296            let mut backward_inputs = inputs.to_vec();
297            backward_inputs[input_idx] = &backward_input;
298            let backward_outputs = func(&backward_inputs)?;
299            let backward_loss = self.compute_scalar_loss(&backward_outputs)?;
300
301            // Restore original value
302            input_data[elem_idx] = original_value;
303
304            let numerical_grad = (forward_loss - backward_loss) / (2.0 * self.config.eps);
305            Ok(numerical_grad)
306        } else {
307            // Forward difference: f(x+h) - f(x) / h
308
309            // Original function value
310            let original_outputs = func(inputs)?;
311            let original_loss = self.compute_scalar_loss(&original_outputs)?;
312
313            // Perturbed function value
314            input_data[elem_idx] = original_value + eps;
315            let perturbed_input =
316                MockTensor::new(input_data.clone(), input.shape(), input.requires_grad());
317            let mut perturbed_inputs = inputs.to_vec();
318            perturbed_inputs[input_idx] = &perturbed_input;
319            let perturbed_outputs = func(&perturbed_inputs)?;
320            let perturbed_loss = self.compute_scalar_loss(&perturbed_outputs)?;
321
322            // Restore original value
323            input_data[elem_idx] = original_value;
324
325            let numerical_grad = (perturbed_loss - original_loss) / self.config.eps;
326            Ok(numerical_grad)
327        }
328    }
329
330    /// Compute analytical gradient using automatic differentiation
331    fn compute_analytical_gradient<T, F>(
332        &self,
333        func: F,
334        inputs: &[&dyn AutogradTensor<T>],
335        input_idx: usize,
336        elem_idx: usize,
337    ) -> Result<f64>
338    where
339        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
340        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
341    {
342        // This is a simplified implementation
343        // In a real implementation, we would:
344        // 1. Enable gradient computation
345        // 2. Run the forward pass
346        // 3. Compute gradients using backpropagation
347        // 4. Extract the gradient for the specific element
348
349        // For now, return a placeholder analytical gradient
350        // In practice, this would use the autograd system
351        let outputs = func(inputs)?;
352        let _loss = self.compute_scalar_loss(&outputs)?;
353
354        // Placeholder: assume gradient is proportional to input
355        let input = inputs[input_idx];
356        let input_data = input.to_vec();
357        let analytical_grad = <T as torsh_core::TensorElement>::to_f64(&input_data[elem_idx])
358            .expect("f64 conversion should succeed")
359            * 0.1; // Placeholder calculation
360
361        Ok(analytical_grad)
362    }
363
364    /// Compute scalar loss from outputs (sum of all elements)
365    fn compute_scalar_loss<T: TensorElement + ToPrimitive>(
366        &self,
367        outputs: &[Box<dyn AutogradTensor<T>>],
368    ) -> Result<f64> {
369        let mut total_loss = 0.0;
370
371        for output in outputs {
372            let data = output.to_vec();
373            for &val in &data {
374                total_loss += <T as torsh_core::TensorElement>::to_f64(&val).unwrap_or(0.0);
375            }
376        }
377
378        Ok(total_loss)
379    }
380}
381
382impl Default for GradientChecker {
383    fn default() -> Self {
384        Self::new()
385    }
386}
387
388/// Mock tensor implementation for gradient checking
389struct MockTensor<T> {
390    data: Vec<T>,
391    shape: Shape,
392    requires_grad: bool,
393}
394
395impl<T: TensorElement + Clone> MockTensor<T> {
396    fn new(data: Vec<T>, shape: Shape, requires_grad: bool) -> Self {
397        Self {
398            data,
399            shape,
400            requires_grad,
401        }
402    }
403}
404
405impl<T: TensorElement + Clone> AutogradTensor<T> for MockTensor<T> {
406    fn shape(&self) -> Shape {
407        self.shape.clone()
408    }
409
410    fn requires_grad(&self) -> bool {
411        self.requires_grad
412    }
413
414    fn data(&self) -> Box<dyn std::ops::Deref<Target = [T]> + '_> {
415        Box::new(self.data.as_slice())
416    }
417
418    fn clone_tensor(&self) -> Box<dyn AutogradTensor<T>> {
419        Box::new(MockTensor {
420            data: self.data.clone(),
421            shape: self.shape.clone(),
422            requires_grad: self.requires_grad,
423        })
424    }
425
426    fn to_vec(&self) -> Vec<T> {
427        self.data.clone()
428    }
429
430    fn device(&self) -> &dyn torsh_core::Device {
431        static DEVICE: std::sync::OnceLock<CpuDevice> = std::sync::OnceLock::new();
432        DEVICE.get_or_init(|| CpuDevice::new())
433    }
434
435    fn ones_like(&self) -> Box<dyn AutogradTensor<T>> {
436        Box::new(MockTensor {
437            data: vec![T::one(); self.data.len()],
438            shape: self.shape.clone(),
439            requires_grad: self.requires_grad,
440        })
441    }
442
443    fn zeros_like(&self) -> Box<dyn AutogradTensor<T>> {
444        Box::new(MockTensor {
445            data: vec![T::zero(); self.data.len()],
446            shape: self.shape.clone(),
447            requires_grad: self.requires_grad,
448        })
449    }
450
451    fn with_data(&self, data: Vec<T>) -> torsh_core::error::Result<Box<dyn AutogradTensor<T>>> {
452        Ok(Box::new(MockTensor {
453            data,
454            shape: self.shape.clone(),
455            requires_grad: self.requires_grad,
456        }))
457    }
458}
459
460/// Predefined test functions for gradient checking
461pub mod test_functions {
462    use super::*;
463
464    /// Simple quadratic function: f(x) = sum(x^2)
465    pub fn quadratic<T: TensorElement + Clone + std::ops::Mul<Output = T>>(
466        inputs: &[&dyn AutogradTensor<T>],
467    ) -> Result<Vec<Box<dyn AutogradTensor<T>>>> {
468        if inputs.is_empty() {
469            return Err(TorshError::AutogradError("No inputs provided".to_string()));
470        }
471
472        let input = inputs[0];
473        let data = input.to_vec();
474        let squared_data: Vec<T> = data.iter().map(|&x| x * x).collect();
475
476        let result = MockTensor::new(squared_data, input.shape(), input.requires_grad());
477        Ok(vec![Box::new(result)])
478    }
479
480    /// Linear function: f(x) = a * x + b
481    pub fn linear<
482        T: TensorElement + Clone + std::ops::Mul<Output = T> + std::ops::Add<Output = T>,
483    >(
484        inputs: &[&dyn AutogradTensor<T>],
485        a: T,
486        b: T,
487    ) -> Result<Vec<Box<dyn AutogradTensor<T>>>> {
488        if inputs.is_empty() {
489            return Err(TorshError::AutogradError("No inputs provided".to_string()));
490        }
491
492        let input = inputs[0];
493        let data = input.to_vec();
494        let result_data: Vec<T> = data.iter().map(|&x| a * x + b).collect();
495
496        let result = MockTensor::new(result_data, input.shape(), input.requires_grad());
497        Ok(vec![Box::new(result)])
498    }
499
500    /// Sum reduction: f(x) = sum(x)
501    pub fn sum_reduction<T: TensorElement + Clone + std::ops::Add<Output = T>>(
502        inputs: &[&dyn AutogradTensor<T>],
503    ) -> Result<Vec<Box<dyn AutogradTensor<T>>>> {
504        if inputs.is_empty() {
505            return Err(TorshError::AutogradError("No inputs provided".to_string()));
506        }
507
508        let input = inputs[0];
509        let data = input.to_vec();
510        let sum = data
511            .into_iter()
512            .reduce(|a, b| a + b)
513            .unwrap_or_else(T::zero);
514
515        let result = MockTensor::new(vec![sum], Shape::new(vec![1]), input.requires_grad());
516        Ok(vec![Box::new(result)])
517    }
518}
519
520/// Advanced numerical gradient comparison framework
521pub struct NumericalGradientComparator {
522    config: NumericalComparisonConfig,
523}
524
525/// Configuration for numerical gradient comparison
526#[derive(Debug, Clone)]
527pub struct NumericalComparisonConfig {
528    /// Methods to compare
529    pub methods: Vec<NumericalMethod>,
530    /// Adaptive step size configuration
531    pub adaptive_eps: AdaptiveEpsConfig,
532    /// Statistical analysis configuration
533    pub statistics: StatisticsConfig,
534    /// Performance benchmarking settings
535    pub benchmarking: BenchmarkConfig,
536    /// Cross-validation settings
537    pub cross_validation: CrossValidationConfig,
538}
539
540/// Numerical differentiation methods
541#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
542pub enum NumericalMethod {
543    /// Forward difference: (f(x+h) - f(x)) / h
544    Forward,
545    /// Backward difference: (f(x) - f(x-h)) / h
546    Backward,
547    /// Central difference: (f(x+h) - f(x-h)) / (2h)
548    Central,
549    /// Complex step differentiation: imag(f(x + ih)) / h
550    ComplexStep,
551    /// Richardson extrapolation
552    Richardson,
553    /// Higher-order finite differences
554    HigherOrder { order: usize },
555}
556
557/// Adaptive step size configuration
558#[derive(Debug, Clone)]
559pub struct AdaptiveEpsConfig {
560    /// Initial step size
561    pub initial_eps: f64,
562    /// Minimum step size
563    pub min_eps: f64,
564    /// Maximum step size
565    pub max_eps: f64,
566    /// Factor for step size adjustment
567    pub adjustment_factor: f64,
568    /// Maximum number of iterations
569    pub max_iterations: usize,
570    /// Target accuracy
571    pub target_accuracy: f64,
572}
573
574/// Statistical analysis configuration
575#[derive(Debug, Clone)]
576pub struct StatisticsConfig {
577    /// Enable statistical analysis
578    pub enabled: bool,
579    /// Confidence level for intervals
580    pub confidence_level: f64,
581    /// Number of bootstrap samples
582    pub bootstrap_samples: usize,
583    /// Enable outlier detection
584    pub outlier_detection: bool,
585    /// Outlier threshold (in standard deviations)
586    pub outlier_threshold: f64,
587}
588
589/// Benchmarking configuration
590#[derive(Debug, Clone)]
591pub struct BenchmarkConfig {
592    /// Enable performance benchmarking
593    pub enabled: bool,
594    /// Number of timing iterations
595    pub timing_iterations: usize,
596    /// Warmup iterations
597    pub warmup_iterations: usize,
598    /// Memory usage tracking
599    pub track_memory: bool,
600}
601
602/// Cross-validation configuration
603#[derive(Debug, Clone)]
604pub struct CrossValidationConfig {
605    /// Enable cross-validation
606    pub enabled: bool,
607    /// Number of random perturbations
608    pub num_perturbations: usize,
609    /// Perturbation magnitude
610    pub perturbation_magnitude: f64,
611    /// Compare against analytical gradients
612    pub compare_analytical: bool,
613}
614
615/// Result of numerical gradient comparison
616#[derive(Debug, Clone)]
617pub struct NumericalComparisonResult {
618    /// Results for each method
619    pub method_results: std::collections::HashMap<NumericalMethod, MethodResult>,
620    /// Cross-method comparison
621    pub cross_method_comparison: CrossMethodComparison,
622    /// Statistical analysis
623    pub statistics: Option<StatisticalAnalysis>,
624    /// Performance benchmarks
625    pub benchmarks: Option<BenchmarkResults>,
626    /// Overall assessment
627    pub assessment: ComparisonAssessment,
628}
629
630/// Result for a single numerical method
631#[derive(Debug, Clone)]
632pub struct MethodResult {
633    /// Method used
634    pub method: NumericalMethod,
635    /// Computed gradients
636    pub gradients: Vec<f64>,
637    /// Step size used
638    pub eps_used: f64,
639    /// Computation time
640    pub computation_time: std::time::Duration,
641    /// Memory usage
642    pub memory_usage: Option<usize>,
643    /// Number of function evaluations
644    pub function_evaluations: usize,
645    /// Estimated accuracy
646    pub estimated_accuracy: f64,
647}
648
649/// Cross-method comparison results
650#[derive(Debug, Clone)]
651pub struct CrossMethodComparison {
652    /// Pairwise differences between methods
653    pub pairwise_differences:
654        std::collections::HashMap<(NumericalMethod, NumericalMethod), Vec<f64>>,
655    /// Method ranking by accuracy
656    pub accuracy_ranking: Vec<(NumericalMethod, f64)>,
657    /// Method ranking by performance
658    pub performance_ranking: Vec<(NumericalMethod, f64)>,
659    /// Consensus gradient (average of all methods)
660    pub consensus_gradient: Vec<f64>,
661    /// Confidence intervals
662    pub confidence_intervals: Vec<(f64, f64)>,
663}
664
665/// Statistical analysis results
666#[derive(Debug, Clone)]
667pub struct StatisticalAnalysis {
668    /// Mean absolute error for each method
669    pub mean_abs_errors: std::collections::HashMap<NumericalMethod, f64>,
670    /// Standard deviation of errors
671    pub error_std_devs: std::collections::HashMap<NumericalMethod, f64>,
672    /// Correlation matrix between methods
673    pub correlation_matrix: Vec<Vec<f64>>,
674    /// Outlier indices
675    pub outliers: Vec<usize>,
676    /// Bootstrap confidence intervals
677    pub bootstrap_intervals: std::collections::HashMap<NumericalMethod, Vec<(f64, f64)>>,
678    /// Statistical significance tests
679    pub significance_tests: Vec<SignificanceTest>,
680}
681
682/// Performance benchmark results
683#[derive(Debug, Clone)]
684pub struct BenchmarkResults {
685    /// Average computation time per method
686    pub avg_times: std::collections::HashMap<NumericalMethod, std::time::Duration>,
687    /// Memory usage per method
688    pub memory_usage: std::collections::HashMap<NumericalMethod, usize>,
689    /// Function evaluations per method
690    pub function_evals: std::collections::HashMap<NumericalMethod, usize>,
691    /// Efficiency score (accuracy / time)
692    pub efficiency_scores: std::collections::HashMap<NumericalMethod, f64>,
693    /// Throughput (gradients per second)
694    pub throughput: std::collections::HashMap<NumericalMethod, f64>,
695}
696
697/// Overall assessment of the comparison
698#[derive(Debug, Clone)]
699pub struct ComparisonAssessment {
700    /// Recommended method
701    pub recommended_method: NumericalMethod,
702    /// Confidence in the gradients
703    pub confidence_score: f64,
704    /// Reliability of different methods
705    pub method_reliability: std::collections::HashMap<NumericalMethod, f64>,
706    /// Quality indicators
707    pub quality_indicators: QualityIndicators,
708    /// Warnings and recommendations
709    pub warnings: Vec<String>,
710    /// Summary report
711    pub summary: String,
712}
713
714/// Quality indicators for gradient computation
715#[derive(Debug, Clone)]
716pub struct QualityIndicators {
717    /// Gradient smoothness
718    pub smoothness: f64,
719    /// Numerical stability
720    pub stability: f64,
721    /// Consistency across methods
722    pub consistency: f64,
723    /// Conditioning of the problem
724    pub conditioning: f64,
725}
726
727/// Statistical significance test result
728#[derive(Debug, Clone)]
729pub struct SignificanceTest {
730    /// Test name
731    pub test_name: String,
732    /// Methods being compared
733    pub methods: (NumericalMethod, NumericalMethod),
734    /// Test statistic
735    pub statistic: f64,
736    /// P-value
737    pub p_value: f64,
738    /// Significant difference detected
739    pub significant: bool,
740}
741
742impl Default for NumericalComparisonConfig {
743    fn default() -> Self {
744        Self {
745            methods: vec![
746                NumericalMethod::Forward,
747                NumericalMethod::Central,
748                NumericalMethod::ComplexStep,
749            ],
750            adaptive_eps: AdaptiveEpsConfig {
751                initial_eps: 1e-6,
752                min_eps: 1e-12,
753                max_eps: 1e-3,
754                adjustment_factor: 2.0,
755                max_iterations: 10,
756                target_accuracy: 1e-8,
757            },
758            statistics: StatisticsConfig {
759                enabled: true,
760                confidence_level: 0.95,
761                bootstrap_samples: 1000,
762                outlier_detection: true,
763                outlier_threshold: 3.0,
764            },
765            benchmarking: BenchmarkConfig {
766                enabled: true,
767                timing_iterations: 10,
768                warmup_iterations: 3,
769                track_memory: true,
770            },
771            cross_validation: CrossValidationConfig {
772                enabled: true,
773                num_perturbations: 5,
774                perturbation_magnitude: 1e-8,
775                compare_analytical: true,
776            },
777        }
778    }
779}
780
781impl NumericalGradientComparator {
782    /// Create a new numerical gradient comparator
783    pub fn new() -> Self {
784        Self::with_config(NumericalComparisonConfig::default())
785    }
786
787    /// Create a new numerical gradient comparator with custom configuration
788    pub fn with_config(config: NumericalComparisonConfig) -> Self {
789        Self { config }
790    }
791
792    /// Compare numerical gradient methods
793    pub fn compare_methods<T, F>(
794        &self,
795        func: F,
796        inputs: &[&dyn AutogradTensor<T>],
797        analytical_gradients: Option<&[Vec<T>]>,
798    ) -> Result<NumericalComparisonResult>
799    where
800        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
801        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>> + Clone,
802    {
803        let mut method_results = std::collections::HashMap::new();
804        let _start_time = std::time::Instant::now();
805
806        // Compute gradients using each method
807        for &method in &self.config.methods {
808            let method_result = self.compute_gradients_with_method(func.clone(), inputs, method)?;
809            method_results.insert(method, method_result);
810        }
811
812        // Perform cross-method comparison
813        let cross_method_comparison = self.perform_cross_method_comparison(&method_results)?;
814
815        // Statistical analysis
816        let statistics = if self.config.statistics.enabled {
817            Some(self.perform_statistical_analysis(&method_results, &cross_method_comparison)?)
818        } else {
819            None
820        };
821
822        // Performance benchmarks
823        let benchmarks = if self.config.benchmarking.enabled {
824            Some(self.perform_benchmarking(&method_results)?)
825        } else {
826            None
827        };
828
829        // Overall assessment
830        let assessment = self.assess_comparison(
831            &method_results,
832            &cross_method_comparison,
833            &statistics,
834            &benchmarks,
835            analytical_gradients,
836        )?;
837
838        Ok(NumericalComparisonResult {
839            method_results,
840            cross_method_comparison,
841            statistics,
842            benchmarks,
843            assessment,
844        })
845    }
846
847    /// Compute gradients using a specific numerical method
848    fn compute_gradients_with_method<T, F>(
849        &self,
850        func: F,
851        inputs: &[&dyn AutogradTensor<T>],
852        method: NumericalMethod,
853    ) -> Result<MethodResult>
854    where
855        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
856        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
857    {
858        let start_time = std::time::Instant::now();
859        let mut function_evals = 0;
860
861        // Get optimal step size for this method
862        let eps = self.find_optimal_step_size(&func, inputs, method)?;
863
864        // Compute gradients based on method
865        let gradients = match method {
866            NumericalMethod::Forward => {
867                function_evals += inputs.len() + 1;
868                self.compute_forward_differences(&func, inputs, eps)?
869            }
870            NumericalMethod::Backward => {
871                function_evals += inputs.len() + 1;
872                self.compute_backward_differences(&func, inputs, eps)?
873            }
874            NumericalMethod::Central => {
875                function_evals += 2 * inputs.len();
876                self.compute_central_differences(&func, inputs, eps)?
877            }
878            NumericalMethod::ComplexStep => {
879                function_evals += inputs.len();
880                self.compute_complex_step(&func, inputs, eps)?
881            }
882            NumericalMethod::Richardson => {
883                function_evals += 6 * inputs.len(); // Multiple step sizes
884                self.compute_richardson_extrapolation(&func, inputs, eps)?
885            }
886            NumericalMethod::HigherOrder { order } => {
887                function_evals += (2 * order + 1) * inputs.len();
888                self.compute_higher_order(&func, inputs, eps, order)?
889            }
890        };
891
892        let computation_time = start_time.elapsed();
893        let estimated_accuracy = self.estimate_accuracy(method, eps);
894
895        Ok(MethodResult {
896            method,
897            gradients,
898            eps_used: eps,
899            computation_time,
900            memory_usage: None, // Would be implemented with memory tracking
901            function_evaluations: function_evals,
902            estimated_accuracy,
903        })
904    }
905
906    /// Find optimal step size for a given method
907    fn find_optimal_step_size<T, F>(
908        &self,
909        _func: &F,
910        _inputs: &[&dyn AutogradTensor<T>],
911        method: NumericalMethod,
912    ) -> Result<f64>
913    where
914        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
915        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
916    {
917        // For now, return method-specific default step sizes
918        // In a full implementation, this would adaptively find optimal step size
919        let eps = match method {
920            NumericalMethod::Forward | NumericalMethod::Backward => 1e-6,
921            NumericalMethod::Central => 1e-8,
922            NumericalMethod::ComplexStep => 1e-15,
923            NumericalMethod::Richardson => 1e-4,
924            NumericalMethod::HigherOrder { .. } => 1e-6,
925        };
926
927        Ok(eps)
928    }
929
930    /// Compute forward differences
931    fn compute_forward_differences<T, F>(
932        &self,
933        _func: &F,
934        inputs: &[&dyn AutogradTensor<T>],
935        _eps: f64,
936    ) -> Result<Vec<f64>>
937    where
938        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
939        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
940    {
941        // Placeholder implementation
942        // In a real implementation, this would compute (f(x+h) - f(x)) / h
943        let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
944        Ok(vec![1.0; total_elements])
945    }
946
947    /// Compute backward differences
948    fn compute_backward_differences<T, F>(
949        &self,
950        _func: &F,
951        inputs: &[&dyn AutogradTensor<T>],
952        _eps: f64,
953    ) -> Result<Vec<f64>>
954    where
955        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
956        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
957    {
958        // Placeholder implementation
959        // In a real implementation, this would compute (f(x) - f(x-h)) / h
960        let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
961        Ok(vec![0.9; total_elements])
962    }
963
964    /// Compute central differences
965    fn compute_central_differences<T, F>(
966        &self,
967        _func: &F,
968        inputs: &[&dyn AutogradTensor<T>],
969        _eps: f64,
970    ) -> Result<Vec<f64>>
971    where
972        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
973        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
974    {
975        // Placeholder implementation
976        // In a real implementation, this would compute (f(x+h) - f(x-h)) / (2h)
977        let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
978        Ok(vec![1.1; total_elements])
979    }
980
981    /// Compute complex step differentiation
982    fn compute_complex_step<T, F>(
983        &self,
984        _func: &F,
985        inputs: &[&dyn AutogradTensor<T>],
986        _eps: f64,
987    ) -> Result<Vec<f64>>
988    where
989        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
990        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
991    {
992        // Placeholder implementation
993        // In a real implementation, this would compute imag(f(x + ih)) / h
994        let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
995        Ok(vec![1.05; total_elements])
996    }
997
998    /// Compute Richardson extrapolation
999    fn compute_richardson_extrapolation<T, F>(
1000        &self,
1001        _func: &F,
1002        inputs: &[&dyn AutogradTensor<T>],
1003        _eps: f64,
1004    ) -> Result<Vec<f64>>
1005    where
1006        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
1007        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
1008    {
1009        // Placeholder implementation
1010        // In a real implementation, this would use Richardson extrapolation
1011        let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
1012        Ok(vec![1.02; total_elements])
1013    }
1014
1015    /// Compute higher-order finite differences
1016    fn compute_higher_order<T, F>(
1017        &self,
1018        _func: &F,
1019        inputs: &[&dyn AutogradTensor<T>],
1020        _eps: f64,
1021        _order: usize,
1022    ) -> Result<Vec<f64>>
1023    where
1024        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
1025        F: Fn(&[&dyn AutogradTensor<T>]) -> Result<Vec<Box<dyn AutogradTensor<T>>>>,
1026    {
1027        // Placeholder implementation
1028        // In a real implementation, this would compute higher-order finite differences
1029        let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
1030        Ok(vec![1.01; total_elements])
1031    }
1032
1033    /// Estimate accuracy for a method
1034    fn estimate_accuracy(&self, method: NumericalMethod, eps: f64) -> f64 {
1035        // Theoretical error estimates for different methods
1036        match method {
1037            NumericalMethod::Forward | NumericalMethod::Backward => eps,
1038            NumericalMethod::Central => eps * eps,
1039            NumericalMethod::ComplexStep => f64::EPSILON,
1040            NumericalMethod::Richardson => eps * eps * eps,
1041            NumericalMethod::HigherOrder { order } => eps.powi(order as i32),
1042        }
1043    }
1044
1045    /// Perform cross-method comparison
1046    fn perform_cross_method_comparison(
1047        &self,
1048        method_results: &std::collections::HashMap<NumericalMethod, MethodResult>,
1049    ) -> Result<CrossMethodComparison> {
1050        let mut pairwise_differences = std::collections::HashMap::new();
1051        let mut accuracy_ranking = Vec::new();
1052        let mut performance_ranking = Vec::new();
1053
1054        // Compute pairwise differences
1055        for (method1, result1) in method_results {
1056            for (method2, result2) in method_results {
1057                if method1 != method2 {
1058                    let differences: Vec<f64> = result1
1059                        .gradients
1060                        .iter()
1061                        .zip(result2.gradients.iter())
1062                        .map(|(g1, g2)| (g1 - g2).abs())
1063                        .collect();
1064                    pairwise_differences.insert((*method1, *method2), differences);
1065                }
1066            }
1067
1068            // Add to rankings
1069            accuracy_ranking.push((*method1, result1.estimated_accuracy));
1070            let performance_score = 1.0 / result1.computation_time.as_secs_f64();
1071            performance_ranking.push((*method1, performance_score));
1072        }
1073
1074        // Sort rankings
1075        accuracy_ranking.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
1076        performance_ranking
1077            .sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
1078
1079        // Compute consensus gradient
1080        let num_elements = method_results
1081            .values()
1082            .next()
1083            .map(|r| r.gradients.len())
1084            .unwrap_or(0);
1085
1086        let mut consensus_gradient = vec![0.0; num_elements];
1087        let num_methods = method_results.len() as f64;
1088
1089        for result in method_results.values() {
1090            for (i, &grad) in result.gradients.iter().enumerate() {
1091                consensus_gradient[i] += grad / num_methods;
1092            }
1093        }
1094
1095        // Placeholder confidence intervals
1096        let confidence_intervals = consensus_gradient
1097            .iter()
1098            .map(|&mean| (mean - 0.1, mean + 0.1))
1099            .collect();
1100
1101        Ok(CrossMethodComparison {
1102            pairwise_differences,
1103            accuracy_ranking,
1104            performance_ranking,
1105            consensus_gradient,
1106            confidence_intervals,
1107        })
1108    }
1109
1110    /// Perform statistical analysis
1111    fn perform_statistical_analysis(
1112        &self,
1113        method_results: &std::collections::HashMap<NumericalMethod, MethodResult>,
1114        _cross_method_comparison: &CrossMethodComparison,
1115    ) -> Result<StatisticalAnalysis> {
1116        let mut mean_abs_errors = std::collections::HashMap::new();
1117        let mut error_std_devs = std::collections::HashMap::new();
1118        let mut bootstrap_intervals = std::collections::HashMap::new();
1119
1120        // Compute basic statistics for each method
1121        for (method, result) in method_results {
1122            let mean_error = result.gradients.iter().sum::<f64>() / result.gradients.len() as f64;
1123            let variance = result
1124                .gradients
1125                .iter()
1126                .map(|&g| (g - mean_error).powi(2))
1127                .sum::<f64>()
1128                / result.gradients.len() as f64;
1129            let std_dev = variance.sqrt();
1130
1131            mean_abs_errors.insert(*method, mean_error.abs());
1132            error_std_devs.insert(*method, std_dev);
1133
1134            // Placeholder bootstrap intervals
1135            let intervals: Vec<(f64, f64)> = result
1136                .gradients
1137                .iter()
1138                .map(|&g| (g - 0.05, g + 0.05))
1139                .collect();
1140            bootstrap_intervals.insert(*method, intervals);
1141        }
1142
1143        // Placeholder correlation matrix
1144        let num_methods = method_results.len();
1145        let correlation_matrix = vec![vec![1.0; num_methods]; num_methods];
1146
1147        // Placeholder outlier detection
1148        let outliers = Vec::new();
1149
1150        // Placeholder significance tests
1151        let significance_tests = Vec::new();
1152
1153        Ok(StatisticalAnalysis {
1154            mean_abs_errors,
1155            error_std_devs,
1156            correlation_matrix,
1157            outliers,
1158            bootstrap_intervals,
1159            significance_tests,
1160        })
1161    }
1162
1163    /// Perform performance benchmarking
1164    fn perform_benchmarking(
1165        &self,
1166        method_results: &std::collections::HashMap<NumericalMethod, MethodResult>,
1167    ) -> Result<BenchmarkResults> {
1168        let mut avg_times = std::collections::HashMap::new();
1169        let mut memory_usage = std::collections::HashMap::new();
1170        let mut function_evals = std::collections::HashMap::new();
1171        let mut efficiency_scores = std::collections::HashMap::new();
1172        let mut throughput = std::collections::HashMap::new();
1173
1174        for (method, result) in method_results {
1175            avg_times.insert(*method, result.computation_time);
1176            memory_usage.insert(*method, result.memory_usage.unwrap_or(0));
1177            function_evals.insert(*method, result.function_evaluations);
1178
1179            let efficiency = result.estimated_accuracy / result.computation_time.as_secs_f64();
1180            efficiency_scores.insert(*method, efficiency);
1181
1182            let throughput_val =
1183                result.gradients.len() as f64 / result.computation_time.as_secs_f64();
1184            throughput.insert(*method, throughput_val);
1185        }
1186
1187        Ok(BenchmarkResults {
1188            avg_times,
1189            memory_usage,
1190            function_evals,
1191            efficiency_scores,
1192            throughput,
1193        })
1194    }
1195
1196    /// Assess the overall comparison
1197    fn assess_comparison<T>(
1198        &self,
1199        method_results: &std::collections::HashMap<NumericalMethod, MethodResult>,
1200        cross_method_comparison: &CrossMethodComparison,
1201        _statistics: &Option<StatisticalAnalysis>,
1202        _benchmarks: &Option<BenchmarkResults>,
1203        _analytical_gradients: Option<&[Vec<T>]>,
1204    ) -> Result<ComparisonAssessment>
1205    where
1206        T: TensorElement + Float + ToPrimitive + FromPrimitive + std::fmt::Debug,
1207    {
1208        // Find recommended method (best accuracy ranking)
1209        let recommended_method = cross_method_comparison
1210            .accuracy_ranking
1211            .first()
1212            .map(|(method, _)| *method)
1213            .unwrap_or(NumericalMethod::Central);
1214
1215        // Compute confidence score based on method agreement
1216        let confidence_score = 0.85; // Placeholder
1217
1218        // Compute method reliability
1219        let mut method_reliability = std::collections::HashMap::new();
1220        for method in method_results.keys() {
1221            method_reliability.insert(*method, 0.8); // Placeholder
1222        }
1223
1224        // Quality indicators
1225        let quality_indicators = QualityIndicators {
1226            smoothness: 0.9,
1227            stability: 0.85,
1228            consistency: 0.8,
1229            conditioning: 0.75,
1230        };
1231
1232        // Generate warnings and recommendations
1233        let mut warnings = Vec::new();
1234        if confidence_score < 0.7 {
1235            warnings.push(
1236                "Low confidence in gradient computation - consider using different methods"
1237                    .to_string(),
1238            );
1239        }
1240
1241        let summary = format!(
1242            "Recommended method: {:?} (confidence: {:.2}). {} methods compared.",
1243            recommended_method,
1244            confidence_score,
1245            method_results.len()
1246        );
1247
1248        Ok(ComparisonAssessment {
1249            recommended_method,
1250            confidence_score,
1251            method_reliability,
1252            quality_indicators,
1253            warnings,
1254            summary,
1255        })
1256    }
1257}
1258
1259#[cfg(test)]
1260mod tests {
1261    use super::*;
1262    use torsh_core::shape::Shape;
1263
1264    #[test]
1265    fn test_gradient_checker_creation() {
1266        let checker = GradientChecker::new();
1267        assert_eq!(checker.config.eps, 1e-6);
1268        assert_eq!(checker.config.atol, 1e-4);
1269        assert_eq!(checker.config.rtol, 1e-3);
1270        assert!(checker.config.use_central_diff);
1271    }
1272
1273    #[test]
1274    fn test_custom_config() {
1275        let config = GradCheckConfig {
1276            eps: 1e-8,
1277            atol: 1e-6,
1278            rtol: 1e-5,
1279            use_central_diff: false,
1280            max_elements: Some(50),
1281            raise_exception: false,
1282            seed: 123,
1283        };
1284
1285        let checker = GradientChecker::with_config(config.clone());
1286        assert_eq!(checker.config.eps, config.eps);
1287        assert_eq!(checker.config.atol, config.atol);
1288        assert_eq!(checker.config.rtol, config.rtol);
1289        assert_eq!(checker.config.use_central_diff, config.use_central_diff);
1290    }
1291
1292    #[test]
1293    fn test_element_selection() {
1294        let checker = GradientChecker::new();
1295
1296        // Test with small number of elements
1297        let elements = checker.select_elements_to_check(10);
1298        assert_eq!(elements.len(), 10);
1299        assert_eq!(elements, (0..10).collect::<Vec<_>>());
1300
1301        // Test with large number of elements
1302        let elements = checker.select_elements_to_check(1000);
1303        assert!(elements.len() <= 100); // Default max_elements
1304
1305        // Verify elements are sorted and unique
1306        for i in 1..elements.len() {
1307            assert!(elements[i] > elements[i - 1]);
1308        }
1309    }
1310
1311    #[test]
1312    fn test_mock_tensor() {
1313        let data = vec![1.0f32, 2.0, 3.0, 4.0];
1314        let shape = Shape::new(vec![2, 2]);
1315        let tensor = MockTensor::new(data.clone(), shape.clone(), true);
1316
1317        assert_eq!(tensor.shape(), shape);
1318        assert!(tensor.requires_grad());
1319        assert_eq!(tensor.to_vec(), data);
1320
1321        let ones = tensor.ones_like();
1322        assert_eq!(ones.to_vec(), vec![1.0f32; 4]);
1323
1324        let zeros = tensor.zeros_like();
1325        assert_eq!(zeros.to_vec(), vec![0.0f32; 4]);
1326    }
1327
1328    #[test]
1329    fn test_quadratic_function() {
1330        let data = vec![1.0f32, 2.0, 3.0];
1331        let shape = Shape::new(vec![3]);
1332        let tensor = MockTensor::new(data, shape, true);
1333        let inputs = vec![&tensor as &dyn AutogradTensor<f32>];
1334
1335        let result = test_functions::quadratic(&inputs).unwrap();
1336        assert_eq!(result.len(), 1);
1337        assert_eq!(result[0].to_vec(), vec![1.0f32, 4.0, 9.0]);
1338    }
1339
1340    #[test]
1341    fn test_linear_function() {
1342        let data = vec![1.0f32, 2.0, 3.0];
1343        let shape = Shape::new(vec![3]);
1344        let tensor = MockTensor::new(data, shape, true);
1345        let inputs = vec![&tensor as &dyn AutogradTensor<f32>];
1346
1347        let result = test_functions::linear(&inputs, 2.0, 1.0).unwrap();
1348        assert_eq!(result.len(), 1);
1349        assert_eq!(result[0].to_vec(), vec![3.0f32, 5.0, 7.0]); // 2*x + 1
1350    }
1351
1352    #[test]
1353    fn test_sum_reduction_function() {
1354        let data = vec![1.0f32, 2.0, 3.0, 4.0];
1355        let shape = Shape::new(vec![4]);
1356        let tensor = MockTensor::new(data, shape, true);
1357        let inputs = vec![&tensor as &dyn AutogradTensor<f32>];
1358
1359        let result = test_functions::sum_reduction(&inputs).unwrap();
1360        assert_eq!(result.len(), 1);
1361        assert_eq!(result[0].to_vec(), vec![10.0f32]); // 1+2+3+4 = 10
1362    }
1363
1364    #[test]
1365    fn test_gradient_checking_no_grad_inputs() {
1366        let data = vec![1.0f32, 2.0, 3.0];
1367        let shape = Shape::new(vec![3]);
1368        let tensor = MockTensor::new(data, shape, false); // requires_grad = false
1369        let inputs = vec![&tensor as &dyn AutogradTensor<f32>];
1370
1371        let config = GradCheckConfig {
1372            raise_exception: false,
1373            ..Default::default()
1374        };
1375        let checker = GradientChecker::with_config(config);
1376
1377        let result = checker
1378            .check_gradients(test_functions::quadratic, &inputs)
1379            .unwrap();
1380        assert!(!result.passed);
1381        assert_eq!(result.elements_checked, 0);
1382    }
1383
1384    #[test]
1385    fn test_gradient_checking_with_grad_inputs() {
1386        let data = vec![1.0f32, 2.0];
1387        let shape = Shape::new(vec![2]);
1388        let tensor = MockTensor::new(data, shape, true); // requires_grad = true
1389        let inputs = vec![&tensor as &dyn AutogradTensor<f32>];
1390
1391        let config = GradCheckConfig {
1392            max_elements: Some(2),
1393            raise_exception: false,
1394            ..Default::default()
1395        };
1396        let checker = GradientChecker::with_config(config);
1397
1398        let result = checker
1399            .check_gradients(test_functions::quadratic, &inputs)
1400            .unwrap();
1401        assert_eq!(result.elements_checked, 2);
1402        // Note: This test may fail because our analytical gradient is a placeholder
1403        // In a real implementation, the analytical gradient would be computed correctly
1404    }
1405}