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}