cobyla/lib.rs
1#![doc = include_str!("../README.md")]
2
3use nlopt_cobyla::nlopt_constraint;
4
5mod nlopt_cobyla;
6pub use crate::nlopt_cobyla::Func;
7
8use crate::nlopt_cobyla::{
9 NLoptConstraintCfg,
10 NLoptFunctionCfg,
11 cobyla_minimize,
12 nlopt_constraint_raw_callback, // nlopt_eval_constraint,
13 nlopt_function_raw_callback,
14 nlopt_stopping,
15};
16
17use std::os::raw::c_void;
18
19/// Failed termination status of the optimization process
20#[derive(Debug, Clone, Copy)]
21pub enum FailStatus {
22 Failure,
23 InvalidArgs,
24 OutOfMemory,
25 RoundoffLimited,
26 ForcedStop,
27 UnexpectedError,
28}
29
30/// Successful termination status of the optimization process
31#[derive(Debug, Clone, Copy)]
32pub enum SuccessStatus {
33 Success,
34 StopValReached,
35 FtolReached,
36 XtolReached,
37 MaxEvalReached,
38 MaxTimeReached,
39}
40
41/// Outcome when optimization process fails
42type FailOutcome = (FailStatus, Vec<f64>, f64);
43/// Outcome when optimization process succeeds
44type SuccessOutcome = (SuccessStatus, Vec<f64>, f64);
45
46/// Tolerances used as termination criteria.
47/// For all, condition is disabled if value is not strictly positive.
48/// ```rust
49/// # use crate::cobyla::StopTols;
50/// let stop_tol = StopTols {
51/// ftol_rel: 1e-4,
52/// xtol_abs: vec![1e-3; 3], // size should be equal to x dim
53/// ..StopTols::default() // default stop conditions are disabled
54/// };
55/// ```
56#[derive(Debug, Clone, Default)]
57pub struct StopTols {
58 /// Relative tolerance on function value, algorithm stops when `func(x)` changes by less than `ftol_rel * func(x)`
59 pub ftol_rel: f64,
60 /// Absolute tolerance on function value, algorithm stops when `func(x)` change is less than `ftol_rel`
61 pub ftol_abs: f64,
62 /// Relative tolerance on optimization parameters, algorithm stops when all `x[i]` changes by less than `xtol_rel * x[i]`
63 pub xtol_rel: f64,
64 /// Relative tolerance on optimization parameters, algorithm stops when `x[i]` changes by less than `xtol_abs[i]`
65 pub xtol_abs: Vec<f64>,
66}
67
68/// An enum for specifying the initial change of x which correspond to the `rhobeg`
69/// argument of the original Powell's algorithm (hence the name)
70pub enum RhoBeg {
71 /// Used when all x components changes are specified with a single given value
72 All(f64),
73 /// Used to set the components with the given x-dim-sized vector
74 Set(Vec<f64>),
75}
76
77/// Minimizes a function using the Constrained Optimization By Linear Approximation (COBYLA) method.
78///
79/// ## Arguments
80///
81/// * `func` - the function to minimize
82/// * `xinit` - n-vector the initial guess
83/// * `bounds` - x domain specified as a n-vector of tuple `(lower bound, upper bound)`
84/// * `cons` - slice of constraint function intended to be negative at the end
85/// * `args` - user data pass to objective and constraint functions
86/// * `maxeval` - maximum number of objective function evaluation
87/// * `rhobeg`- initial changes to the x component
88///
89/// ## Returns
90///
91/// The status of the optimization process, the argmin value and the objective function value
92///
93/// ## Panics
94///
95/// When some vector arguments like `bounds`, `xtol_abs` do not have the same size as `xinit`
96///
97/// ## Implementation note:
98///
99/// This implementation is a translation of [NLopt](https://github.com/stevengj/nlopt) 2.7.1
100/// See also [NLopt SLSQP](https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#slsqp) documentation.
101///
102/// ## Example
103/// ```
104/// # use approx::assert_abs_diff_eq;
105/// use cobyla::{minimize, Func, RhoBeg};
106///
107/// fn paraboloid(x: &[f64], _data: &mut ()) -> f64 {
108/// 10. * (x[0] + 1.).powf(2.) + x[1].powf(2.)
109/// }
110///
111/// let mut x = vec![1., 1.];
112///
113/// // Constraints definition to be positive eventually: here `x_0 > 0`
114/// let cstr1 = |x: &[f64], _user_data: &mut ()| x[0];
115/// let cons: Vec<&dyn Func<()>> = vec![&cstr1];
116///
117/// match minimize(
118/// paraboloid,
119/// &mut x,
120/// &[(-10., 10.), (-10., 10.)],
121/// &cons,
122/// (),
123/// 200,
124/// RhoBeg::All(0.5),
125/// None
126/// ) {
127/// Ok((status, x_opt, y_opt)) => {
128/// println!("status = {:?}", status);
129/// println!("x_opt = {:?}", x_opt);
130/// println!("y_opt = {}", y_opt);
131/// # assert_abs_diff_eq!(y_opt, 10.0);
132/// }
133/// Err((e, _, _)) => println!("Optim error: {:?}", e),
134/// }
135/// ```
136///
137/// ## Algorithm description:
138///
139/// COBYLA minimizes an objective function F(X) subject to M inequality
140/// constraints on X, where X is a vector of variables that has N components.
141///
142/// The algorithm employs linear approximations to the objective and constraint
143/// functions, the approximations being formed by linear interpolation at N+1
144/// points in the space of the variables. We regard these interpolation points
145/// as vertices of a simplex.
146///
147/// The parameter RHO controls the size of the simplex and it is reduced
148/// automatically from RHOBEG to RHOEND. For each RHO the subroutine tries
149/// to achieve a good vector of variables for the current size, and then RHO
150/// is reduced until the value RHOEND is reached.
151///
152/// Therefore RHOBEG and RHOEND should be set to reasonable initial changes to and the
153/// required accuracy in the variables respectively, but this accuracy should be
154/// viewed as a subject for experimentation because it is not guaranteed.
155///
156/// The subroutine has an advantage over many of its competitors, however, which is
157/// that it treats each constraint individually when calculating a change to the
158/// variables, instead of lumping the constraints together into a single penalty
159/// function.
160///
161/// The name of the algorithm is derived from the phrase Constrained
162/// Optimization BY Linear Approximations.
163///
164/// The user can set the values of RHOBEG and RHOEND, and must provide an
165/// initial vector of variables in X. Further, the value of IPRINT should be
166/// set to 0, 1, 2 or 3, which controls the amount of printing during the
167/// calculation. Specifically, there is no output if IPRINT=0 and there is
168/// output only at the end of the calculation if IPRINT=1.
169/// Otherwise each new value of RHO and SIGMA is printed.
170///
171/// Further, the vector of variables and some function information are
172/// given either when RHO is reduced or when each
173/// new value of F(X) is computed in the cases IPRINT=2 or IPRINT=3
174/// respectively. Here SIGMA is a penalty parameter, it being assumed that a
175/// change to X is an improvement if it reduces the merit function:
176///
177/// F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)),
178///
179/// where C1, C2, ..., CM denote the constraint functions that should become
180/// nonnegative eventually, at least to the precision of RHOEND. In the
181/// printed output the displayed term that is multiplied by SIGMA is
182/// called MAXCV, which stands for 'MAXimum Constraint Violation'.
183///
184/// This implementation is a translation/adaptation of [NLopt](https://github.com/stevengj/nlopt) 2.7.1
185/// See [NLopt COBYLA](https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#cobyla-constrained-optimization-by-linear-approximations) documentation.
186#[allow(clippy::useless_conversion)]
187#[allow(clippy::too_many_arguments)]
188pub fn minimize<F: Func<U>, G: Func<U>, U: Clone>(
189 func: F,
190 xinit: &[f64],
191 bounds: &[(f64, f64)],
192 cons: &[G],
193 args: U,
194 maxeval: usize,
195 rhobeg: RhoBeg,
196 stop_tol: Option<StopTols>,
197) -> Result<SuccessOutcome, FailOutcome> {
198 let fn_cfg = Box::new(NLoptFunctionCfg {
199 objective_fn: func,
200 user_data: args.clone(),
201 });
202 let fn_cfg_ptr = Box::into_raw(fn_cfg) as *mut c_void;
203 let mut cstr_tol = 0.0; // no cstr tolerance
204
205 let mut cstr_cfg = cons
206 .iter()
207 .map(|c| {
208 let c_cfg = Box::new(NLoptConstraintCfg {
209 constraint_fn: c as &dyn Func<U>,
210 user_data: args.clone(),
211 });
212 let c_cfg_ptr = Box::into_raw(c_cfg) as *mut c_void;
213
214 nlopt_constraint {
215 m: 1,
216 f: Some(nlopt_constraint_raw_callback::<F, U>),
217 pre: None,
218 mf: None,
219 f_data: c_cfg_ptr,
220 tol: &mut cstr_tol,
221 }
222 })
223 .collect::<Vec<_>>();
224
225 let mut x = vec![0.; xinit.len()];
226 x.copy_from_slice(xinit);
227 let n = x.len() as u32;
228 let m = cons.len() as u32;
229
230 if bounds.len() != x.len() {
231 panic!(
232 "{}",
233 format!(
234 "Minimize Error: Bounds and x should have same size! Got {} for bounds and {} for x.",
235 bounds.len(),
236 x.len()
237 )
238 )
239 }
240 let lbs: Vec<f64> = bounds.iter().map(|b| b.0).collect();
241 let ubs: Vec<f64> = bounds.iter().map(|b| b.1).collect();
242 let x_weights = vec![0.; n as usize];
243 let mut dx = match rhobeg {
244 RhoBeg::All(val) => vec![val; n as usize],
245 RhoBeg::Set(val) => {
246 if val.len() != n as usize {
247 panic!(
248 "{}",
249 format!(
250 "Minimize Error: xtol_abs should have x dim size ({}), got {}",
251 n,
252 val.len()
253 )
254 )
255 } else {
256 val
257 }
258 }
259 };
260 let mut minf = f64::INFINITY;
261 let mut nevals_p = 0;
262 let mut force_stop = 0;
263
264 let stop_tol = stop_tol.unwrap_or_default();
265 let xtol_abs = if stop_tol.xtol_abs.is_empty() {
266 std::ptr::null()
267 } else if stop_tol.xtol_abs.len() != n as usize {
268 panic!(
269 "{}",
270 format!(
271 "Minimize Error: xtol_abs should have x dim size ({}), got {}",
272 n,
273 stop_tol.xtol_abs.len()
274 )
275 );
276 } else {
277 stop_tol.xtol_abs.as_ptr()
278 };
279 let mut stop = nlopt_stopping {
280 n,
281 minf_max: -f64::INFINITY,
282 ftol_rel: stop_tol.ftol_rel,
283 ftol_abs: stop_tol.ftol_abs,
284 xtol_rel: stop_tol.xtol_rel,
285 xtol_abs,
286 x_weights: x_weights.as_ptr(),
287 nevals_p: &mut nevals_p,
288 maxeval: maxeval as i32,
289 maxtime: 0.0,
290 start: 0.0,
291 force_stop: &mut force_stop,
292 stop_msg: "".to_string(),
293 };
294
295 let status = unsafe {
296 cobyla_minimize::<U>(
297 n.into(),
298 Some(nlopt_function_raw_callback::<F, U>),
299 fn_cfg_ptr,
300 m.into(),
301 cstr_cfg.as_mut_ptr(),
302 0,
303 std::ptr::null_mut(),
304 lbs.as_ptr(),
305 ubs.as_ptr(),
306 x.as_mut_ptr(),
307 &mut minf,
308 &mut stop,
309 dx.as_mut_ptr(),
310 )
311 };
312
313 // Convert the raw pointer back into a Box with the B::from_raw function,
314 // allowing the Box destructor to perform the cleanup.
315 unsafe {
316 let _ = Box::from_raw(fn_cfg_ptr as *mut NLoptFunctionCfg<F, U>);
317 };
318
319 match status {
320 -1 => Err((FailStatus::Failure, x, minf)),
321 -2 => Err((FailStatus::InvalidArgs, x, minf)),
322 -3 => Err((FailStatus::OutOfMemory, x, minf)),
323 -4 => Err((FailStatus::RoundoffLimited, x, minf)),
324 -5 => Err((FailStatus::ForcedStop, x, minf)),
325 1 => Ok((SuccessStatus::Success, x, minf)),
326 2 => Ok((SuccessStatus::StopValReached, x, minf)),
327 3 => Ok((SuccessStatus::FtolReached, x, minf)),
328 4 => Ok((SuccessStatus::XtolReached, x, minf)),
329 5 => Ok((SuccessStatus::MaxEvalReached, x, minf)),
330 6 => Ok((SuccessStatus::MaxTimeReached, x, minf)),
331 _ => Err((FailStatus::UnexpectedError, x, minf)),
332 }
333}
334
335#[cfg(test)]
336mod tests {
337 use super::*;
338 use approx::assert_abs_diff_eq;
339
340 use crate::nlopt_cobyla::cobyla_minimize;
341 use crate::nlopt_cobyla::nlopt_stopping;
342
343 /////////////////////////////////////////////////////////////////////////
344 // Second problem (see cobyla.c case 6)
345
346 fn raw_paraboloid(
347 _n: libc::c_uint,
348 x: *const libc::c_double,
349 _gradient: *mut libc::c_double,
350 _func_data: *mut libc::c_void,
351 ) -> libc::c_double {
352 unsafe {
353 let r1 = *x.offset(0) + 1.0;
354 let r2 = *x.offset(1);
355 10.0 * (r1 * r1) + (r2 * r2) as libc::c_double
356 }
357 }
358
359 #[test]
360 fn test_cobyla_minimize() {
361 let mut x = vec![1., 1.];
362 let mut lb = vec![-10.0, -10.0];
363 let mut ub = vec![10.0, 10.0];
364 let mut dx = vec![0.5, 0.5];
365 let mut minf = f64::INFINITY;
366 let mut nevals_p = 0;
367 let mut force_stop = 0;
368
369 let mut stop = nlopt_stopping {
370 n: 2,
371 minf_max: -f64::INFINITY,
372 ftol_rel: 0.0,
373 ftol_abs: 0.0,
374 xtol_rel: 0.0,
375 xtol_abs: std::ptr::null(),
376 x_weights: std::ptr::null(),
377 nevals_p: &mut nevals_p,
378 maxeval: 1000,
379 maxtime: 0.0,
380 start: 0.0,
381 force_stop: &mut force_stop,
382 stop_msg: "".to_string(),
383 };
384
385 let res = unsafe {
386 cobyla_minimize::<()>(
387 2,
388 Some(raw_paraboloid),
389 std::ptr::null_mut(),
390 0,
391 std::ptr::null_mut(),
392 0,
393 std::ptr::null_mut(),
394 lb.as_mut_ptr(),
395 ub.as_mut_ptr(),
396 x.as_mut_ptr(),
397 &mut minf,
398 &mut stop,
399 dx.as_mut_ptr(),
400 )
401 };
402
403 println!("status = {:?}", res);
404 println!("x = {:?}", x);
405
406 assert_abs_diff_eq!(x.as_slice(), [-1., 0.].as_slice(), epsilon = 1e-4);
407 }
408
409 fn paraboloid(x: &[f64], _data: &mut ()) -> f64 {
410 10. * (x[0] + 1.).powf(2.) + x[1].powf(2.)
411 }
412
413 #[test]
414 fn test_paraboloid() {
415 let xinit = vec![1., 1.];
416
417 // let mut cons: Vec<&dyn Func<()>> = vec![];
418 let mut cons: Vec<&dyn Func<()>> = vec![];
419 let cstr1 = |x: &[f64], _user_data: &mut ()| x[0];
420 cons.push(&cstr1 as &dyn Func<()>);
421
422 // x_opt = [0, 0]
423 match minimize(
424 paraboloid,
425 &xinit,
426 &[(-10., 10.), (-10., 10.)],
427 &cons,
428 (),
429 200,
430 RhoBeg::All(0.5),
431 None,
432 ) {
433 Ok((_, x, _)) => {
434 let exp = [0., 0.];
435 for (act, exp) in x.iter().zip(exp.iter()) {
436 assert_abs_diff_eq!(act, exp, epsilon = 1e-3);
437 }
438 }
439 Err((status, _, _)) => {
440 panic!("{}", format!("Error status : {:?}", status));
441 }
442 }
443 }
444
445 fn fletcher9115(x: &[f64], _user_data: &mut ()) -> f64 {
446 -x[0] - x[1]
447 }
448
449 fn cstr1(x: &[f64], _user_data: &mut ()) -> f64 {
450 x[1] - x[0] * x[0]
451 }
452 fn cstr2(x: &[f64], _user_data: &mut ()) -> f64 {
453 1. - x[0] * x[0] - x[1] * x[1]
454 }
455
456 #[test]
457 fn test_fletcher9115() {
458 let xinit = vec![1., 1.];
459
460 let cons = vec![&cstr1 as &dyn Func<()>, &cstr2 as &dyn Func<()>];
461
462 let stop_tol = StopTols {
463 ftol_rel: 1e-4,
464 xtol_rel: 1e-4,
465 ..StopTols::default()
466 };
467
468 match minimize(
469 fletcher9115,
470 &xinit,
471 &[(-10., 10.), (-10., 10.)],
472 &cons,
473 (),
474 200,
475 RhoBeg::All(0.5),
476 Some(stop_tol),
477 ) {
478 Ok((_, x, _)) => {
479 let sqrt_0_5: f64 = 0.5_f64.sqrt();
480 let exp = [sqrt_0_5, sqrt_0_5];
481 for (act, exp) in x.iter().zip(exp.iter()) {
482 assert_abs_diff_eq!(act, exp, epsilon = 1e-3);
483 }
484 }
485 Err((status, _, _)) => {
486 println!("Error status : {:?}", status);
487 panic!("Test fail");
488 }
489 }
490 }
491
492 fn xsinx(x: &[f64], _user_data: &mut ()) -> f64 {
493 //(x - 3.5) * ((x - 3.5) / std::f64::consts::PI).mapv(|v| v.sin())
494 (x[0] - 3.5) * f64::sin((x[0] - 3.5) / std::f64::consts::PI)
495 }
496
497 #[test]
498 fn test_xsinx() {
499 let xinit = vec![10.];
500
501 // let mut cons: Vec<&dyn Func<()>> = vec![];
502 let mut cons: Vec<&dyn Func<()>> = vec![];
503 let cstr1 = |x: &[f64], _user_data: &mut ()| 17. - x[0];
504 cons.push(&cstr1 as &dyn Func<()>);
505
506 // x_opt = [0, 0]
507 match minimize(
508 xsinx,
509 &xinit,
510 &[(0., 25.)],
511 &cons,
512 (),
513 200,
514 RhoBeg::All(0.5),
515 None,
516 ) {
517 Ok((_, x, _)) => {
518 //let exp = [18.935];
519 let exp = [17.];
520 for (act, exp) in x.iter().zip(exp.iter()) {
521 assert_abs_diff_eq!(act, exp, epsilon = 1e-2);
522 }
523 }
524 Err((status, _, _)) => {
525 panic!("{}", format!("Error status : {:?}", status));
526 }
527 }
528 }
529}