1use crate::{AutogradTensor, Result};
8use scirs2_core::numeric::{Float, FromPrimitive, ToPrimitive};
9use scirs2_core::random::thread_rng; use torsh_core::device::CpuDevice;
11use torsh_core::dtype::TensorElement;
12use torsh_core::error::TorshError;
13use torsh_core::shape::Shape;
14
15#[derive(Debug, Clone)]
17pub struct GradCheckConfig {
18 pub eps: f64,
20 pub atol: f64,
22 pub rtol: f64,
24 pub use_central_diff: bool,
26 pub max_elements: Option<usize>,
28 pub raise_exception: bool,
30 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#[derive(Debug, Clone)]
50pub struct GradCheckResult {
51 pub passed: bool,
53 pub max_abs_error: f64,
55 pub max_rel_error: f64,
57 pub elements_checked: usize,
59 pub failed_elements: usize,
61 pub failure_details: Vec<GradCheckFailure>,
63}
64
65#[derive(Debug, Clone)]
67pub struct GradCheckFailure {
68 pub element_index: usize,
70 pub analytical_grad: f64,
72 pub numerical_grad: f64,
74 pub abs_error: f64,
76 pub rel_error: f64,
78}
79
80pub struct GradientChecker {
82 config: GradCheckConfig,
83 }
85
86impl GradientChecker {
87 pub fn new() -> Self {
89 Self::with_config(GradCheckConfig::default())
90 }
91
92 pub fn with_config(config: GradCheckConfig) -> Self {
94 Self {
95 config,
96 }
98 }
99
100 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 if inputs.is_empty() {
112 return Err(TorshError::AutogradError(
113 "No inputs provided for gradient checking".to_string(),
114 ));
115 }
116
117 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 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 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 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 let numerical_grad =
199 self.compute_numerical_gradient(&func, inputs, input_idx, elem_idx)?;
200
201 let analytical_grad =
203 self.compute_analytical_gradient(&func, inputs, input_idx, elem_idx)?;
204
205 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 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 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 (0..total_elements).collect()
245 } else {
246 let mut elements = Vec::new();
248 let mut rng = thread_rng(); 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 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 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 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 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 let original_outputs = func(inputs)?;
311 let original_loss = self.compute_scalar_loss(&original_outputs)?;
312
313 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 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 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 let outputs = func(inputs)?;
352 let _loss = self.compute_scalar_loss(&outputs)?;
353
354 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; Ok(analytical_grad)
362 }
363
364 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
388struct 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
460pub mod test_functions {
462 use super::*;
463
464 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 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 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
520pub struct NumericalGradientComparator {
522 config: NumericalComparisonConfig,
523}
524
525#[derive(Debug, Clone)]
527pub struct NumericalComparisonConfig {
528 pub methods: Vec<NumericalMethod>,
530 pub adaptive_eps: AdaptiveEpsConfig,
532 pub statistics: StatisticsConfig,
534 pub benchmarking: BenchmarkConfig,
536 pub cross_validation: CrossValidationConfig,
538}
539
540#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
542pub enum NumericalMethod {
543 Forward,
545 Backward,
547 Central,
549 ComplexStep,
551 Richardson,
553 HigherOrder { order: usize },
555}
556
557#[derive(Debug, Clone)]
559pub struct AdaptiveEpsConfig {
560 pub initial_eps: f64,
562 pub min_eps: f64,
564 pub max_eps: f64,
566 pub adjustment_factor: f64,
568 pub max_iterations: usize,
570 pub target_accuracy: f64,
572}
573
574#[derive(Debug, Clone)]
576pub struct StatisticsConfig {
577 pub enabled: bool,
579 pub confidence_level: f64,
581 pub bootstrap_samples: usize,
583 pub outlier_detection: bool,
585 pub outlier_threshold: f64,
587}
588
589#[derive(Debug, Clone)]
591pub struct BenchmarkConfig {
592 pub enabled: bool,
594 pub timing_iterations: usize,
596 pub warmup_iterations: usize,
598 pub track_memory: bool,
600}
601
602#[derive(Debug, Clone)]
604pub struct CrossValidationConfig {
605 pub enabled: bool,
607 pub num_perturbations: usize,
609 pub perturbation_magnitude: f64,
611 pub compare_analytical: bool,
613}
614
615#[derive(Debug, Clone)]
617pub struct NumericalComparisonResult {
618 pub method_results: std::collections::HashMap<NumericalMethod, MethodResult>,
620 pub cross_method_comparison: CrossMethodComparison,
622 pub statistics: Option<StatisticalAnalysis>,
624 pub benchmarks: Option<BenchmarkResults>,
626 pub assessment: ComparisonAssessment,
628}
629
630#[derive(Debug, Clone)]
632pub struct MethodResult {
633 pub method: NumericalMethod,
635 pub gradients: Vec<f64>,
637 pub eps_used: f64,
639 pub computation_time: std::time::Duration,
641 pub memory_usage: Option<usize>,
643 pub function_evaluations: usize,
645 pub estimated_accuracy: f64,
647}
648
649#[derive(Debug, Clone)]
651pub struct CrossMethodComparison {
652 pub pairwise_differences:
654 std::collections::HashMap<(NumericalMethod, NumericalMethod), Vec<f64>>,
655 pub accuracy_ranking: Vec<(NumericalMethod, f64)>,
657 pub performance_ranking: Vec<(NumericalMethod, f64)>,
659 pub consensus_gradient: Vec<f64>,
661 pub confidence_intervals: Vec<(f64, f64)>,
663}
664
665#[derive(Debug, Clone)]
667pub struct StatisticalAnalysis {
668 pub mean_abs_errors: std::collections::HashMap<NumericalMethod, f64>,
670 pub error_std_devs: std::collections::HashMap<NumericalMethod, f64>,
672 pub correlation_matrix: Vec<Vec<f64>>,
674 pub outliers: Vec<usize>,
676 pub bootstrap_intervals: std::collections::HashMap<NumericalMethod, Vec<(f64, f64)>>,
678 pub significance_tests: Vec<SignificanceTest>,
680}
681
682#[derive(Debug, Clone)]
684pub struct BenchmarkResults {
685 pub avg_times: std::collections::HashMap<NumericalMethod, std::time::Duration>,
687 pub memory_usage: std::collections::HashMap<NumericalMethod, usize>,
689 pub function_evals: std::collections::HashMap<NumericalMethod, usize>,
691 pub efficiency_scores: std::collections::HashMap<NumericalMethod, f64>,
693 pub throughput: std::collections::HashMap<NumericalMethod, f64>,
695}
696
697#[derive(Debug, Clone)]
699pub struct ComparisonAssessment {
700 pub recommended_method: NumericalMethod,
702 pub confidence_score: f64,
704 pub method_reliability: std::collections::HashMap<NumericalMethod, f64>,
706 pub quality_indicators: QualityIndicators,
708 pub warnings: Vec<String>,
710 pub summary: String,
712}
713
714#[derive(Debug, Clone)]
716pub struct QualityIndicators {
717 pub smoothness: f64,
719 pub stability: f64,
721 pub consistency: f64,
723 pub conditioning: f64,
725}
726
727#[derive(Debug, Clone)]
729pub struct SignificanceTest {
730 pub test_name: String,
732 pub methods: (NumericalMethod, NumericalMethod),
734 pub statistic: f64,
736 pub p_value: f64,
738 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 pub fn new() -> Self {
784 Self::with_config(NumericalComparisonConfig::default())
785 }
786
787 pub fn with_config(config: NumericalComparisonConfig) -> Self {
789 Self { config }
790 }
791
792 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 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 let cross_method_comparison = self.perform_cross_method_comparison(&method_results)?;
814
815 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 let benchmarks = if self.config.benchmarking.enabled {
824 Some(self.perform_benchmarking(&method_results)?)
825 } else {
826 None
827 };
828
829 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 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 let eps = self.find_optimal_step_size(&func, inputs, method)?;
863
864 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(); 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, function_evaluations: function_evals,
902 estimated_accuracy,
903 })
904 }
905
906 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 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 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 let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
944 Ok(vec![1.0; total_elements])
945 }
946
947 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 let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
961 Ok(vec![0.9; total_elements])
962 }
963
964 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 let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
978 Ok(vec![1.1; total_elements])
979 }
980
981 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 let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
995 Ok(vec![1.05; total_elements])
996 }
997
998 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 let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
1012 Ok(vec![1.02; total_elements])
1013 }
1014
1015 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 let total_elements: usize = inputs.iter().map(|t| t.to_vec().len()).sum();
1030 Ok(vec![1.01; total_elements])
1031 }
1032
1033 fn estimate_accuracy(&self, method: NumericalMethod, eps: f64) -> f64 {
1035 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 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 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 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 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 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 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 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 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 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 let num_methods = method_results.len();
1145 let correlation_matrix = vec![vec![1.0; num_methods]; num_methods];
1146
1147 let outliers = Vec::new();
1149
1150 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 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 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 let recommended_method = cross_method_comparison
1210 .accuracy_ranking
1211 .first()
1212 .map(|(method, _)| *method)
1213 .unwrap_or(NumericalMethod::Central);
1214
1215 let confidence_score = 0.85; let mut method_reliability = std::collections::HashMap::new();
1220 for method in method_results.keys() {
1221 method_reliability.insert(*method, 0.8); }
1223
1224 let quality_indicators = QualityIndicators {
1226 smoothness: 0.9,
1227 stability: 0.85,
1228 consistency: 0.8,
1229 conditioning: 0.75,
1230 };
1231
1232 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 let elements = checker.select_elements_to_check(10);
1298 assert_eq!(elements.len(), 10);
1299 assert_eq!(elements, (0..10).collect::<Vec<_>>());
1300
1301 let elements = checker.select_elements_to_check(1000);
1303 assert!(elements.len() <= 100); 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]); }
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]); }
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); 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); 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 }
1405}