use nlopt::*;
use std::f64::consts::PI;
fn main() {
bobyqa_example()
}
fn objfn(x: &[f64], _gradient: Option<&mut [f64]>, _params: &mut ()) -> f64 {
let mut f = 0.0;
for i in (4..=x.len()).step_by(2) {
let j1 = i - 2;
for j in (2..=j1).step_by(2) {
let tmpa = x[i - 2] - x[j - 2];
let tmpb = x[i - 1] - x[j - 1];
let tmp = tmpa * tmpa + tmpb * tmpb;
f += 1.0 / tmp.max(1e-6).sqrt()
}
}
f
}
fn bobyqa_example() {
let twopi = 2.0 * PI;
let mut x = [0.0; 100];
let mut xl = [0.0; 100];
let mut xu = [0.0; 100];
let bdl = -1.0;
let bdu = 1.0;
for &m in &[5, 10] {
let q = twopi / m as f64;
let n = 2 * m;
for i in 1..=n {
xl[i - 1] = bdl;
xu[i - 1] = bdu;
}
for &jcase in &[1, 2] {
let npt = if jcase == 2 { 2 * n + 1 } else { n + 6 };
println!("2D output with M={}, N={} and NPT={}", m, n, npt);
for j in 1..=m {
let temp = (j as f64) * q;
x[2 * j - 2] = temp.cos();
x[2 * j - 1] = temp.sin();
}
let mut opt = Nlopt::new(Algorithm::Bobyqa, n, objfn, Target::Minimize, ());
opt.set_lower_bounds(&xl).unwrap();
opt.set_upper_bounds(&xu).unwrap();
opt.set_xtol_rel(1e-6).unwrap();
let res = opt.optimize(&mut x);
println!("Result: {:?}", res);
println!("X vals: {:?}\n", &x[..n]);
}
}
}