scirs2_optimize/sparse_numdiff/finite_diff.rs
1//! Finite difference methods for sparse numerical differentiation
2//!
3//! This module provides utilities and options for finite difference
4//! methods used in sparse numerical differentiation.
5
6use crate::parallel::ParallelOptions;
7use ndarray::ArrayView1;
8
9/// Options for sparse numerical differentiation
10#[derive(Debug, Clone)]
11pub struct SparseFiniteDiffOptions {
12 /// Method to use for finite differences ('2-point', '3-point', 'cs')
13 pub method: String,
14 /// Relative step size (if None, determined automatically)
15 pub rel_step: Option<f64>,
16 /// Absolute step size (if None, determined automatically)
17 pub abs_step: Option<f64>,
18 /// Bounds on the variables
19 pub bounds: Option<Vec<(f64, f64)>>,
20 /// Parallel computation options
21 pub parallel: Option<ParallelOptions>,
22 /// Random seed for coloring algorithm
23 pub seed: Option<u64>,
24 /// Maximum number of columns to group together
25 pub max_group_size: usize,
26}
27
28impl Default for SparseFiniteDiffOptions {
29 fn default() -> Self {
30 Self {
31 method: "2-point".to_string(),
32 rel_step: None,
33 abs_step: None,
34 bounds: None,
35 parallel: None,
36 seed: None,
37 max_group_size: 100,
38 }
39 }
40}
41
42/// Calculates appropriate step sizes for finite differences
43///
44/// # Arguments
45///
46/// * `x` - Point at which to compute derivatives
47/// * `options` - Finite difference options
48///
49/// # Returns
50///
51/// Vector of step sizes for each dimension
52pub fn compute_step_sizes(x: &ArrayView1<f64>, options: &SparseFiniteDiffOptions) -> Vec<f64> {
53 let n = x.len();
54 let mut h = vec![0.0; n];
55
56 // Typical values for different methods
57 let typical_rel_step = match options.method.as_str() {
58 "2-point" => 1e-8,
59 "3-point" => 1e-5,
60 "cs" => 1e-14, // Complex step
61 _ => 1e-8, // Default to 2-point
62 };
63
64 let rel_step = options.rel_step.unwrap_or(typical_rel_step);
65 let abs_step = options.abs_step.unwrap_or(1e-8);
66
67 // Calculate step size for each dimension
68 for i in 0..n {
69 let xi = x[i];
70
71 // Base step size on relative and absolute components
72 let mut hi = abs_step.max(rel_step * xi.abs());
73
74 // Ensure step moves in the correct direction for zero elements
75 if xi == 0.0 {
76 hi = abs_step;
77 }
78
79 // Adjust for bounds if provided
80 if let Some(ref bounds) = options.bounds {
81 if i < bounds.len() {
82 let (lower, upper) = bounds[i];
83 if xi + hi > upper {
84 hi = -hi; // Reverse direction
85 }
86 if xi + hi < lower {
87 hi = abs_step.min((upper - xi) / 2.0); // Reduce size
88 }
89 }
90 }
91
92 h[i] = hi;
93 }
94
95 h
96}