scirs2_optimize/unconstrained/
utils.rs

1//! Common utilities for unconstrained optimization algorithms
2
3use crate::error::OptimizeError;
4use ndarray::{Array1, Array2, ArrayView1};
5
6/// Computes finite difference gradient
7#[allow(dead_code)]
8pub fn finite_difference_gradient<F, S>(
9    fun: &mut F,
10    x: &ArrayView1<f64>,
11    step: f64,
12) -> Result<Array1<f64>, OptimizeError>
13where
14    F: FnMut(&ArrayView1<f64>) -> S,
15    S: Into<f64>,
16{
17    let n = x.len();
18    let mut grad = Array1::<f64>::zeros(n);
19    let mut x_plus = x.to_owned();
20    let mut x_minus = x.to_owned();
21
22    for i in 0..n {
23        let h = step * (1.0 + x[i].abs());
24        x_plus[i] = x[i] + h;
25        x_minus[i] = x[i] - h;
26
27        let f_plus = fun(&x_plus.view()).into();
28        let f_minus = fun(&x_minus.view()).into();
29
30        if !f_plus.is_finite() || !f_minus.is_finite() {
31            return Err(OptimizeError::ComputationError(
32                "Function returned non-finite value during gradient computation".to_string(),
33            ));
34        }
35
36        grad[i] = (f_plus - f_minus) / (2.0 * h);
37
38        // Reset for next iteration
39        x_plus[i] = x[i];
40        x_minus[i] = x[i];
41    }
42
43    Ok(grad)
44}
45
46/// Computes finite difference Hessian
47#[allow(dead_code)]
48pub fn finite_difference_hessian<F, S>(
49    fun: &mut F,
50    x: &ArrayView1<f64>,
51    step: f64,
52) -> Result<Array2<f64>, OptimizeError>
53where
54    F: FnMut(&ArrayView1<f64>) -> S,
55    S: Into<f64>,
56{
57    let n = x.len();
58    let mut hess = Array2::<f64>::zeros((n, n));
59    let mut x_temp = x.to_owned();
60
61    let f0 = fun(&x.view()).into();
62
63    for i in 0..n {
64        let hi = step * (1.0 + x[i].abs());
65
66        // Diagonal elements
67        x_temp[i] = x[i] + hi;
68        let fp = fun(&x_temp.view()).into();
69        x_temp[i] = x[i] - hi;
70        let fm = fun(&x_temp.view()).into();
71        x_temp[i] = x[i];
72
73        hess[[i, i]] = (fp - 2.0 * f0 + fm) / (hi * hi);
74
75        // Off-diagonal elements
76        for j in (i + 1)..n {
77            let hj = step * (1.0 + x[j].abs());
78
79            x_temp[i] = x[i] + hi;
80            x_temp[j] = x[j] + hj;
81            let fpp = fun(&x_temp.view()).into();
82
83            x_temp[i] = x[i] + hi;
84            x_temp[j] = x[j] - hj;
85            let fpm = fun(&x_temp.view()).into();
86
87            x_temp[i] = x[i] - hi;
88            x_temp[j] = x[j] + hj;
89            let fmp = fun(&x_temp.view()).into();
90
91            x_temp[i] = x[i] - hi;
92            x_temp[j] = x[j] - hj;
93            let fmm = fun(&x_temp.view()).into();
94
95            x_temp[i] = x[i];
96            x_temp[j] = x[j];
97
98            let hess_ij = (fpp - fpm - fmp + fmm) / (4.0 * hi * hj);
99            hess[[i, j]] = hess_ij;
100            hess[[j, i]] = hess_ij;
101        }
102    }
103
104    Ok(hess)
105}
106
107/// Check convergence criteria
108#[allow(dead_code)]
109pub fn check_convergence(
110    f_delta: f64,
111    x_delta: f64,
112    g_norm: f64,
113    ftol: f64,
114    xtol: f64,
115    gtol: f64,
116) -> bool {
117    f_delta.abs() < ftol || x_delta < xtol || g_norm < gtol
118}
119
120/// Compute the norm of the difference between two arrays
121#[allow(dead_code)]
122pub fn array_diff_norm(x1: &ArrayView1<f64>, x2: &ArrayView1<f64>) -> f64 {
123    (x1 - x2).mapv(|x| x.powi(2)).sum().sqrt()
124}
125
126/// Clip step size to satisfy bounds
127#[allow(dead_code)]
128pub fn clip_step(
129    x: &ArrayView1<f64>,
130    direction: &ArrayView1<f64>,
131    alpha: f64,
132    lower: &[Option<f64>],
133    upper: &[Option<f64>],
134) -> f64 {
135    let mut clipped_alpha = alpha;
136
137    for i in 0..x.len() {
138        if direction[i] != 0.0 {
139            // Check lower bound
140            if let Some(lb) = lower[i] {
141                if direction[i] < 0.0 {
142                    let max_step = (lb - x[i]) / direction[i];
143                    if max_step >= 0.0 {
144                        clipped_alpha = clipped_alpha.min(max_step);
145                    }
146                }
147            }
148
149            // Check upper bound
150            if let Some(ub) = upper[i] {
151                if direction[i] > 0.0 {
152                    let max_step = (ub - x[i]) / direction[i];
153                    if max_step >= 0.0 {
154                        clipped_alpha = clipped_alpha.min(max_step);
155                    }
156                }
157            }
158        }
159    }
160
161    clipped_alpha.max(0.0)
162}
163
164/// Convert between Array1 and ArrayView1 consistently
165#[allow(dead_code)]
166pub fn to_array_view<T>(arr: &Array1<T>) -> ArrayView1<T> {
167    arr.view()
168}
169
170/// Initialize step size for line search
171#[allow(dead_code)]
172pub fn initial_step_size(_grad_norm: f64, max_step: Option<f64>) -> f64 {
173    let default_step = if _grad_norm > 0.0 {
174        1.0 / _grad_norm
175    } else {
176        1.0
177    };
178
179    if let Some(max_s) = max_step {
180        default_step.min(max_s)
181    } else {
182        default_step
183    }
184}
185
186#[cfg(test)]
187mod tests {
188    use super::*;
189    use approx::assert_abs_diff_eq;
190
191    #[test]
192    fn test_finite_difference_gradient() {
193        let mut quadratic = |x: &ArrayView1<f64>| -> f64 { x[0] * x[0] + 2.0 * x[1] * x[1] };
194
195        let x = Array1::from_vec(vec![1.0, 2.0]);
196        let grad = finite_difference_gradient(&mut quadratic, &x.view(), 1e-8).unwrap();
197
198        assert_abs_diff_eq!(grad[0], 2.0, epsilon = 1e-6);
199        assert_abs_diff_eq!(grad[1], 8.0, epsilon = 1e-6);
200    }
201}
scirs2_optimize/unconstrained/utils.rs

scirs2_optimize/unconstrained/
utils.rs
