fcmaes_core/
da.rs

1//! Dual Annealing — Rust port of the C++ `daoptimizer.cpp`.
2//!
3//! Generalized simulated annealing (derived from SciPy's `_dual_annealing`):
4//! a distorted Cauchy-Lorentz visiting distribution, a Markov strategy chain
5//! with generalized accept/reject, re-annealing, and an optional local search.
6//! The C++ used LBFGSpp's L-BFGS-B for the local search; this port uses a
7//! self-contained bounded limited-memory quasi-Newton (projected L-BFGS on the
8//! `[0,1]` box with a finite-difference gradient), avoiding a Fortran
9//! dependency. C++-only in the original (no pure-Python twin); validated by
10//! convergence rather than a reference distribution.
11
12use crate::fitness::Objective;
13use crate::rng::Rng;
14
15/// Outcome of a Dual Annealing run (mirrors the C++ `DaResult`).
16#[derive(Clone, Debug)]
17pub struct DaResult {
18    pub x: Vec<f64>,
19    pub y: f64,
20    pub evaluations: u64,
21    pub iterations: i32,
22    pub stop: i32,
23}
24
25/// Tunable inputs for [`optimize_da`].
26#[derive(Clone, Debug)]
27pub struct DaParams {
28    pub max_evaluations: u64,
29    pub use_local_search: bool,
30    pub seed: u64,
31    pub runid: i64,
32}
33
34impl Default for DaParams {
35    fn default() -> Self {
36        Self {
37            max_evaluations: 100_000,
38            use_local_search: true,
39            seed: 0,
40            runid: 0,
41        }
42    }
43}
44
45const BIG_VALUE: f64 = 1e16;
46const TAIL_LIMIT: f64 = 1e8;
47const MIN_VISIT_BOUND: f64 = 1e-10;
48const MAX_REINIT_COUNT: i32 = 1000;
49
50const TEMPERATURE_START: f64 = 5230.0;
51const QV: f64 = 2.62;
52const QA: f64 = -5.0;
53const MAXSTEPS: i32 = 1000;
54const TEMPERATURE_RESTART: f64 = 0.1;
55
56struct Da<'a, O: Objective> {
57    obj: &'a O,
58    dim: usize,
59    has_bounds: bool,
60    lower: Vec<f64>,
61    scale: Vec<f64>,
62    max_evals: u64,
63    eval_counter: u64,
64    use_local_search: bool,
65    rng: Rng,
66
67    // fitness best (reset per local search)
68    fit_best_y: f64,
69    fit_best_x: Vec<f64>,
70
71    // visiting distribution factors
72    factor4_p: f64,
73    factor6: f64,
74
75    // energy state
76    ebest: f64,
77    xbest: Vec<f64>,
78    current_energy: f64,
79    current_location: Vec<f64>,
80
81    // strategy chain
82    emin: f64,
83    xmin: Vec<f64>,
84    not_improved_idx: i32,
85    not_improved_max_idx: i32,
86    temperature_step: f64,
87    k: f64,
88    state_improved: bool,
89
90    reinit_failed: bool,
91}
92
93impl<'a, O: Objective> Da<'a, O> {
94    fn new(
95        obj: &'a O,
96        dim: usize,
97        lower: Vec<f64>,
98        upper: Vec<f64>,
99        max_evals: u64,
100        use_local_search: bool,
101        seed: u64,
102    ) -> Self {
103        let has_bounds = !lower.is_empty();
104        let scale: Vec<f64> = if has_bounds {
105            upper.iter().zip(&lower).map(|(u, l)| u - l).collect()
106        } else {
107            vec![1.0; dim]
108        };
109        // visiting-distribution invariants
110        let factor2 = ((4.0 - QV) * (QV - 1.0).ln()).exp();
111        let factor3 = ((2.0 - QV) * 2.0_f64.ln() / (QV - 1.0)).exp();
112        let factor4_p = std::f64::consts::PI.sqrt() * factor2 / (factor3 * (3.0 - QV));
113        let factor5 = 1.0 / (QV - 1.0) - 0.5;
114        let d1 = 2.0 - factor5;
115        let factor6 = std::f64::consts::PI * (1.0 - factor5)
116            / (std::f64::consts::PI * (1.0 - factor5)).sin()
117            / libm::lgamma(d1).exp();
118        Da {
119            obj,
120            dim,
121            has_bounds,
122            lower,
123            scale,
124            max_evals,
125            eval_counter: 0,
126            use_local_search,
127            rng: Rng::new(seed),
128            fit_best_y: f64::MAX,
129            fit_best_x: vec![],
130            factor4_p,
131            factor6,
132            ebest: f64::MAX,
133            xbest: vec![],
134            current_energy: f64::MAX,
135            current_location: vec![],
136            emin: f64::MAX,
137            xmin: vec![],
138            not_improved_idx: 0,
139            not_improved_max_idx: 1000,
140            temperature_step: 0.0,
141            k: 100.0 * dim as f64,
142            state_improved: false,
143            reinit_failed: false,
144        }
145    }
146
147    fn closest_feasible(&self, x: &[f64]) -> Vec<f64> {
148        if self.has_bounds {
149            x.iter().map(|&v| v.clamp(-1.0, 1.0)).collect()
150        } else {
151            x.to_vec()
152        }
153    }
154
155    fn encode(&self, x: &[f64]) -> Vec<f64> {
156        if self.has_bounds {
157            (0..self.dim)
158                .map(|i| (x[i] - self.lower[i]) / self.scale[i])
159                .collect()
160        } else {
161            x.to_vec()
162        }
163    }
164
165    fn decode(&self, x: &[f64]) -> Vec<f64> {
166        if self.has_bounds {
167            (0..self.dim)
168                .map(|i| x[i] * self.scale[i] + self.lower[i])
169                .collect()
170        } else {
171            x.to_vec()
172        }
173    }
174
175    fn raw_eval(&mut self, x_decoded: &[f64]) -> f64 {
176        self.eval_counter += 1;
177        self.obj.eval_scalar(x_decoded)
178    }
179
180    /// Evaluate an *encoded* point, tracking the local-search best.
181    fn value(&mut self, x: &[f64]) -> f64 {
182        let res = if self.has_bounds {
183            let feas = self.closest_feasible(x);
184            let dec = self.decode(&feas);
185            self.raw_eval(&dec)
186        } else {
187            self.raw_eval(x)
188        };
189        if res < self.fit_best_y {
190            self.fit_best_y = res;
191            self.fit_best_x = x.to_vec();
192        }
193        res
194    }
195
196    fn max_eval_reached(&self) -> bool {
197        self.eval_counter >= self.max_evals
198    }
199
200    fn normal_vec(&mut self) -> Vec<f64> {
201        (0..self.dim).map(|_| self.rng.gaussian()).collect()
202    }
203    fn uniform_vec(&mut self) -> Vec<f64> {
204        (0..self.dim).map(|_| self.rng.uniform01()).collect()
205    }
206
207    // ---- visiting distribution ----
208
209    fn visit_fn(&mut self, temperature: f64, n: usize) -> Vec<f64> {
210        let x: Vec<f64> = (0..n).map(|_| self.rng.gaussian()).collect();
211        let y: Vec<f64> = (0..n).map(|_| self.rng.gaussian()).collect();
212        let factor1 = (temperature.ln() / (QV - 1.0)).exp();
213        let factor4 = self.factor4_p * factor1;
214        let sigmax = (-(QV - 1.0) * (self.factor6 / factor4).ln() / (3.0 - QV)).exp();
215        (0..n)
216            .map(|i| {
217                let xi = x[i] * sigmax;
218                let den = ((y[i].abs() * (QV - 1.0)).ln() / (3.0 - QV)).exp();
219                xi / den
220            })
221            .collect()
222    }
223
224    fn visiting(&mut self, x: &[f64], step: usize, temperature: f64) -> Vec<f64> {
225        if step < self.dim {
226            let upper_sample = self.rng.uniform01();
227            let lower_sample = self.rng.uniform01();
228            let mut visits = self.visit_fn(temperature, self.dim);
229            for v in visits.iter_mut() {
230                if *v > TAIL_LIMIT {
231                    *v = TAIL_LIMIT * upper_sample;
232                } else if *v < -TAIL_LIMIT {
233                    *v = -TAIL_LIMIT * lower_sample;
234                }
235            }
236            let mut x_visit: Vec<f64> = (0..self.dim).map(|i| visits[i] + x[i]).collect();
237            for xv in x_visit.iter_mut() {
238                let b = (*xv % 1.0) + 1.0;
239                *xv = b % 1.0;
240                if xv.abs() < MIN_VISIT_BOUND {
241                    *xv += 1e-10;
242                }
243            }
244            x_visit
245        } else {
246            let mut x_visit = x.to_vec();
247            let mut visit = self.visit_fn(temperature, 1)[0];
248            if visit > TAIL_LIMIT {
249                visit = TAIL_LIMIT * self.rng.uniform01();
250            } else if visit < -TAIL_LIMIT {
251                visit = -TAIL_LIMIT * self.rng.uniform01();
252            }
253            let index = step - self.dim;
254            x_visit[index] = visit + x[index];
255            let b = (x_visit[index] % 1.0) + 1.0;
256            x_visit[index] = b % 1.0;
257            if x_visit[index].abs() < MIN_VISIT_BOUND {
258                x_visit[index] += MIN_VISIT_BOUND;
259            }
260            x_visit
261        }
262    }
263
264    // ---- energy state ----
265
266    fn reset_energy(&mut self, x0: &[f64]) {
267        self.current_location = if x0.is_empty() {
268            self.normal_vec()
269        } else {
270            x0.to_vec()
271        };
272        let mut reinit_counter = 0;
273        loop {
274            self.current_energy = self.value(&self.current_location.clone());
275            if self.current_energy >= BIG_VALUE || self.current_energy.is_nan() {
276                if reinit_counter >= MAX_REINIT_COUNT {
277                    self.reinit_failed = true;
278                    return;
279                }
280                self.current_location = self.uniform_vec();
281                reinit_counter += 1;
282            } else {
283                if self.ebest == f64::MAX && self.xbest.is_empty() {
284                    self.ebest = self.current_energy;
285                    self.xbest = self.current_location.clone();
286                }
287                return;
288            }
289        }
290    }
291
292    // ---- strategy chain ----
293
294    fn accept_reject(&mut self, j: usize, e: f64, x_visit: &[f64]) {
295        let r = self.rng.uniform01();
296        let pqv_temp = (QA - 1.0) * (e - self.current_energy) / (self.temperature_step + 1.0);
297        let pqv = if pqv_temp < 0.0 {
298            0.0
299        } else {
300            (pqv_temp.ln() / (1.0 - QA)).exp()
301        };
302        if r <= pqv {
303            self.current_energy = e;
304            self.current_location = x_visit.to_vec();
305            self.xmin = self.current_location.clone();
306        }
307        if self.not_improved_idx >= self.not_improved_max_idx
308            && (j == 0 || self.current_energy < self.emin)
309        {
310            self.emin = self.current_energy;
311            self.xmin = self.current_location.clone();
312        }
313    }
314
315    fn run_chain(&mut self, step: usize, temperature: f64) {
316        self.temperature_step = temperature / (step as f64 + 1.0);
317        self.not_improved_idx += 1;
318        let iters = self.current_location.len() * 2;
319        for j in 0..iters {
320            if j == 0 {
321                self.state_improved = false;
322            }
323            if step == 0 && j == 0 {
324                self.state_improved = true;
325            }
326            let x_visit = self.visiting(&self.current_location.clone(), j, temperature);
327            let e = self.value(&x_visit);
328            if e < self.current_energy {
329                self.current_energy = e;
330                self.current_location = x_visit.clone();
331                if e < self.ebest {
332                    self.ebest = e;
333                    self.xbest = x_visit.clone();
334                    self.state_improved = true;
335                    self.not_improved_idx = 0;
336                }
337            } else {
338                self.accept_reject(j, e, &x_visit);
339            }
340            if self.max_eval_reached() {
341                return;
342            }
343        }
344    }
345
346    fn chain_local_search(&mut self) {
347        if self.state_improved {
348            let (e, x) = self.local_search(&self.xbest.clone());
349            if e < self.ebest {
350                self.not_improved_idx = 0;
351                self.ebest = e;
352                self.xbest = x.clone();
353                self.current_energy = e;
354                self.current_location = x;
355                if self.max_eval_reached() {
356                    return;
357                }
358            }
359        }
360        let mut do_ls = false;
361        if self.k < 90.0 * self.dim as f64 {
362            let pls = (self.k * (self.ebest - self.current_energy) / self.temperature_step).exp();
363            if pls >= self.rng.uniform01() {
364                do_ls = true;
365            }
366        }
367        if self.not_improved_idx >= self.not_improved_max_idx {
368            do_ls = true;
369        }
370        if do_ls {
371            let (e, x) = self.local_search(&self.xmin.clone());
372            self.xmin = x.clone();
373            self.emin = e;
374            self.not_improved_idx = 0;
375            self.not_improved_max_idx = self.current_location.len() as i32;
376            if e < self.ebest {
377                self.ebest = e;
378                self.xbest = x.clone();
379                self.current_energy = e;
380                self.current_location = x;
381            }
382        }
383    }
384
385    // ---- bounded L-BFGS local search on [0,1]^dim ----
386
387    /// Finite-difference gradient matching the C++ `LBFGSFunc` (per-coordinate
388    /// forward/backward difference with boundary handling); returns `(grad, f)`.
389    fn fd_grad(&mut self, arg: &[f64]) -> (Vec<f64>, f64) {
390        let eps = 1e-6;
391        let mut grad = vec![0.0; self.dim];
392        for i in 0..self.dim {
393            let mut x1 = arg.to_vec();
394            let mut x2 = arg.to_vec();
395            let mut e1 = eps;
396            let mut e2 = eps;
397            x1[i] += eps;
398            if x1[i] > 1.0 {
399                x1[i] = 1.0;
400                e1 = 1.0 - arg[i];
401            }
402            x2[i] -= eps;
403            if x2[i] < 0.0 {
404                x2[i] = 0.0;
405                e2 = arg[i];
406            }
407            let f1 = self.value(&x1);
408            let f2 = self.value(&x2);
409            grad[i] = (f1 - f2) / (e1 + e2);
410        }
411        let f = self.value(arg);
412        (grad, f)
413    }
414
415    fn local_search(&mut self, x0: &[f64]) -> (f64, Vec<f64>) {
416        // reset the fitness-best tracker; the best point seen during the search
417        // (through `value`) is the returned result — matching the C++ contract.
418        self.fit_best_y = f64::MAX;
419        let mut max_iter = (6 * self.dim) as i32;
420        max_iter = max_iter.clamp(100, 1000);
421
422        let clamp01 = |v: &[f64]| -> Vec<f64> { v.iter().map(|x| x.clamp(0.0, 1.0)).collect() };
423        let mut x = clamp01(&self.closest_feasible(x0));
424        let m = 6usize;
425        let mut s_hist: Vec<Vec<f64>> = Vec::new();
426        let mut y_hist: Vec<Vec<f64>> = Vec::new();
427        let mut rho: Vec<f64> = Vec::new();
428
429        let (mut g, mut f) = self.fd_grad(&x);
430        for _ in 0..max_iter {
431            // projected-gradient stopping test
432            let pg_norm: f64 = (0..self.dim)
433                .map(|i| {
434                    let step = (x[i] - g[i]).clamp(0.0, 1.0) - x[i];
435                    step * step
436                })
437                .sum::<f64>()
438                .sqrt();
439            if pg_norm < 1e-10 {
440                break;
441            }
442            // two-loop recursion for the direction d = -H g
443            let mut q = g.clone();
444            let kh = s_hist.len();
445            let mut alpha = vec![0.0; kh];
446            for i in (0..kh).rev() {
447                let a = rho[i] * dot(&s_hist[i], &q);
448                alpha[i] = a;
449                for j in 0..self.dim {
450                    q[j] -= a * y_hist[i][j];
451                }
452            }
453            let gamma = if kh > 0 {
454                let last = kh - 1;
455                dot(&s_hist[last], &y_hist[last]) / dot(&y_hist[last], &y_hist[last])
456            } else {
457                1.0
458            };
459            for qi in q.iter_mut() {
460                *qi *= gamma;
461            }
462            for i in 0..kh {
463                let beta = rho[i] * dot(&y_hist[i], &q);
464                for j in 0..self.dim {
465                    q[j] += (alpha[i] - beta) * s_hist[i][j];
466                }
467            }
468            let d: Vec<f64> = q.iter().map(|v| -v).collect();
469
470            // projected backtracking line search (Armijo)
471            let gd = dot(&g, &d);
472            let mut step = 1.0;
473            let mut x_new = x.clone();
474            let mut f_new = f;
475            let mut ok = false;
476            for _ in 0..20 {
477                let cand: Vec<f64> = (0..self.dim)
478                    .map(|i| (x[i] + step * d[i]).clamp(0.0, 1.0))
479                    .collect();
480                let fc = self.value(&cand);
481                if fc.is_finite() && fc <= f + 1e-4 * step * gd {
482                    x_new = cand;
483                    f_new = fc;
484                    ok = true;
485                    break;
486                }
487                step *= 0.5;
488                if self.max_eval_reached() {
489                    break;
490                }
491            }
492            if !ok || self.max_eval_reached() {
493                break;
494            }
495
496            let (g_new, _) = self.fd_grad(&x_new);
497            let s: Vec<f64> = (0..self.dim).map(|i| x_new[i] - x[i]).collect();
498            let yv: Vec<f64> = (0..self.dim).map(|i| g_new[i] - g[i]).collect();
499            let sy = dot(&s, &yv);
500            if sy > 1e-12 {
501                if s_hist.len() == m {
502                    s_hist.remove(0);
503                    y_hist.remove(0);
504                    rho.remove(0);
505                }
506                s_hist.push(s);
507                y_hist.push(yv);
508                rho.push(1.0 / sy);
509            }
510            x = x_new;
511            g = g_new;
512            f = f_new;
513            if self.max_eval_reached() {
514                break;
515            }
516        }
517        (self.fit_best_y, self.fit_best_x.clone())
518    }
519
520    // ---- main annealing loop ----
521
522    fn search(&mut self) {
523        let mut iter = 0i32;
524        let t1 = ((QV - 1.0) * 2.0_f64.ln()).exp() - 1.0;
525        loop {
526            for i in 0..MAXSTEPS {
527                let s = i as f64 + 2.0;
528                let t2 = ((QV - 1.0) * s.ln()).exp() - 1.0;
529                let temperature = TEMPERATURE_START * t1 / t2;
530                iter += 1;
531                if iter >= MAXSTEPS {
532                    return;
533                }
534                if temperature < TEMPERATURE_RESTART {
535                    self.reset_energy(&[]);
536                    if self.reinit_failed {
537                        return;
538                    }
539                    break;
540                }
541                self.run_chain(i as usize, temperature);
542                if self.max_eval_reached() {
543                    return;
544                }
545                if self.use_local_search {
546                    self.chain_local_search();
547                    if self.max_eval_reached() {
548                        return;
549                    }
550                }
551            }
552        }
553    }
554
555    fn optimize(&mut self, guess: &[f64]) -> DaResult {
556        let enc = self.encode(guess);
557        self.reset_energy(&enc);
558        self.emin = self.current_energy;
559        self.xmin = self.current_location.clone();
560        self.not_improved_max_idx = 1000;
561        if !self.reinit_failed {
562            self.search();
563        }
564        DaResult {
565            x: self.decode(&self.xbest),
566            y: self.ebest,
567            evaluations: self.eval_counter,
568            iterations: 0,
569            stop: if self.reinit_failed { -1 } else { 0 },
570        }
571    }
572}
573
574fn dot(a: &[f64], b: &[f64]) -> f64 {
575    a.iter().zip(b).map(|(x, y)| x * y).sum()
576}
577
578/// Run Dual Annealing. `lower`/`upper` empty ⇒ unbounded.
579pub fn optimize_da(
580    obj: &impl Objective,
581    guess: &[f64],
582    lower: Vec<f64>,
583    upper: Vec<f64>,
584    p: &DaParams,
585) -> DaResult {
586    let dim = guess.len();
587    let max_evals = if p.max_evaluations == 0 {
588        10_000_000
589    } else {
590        p.max_evaluations
591    };
592    let mut da = Da::new(
593        obj,
594        dim,
595        lower,
596        upper,
597        max_evals,
598        p.use_local_search,
599        p.seed.wrapping_add(p.runid as u64),
600    );
601    da.optimize(guess)
602}
603
604#[cfg(test)]
605mod tests {
606    use super::*;
607
608    fn sphere(x: &[f64]) -> f64 {
609        x.iter().map(|v| v * v).sum()
610    }
611    fn rosen(x: &[f64]) -> f64 {
612        (0..x.len() - 1)
613            .map(|i| 100.0 * (x[i + 1] - x[i] * x[i]).powi(2) + (1.0 - x[i]).powi(2))
614            .sum()
615    }
616
617    fn run(obj: impl Objective, dim: usize, seed: u64, ls: bool) -> DaResult {
618        let params = DaParams {
619            max_evaluations: 40_000,
620            use_local_search: ls,
621            seed,
622            ..Default::default()
623        };
624        optimize_da(
625            &obj,
626            &vec![0.0; dim],
627            vec![-5.0; dim],
628            vec![5.0; dim],
629            &params,
630        )
631    }
632
633    #[test]
634    fn minimizes_sphere_with_local_search() {
635        let r = run(sphere as fn(&[f64]) -> f64, 4, 1, true);
636        assert!(r.y < 1e-6, "sphere not solved: {}", r.y);
637        assert!((sphere(&r.x) - r.y).abs() < 1e-6);
638    }
639
640    #[test]
641    fn minimizes_sphere_without_local_search() {
642        let r = run(sphere as fn(&[f64]) -> f64, 4, 2, false);
643        assert!(r.y < 1e-2, "sphere (no ls) too large: {}", r.y);
644    }
645
646    #[test]
647    fn minimizes_rosenbrock() {
648        let r = run(rosen as fn(&[f64]) -> f64, 3, 3, true);
649        assert!(r.y < 1e-2, "rosenbrock not solved: {}", r.y);
650    }
651
652    #[test]
653    fn unbounded_nonfinite_objective_reports_reinitialization_failure() {
654        let result = optimize_da(
655            &(|_: &[f64]| f64::NAN),
656            &[0.0],
657            Vec::new(),
658            Vec::new(),
659            &DaParams {
660                max_evaluations: 0,
661                use_local_search: false,
662                seed: 4,
663                runid: -1,
664            },
665        );
666        assert_eq!(result.stop, -1);
667        assert!(result.evaluations > 1_000);
668    }
669
670    #[test]
671    fn finite_difference_handles_both_box_boundaries() {
672        let objective = sphere as fn(&[f64]) -> f64;
673        let mut optimizer = Da::new(&objective, 2, vec![0.0; 2], vec![1.0; 2], 100, false, 5);
674        let (gradient, value) = optimizer.fd_grad(&[0.0, 1.0]);
675        assert!(gradient.iter().all(|component| component.is_finite()));
676        assert_eq!(value, 1.0);
677    }
678}
fcmaes_core/da.rs

fcmaes_core/
da.rs
