scirs2_optimize/unconstrained/
mod.rs

1//! Unconstrained optimization algorithms
2//!
3//! This module provides various algorithms for unconstrained minimization problems.
4
5use crate::error::OptimizeError;
6use ndarray::{Array1, ArrayView1};
7use std::fmt;
8
9// Sub-modules
10pub mod bfgs;
11pub mod conjugate_gradient;
12pub mod lbfgs;
13pub mod line_search;
14pub mod nelder_mead;
15pub mod newton;
16pub mod powell;
17pub mod result;
18pub mod trust_region;
19pub mod utils;
20
21// Import result type
22pub use result::OptimizeResult;
23
24// Re-export commonly used items
25pub use bfgs::minimize_bfgs;
26pub use conjugate_gradient::minimize_conjugate_gradient;
27pub use lbfgs::{minimize_lbfgs, minimize_lbfgsb};
28pub use nelder_mead::minimize_nelder_mead;
29pub use newton::minimize_newton_cg;
30pub use powell::minimize_powell;
31pub use trust_region::{minimize_trust_exact, minimize_trust_krylov, minimize_trust_ncg};
32
33/// Optimization methods for unconstrained minimization.
34#[derive(Debug, Clone, Copy)]
35pub enum Method {
36    /// Nelder-Mead simplex method
37    NelderMead,
38    /// Powell's method
39    Powell,
40    /// Conjugate gradient method
41    CG,
42    /// BFGS quasi-Newton method
43    BFGS,
44    /// Limited-memory BFGS method
45    LBFGS,
46    /// Limited-memory BFGS method with bounds support
47    LBFGSB,
48    /// Newton's method with conjugate gradient solver
49    NewtonCG,
50    /// Trust-region Newton method with conjugate gradient solver
51    TrustNCG,
52    /// Trust-region method with Krylov subproblem solver
53    TrustKrylov,
54    /// Trust-region method with exact subproblem solver
55    TrustExact,
56}
57
58/// Bounds for optimization variables
59#[derive(Debug, Clone)]
60pub struct Bounds {
61    /// Lower bounds
62    pub lower: Vec<Option<f64>>,
63    /// Upper bounds
64    pub upper: Vec<Option<f64>>,
65}
66
67/// Options for optimization algorithms
68#[derive(Debug, Clone)]
69pub struct Options {
70    /// Maximum number of iterations
71    pub max_iter: usize,
72    /// Maximum number of function evaluations
73    pub max_fev: Option<usize>,
74    /// Function tolerance for convergence
75    pub ftol: f64,
76    /// Change tolerance for convergence
77    pub xtol: f64,
78    /// Gradient tolerance for convergence
79    pub gtol: f64,
80    /// Initial step size
81    pub initial_step: Option<f64>,
82    /// Maximum step size for line search
83    pub maxstep: Option<f64>,
84    /// Whether to use finite differences for gradient
85    pub finite_diff: bool,
86    /// Finite difference step size
87    pub eps: f64,
88    /// Initial trust-region radius for trust-region methods
89    pub trust_radius: Option<f64>,
90    /// Maximum trust-region radius for trust-region methods
91    pub max_trust_radius: Option<f64>,
92    /// Minimum trust-region radius for trust-region methods
93    pub min_trust_radius: Option<f64>,
94    /// Tolerance for the trust-region subproblem
95    pub trust_tol: Option<f64>,
96    /// Maximum iterations for trust-region subproblem
97    pub trust_max_iter: Option<usize>,
98    /// Threshold for accepting a step in the trust-region method
99    pub trust_eta: Option<f64>,
100    /// Bounds constraints for variables
101    pub bounds: Option<Bounds>,
102}
103
104// Implement Display for Method
105impl fmt::Display for Method {
106    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
107        match self {
108            Method::NelderMead => write!(f, "Nelder-Mead"),
109            Method::Powell => write!(f, "Powell"),
110            Method::CG => write!(f, "Conjugate Gradient"),
111            Method::BFGS => write!(f, "BFGS"),
112            Method::LBFGS => write!(f, "L-BFGS"),
113            Method::LBFGSB => write!(f, "L-BFGS-B"),
114            Method::NewtonCG => write!(f, "Newton-CG"),
115            Method::TrustNCG => write!(f, "Trust-NCG"),
116            Method::TrustKrylov => write!(f, "Trust-Krylov"),
117            Method::TrustExact => write!(f, "Trust-Exact"),
118        }
119    }
120}
121
122// Implement Default for Options
123impl Default for Options {
124    fn default() -> Self {
125        Options {
126            max_iter: 1000,
127            max_fev: None,
128            ftol: 1e-8,
129            xtol: 1e-8,
130            gtol: 1e-5,
131            initial_step: None,
132            maxstep: None,
133            finite_diff: false,
134            eps: 1.4901161193847656e-8,
135            trust_radius: Some(1.0),
136            max_trust_radius: Some(100.0),
137            min_trust_radius: Some(1e-10),
138            trust_tol: Some(1e-8),
139            trust_max_iter: Some(100),
140            trust_eta: Some(0.1),
141            bounds: None,
142        }
143    }
144}
145
146// Implement Bounds methods
147impl Bounds {
148    /// Create new bounds from arrays
149    pub fn new(bounds: &[(Option<f64>, Option<f64>)]) -> Self {
150        let (lower, upper): (Vec<_>, Vec<_>) = bounds.iter().cloned().unzip();
151        Self { lower, upper }
152    }
153
154    /// Create bounds from vectors
155    pub fn from_vecs(lb: Vec<Option<f64>>, ub: Vec<Option<f64>>) -> Result<Self, OptimizeError> {
156        if lb.len() != ub.len() {
157            return Err(OptimizeError::ValueError(
158                "Lower and upper bounds must have the same length".to_string(),
159            ));
160        }
161
162        for (l, u) in lb.iter().zip(ub.iter()) {
163            if let (Some(l_val), Some(u_val)) = (l, u) {
164                if l_val > u_val {
165                    return Err(OptimizeError::ValueError(
166                        "Lower bound must be less than or equal to upper bound".to_string(),
167                    ));
168                }
169            }
170        }
171
172        Ok(Self {
173            lower: lb,
174            upper: ub,
175        })
176    }
177
178    /// Check if point is feasible
179    pub fn is_feasible(&self, x: &[f64]) -> bool {
180        if x.len() != self.lower.len() {
181            return false;
182        }
183
184        for (&xi, (&lb, &ub)) in x.iter().zip(self.lower.iter().zip(self.upper.iter())) {
185            if let Some(l) = lb {
186                if xi < l {
187                    return false;
188                }
189            }
190            if let Some(u) = ub {
191                if xi > u {
192                    return false;
193                }
194            }
195        }
196        true
197    }
198
199    /// Project point onto feasible region
200    pub fn project(&self, x: &mut [f64]) {
201        for (xi, (&lb, &ub)) in x.iter_mut().zip(self.lower.iter().zip(self.upper.iter())) {
202            if let Some(l) = lb {
203                if *xi < l {
204                    *xi = l;
205                }
206            }
207            if let Some(u) = ub {
208                if *xi > u {
209                    *xi = u;
210                }
211            }
212        }
213    }
214
215    /// Check if bounds are active
216    pub fn has_bounds(&self) -> bool {
217        self.lower.iter().any(|b| b.is_some()) || self.upper.iter().any(|b| b.is_some())
218    }
219}
220
221/// Main minimize function for unconstrained optimization
222pub fn minimize<F, S>(
223    fun: F,
224    x0: &[f64],
225    method: Method,
226    options: Option<Options>,
227) -> Result<OptimizeResult<S>, OptimizeError>
228where
229    F: FnMut(&ArrayView1<f64>) -> S + Clone,
230    S: Into<f64> + Clone,
231{
232    let options = &options.unwrap_or_default();
233    let x0 = Array1::from_vec(x0.to_vec());
234
235    // Check initial point feasibility if bounds are provided
236    if let Some(ref bounds) = options.bounds {
237        if !bounds.is_feasible(x0.as_slice().unwrap()) {
238            return Err(OptimizeError::ValueError(
239                "Initial point is not feasible".to_string(),
240            ));
241        }
242    }
243
244    match method {
245        Method::NelderMead => nelder_mead::minimize_nelder_mead(fun, x0, options),
246        Method::Powell => powell::minimize_powell(fun, x0, options),
247        Method::CG => conjugate_gradient::minimize_conjugate_gradient(fun, x0, options),
248        Method::BFGS => bfgs::minimize_bfgs(fun, x0, options),
249        Method::LBFGS => lbfgs::minimize_lbfgs(fun, x0, options),
250        Method::LBFGSB => lbfgs::minimize_lbfgsb(fun, x0, options),
251        Method::NewtonCG => newton::minimize_newton_cg(fun, x0, options),
252        Method::TrustNCG => trust_region::minimize_trust_ncg(fun, x0, options),
253        Method::TrustKrylov => trust_region::minimize_trust_krylov(fun, x0, options),
254        Method::TrustExact => trust_region::minimize_trust_exact(fun, x0, options),
255    }
256}
257
258#[cfg(test)]
259mod tests {
260    use super::*;
261    use approx::assert_abs_diff_eq;
262
263    #[test]
264    fn test_simple_quadratic() {
265        let quadratic = |x: &ArrayView1<f64>| -> f64 { x[0] * x[0] + x[1] * x[1] };
266
267        let x0 = vec![1.0, 1.0];
268        let result = minimize(quadratic, &x0, Method::BFGS, None);
269        assert!(result.is_ok());
270
271        let result = result.unwrap();
272        assert!(result.success);
273        assert_abs_diff_eq!(result.x[0], 0.0, epsilon = 1e-6);
274        assert_abs_diff_eq!(result.x[1], 0.0, epsilon = 1e-6);
275    }
276}
scirs2_optimize/unconstrained/mod.rs

scirs2_optimize/unconstrained/
mod.rs
