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
52#[allow(dead_code)]
53pub fn compute_step_sizes(x: &ArrayView1<f64>, options: &SparseFiniteDiffOptions) -> Vec<f64> {
54    let n = x.len();
55    let mut h = vec![0.0; n];
56
57    // Typical values for different methods
58    let typical_rel_step = match options.method.as_str() {
59        "2-point" => 1e-8,
60        "3-point" => 1e-5,
61        "cs" => 1e-14, // Complex step
62        _ => 1e-8,     // Default to 2-point
63    };
64
65    let rel_step = options.rel_step.unwrap_or(typical_rel_step);
66    let abs_step = options.abs_step.unwrap_or(1e-8);
67
68    // Calculate step size for each dimension
69    for i in 0..n {
70        let xi = x[i];
71
72        // Base step size on relative and absolute components
73        let mut hi = abs_step.max(rel_step * xi.abs());
74
75        // Ensure step moves in the correct direction for zero elements
76        if xi == 0.0 {
77            hi = abs_step;
78        }
79
80        // Adjust for bounds if provided
81        if let Some(ref bounds) = options.bounds {
82            if i < bounds.len() {
83                let (lower, upper) = bounds[i];
84                if xi + hi > upper {
85                    hi = -hi; // Reverse direction
86                }
87                if xi + hi < lower {
88                    hi = abs_step.min((upper - xi) / 2.0); // Reduce size
89                }
90            }
91        }
92
93        h[i] = hi;
94    }
95
96    h
97}