math-differential-evolution 0.3.3

Non linear optimisation library with own DE solvers and interface to NLOpt and MetaHeuristics
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127

use crate::{DEReport, DifferentialEvolution};
use ndarray::{Array1, Array2, Zip};

// ------------------------------ Internal helpers ------------------------------

impl<'a, F> DifferentialEvolution<'a, F>
where
    F: Fn(&Array1<f64>) -> f64 + Sync,
{
    pub(crate) fn energy(&self, x: &Array1<f64>) -> f64 {
        let base = (self.func)(x);
        base + self.penalty(x)
    }

    pub(crate) fn penalty(&self, x: &Array1<f64>) -> f64 {
        let mut p = 0.0;
        // Nonlinear ineq: fc(x) <= 0 feasible
        for (f, w) in &self.config.penalty_ineq {
            let v = f(x);
            let viol = v.max(0.0);
            p += w * viol * viol;
        }
        // Nonlinear eq: h(x) = 0
        for (h, w) in &self.config.penalty_eq {
            let v = h(x);
            p += w * v * v;
        }
        // Linear penalties: lb <= A x <= ub
        if let Some(lp) = &self.config.linear_penalty {
            let ax = lp.a.dot(&x.view());
            Zip::from(&ax)
                .and(&lp.lb)
                .and(&lp.ub)
                .for_each(|&v, &lo, &hi| {
                    if v < lo {
                        let d = lo - v;
                        p += lp.weight * d * d;
                    } else if v > hi {
                        let d = v - hi;
                        p += lp.weight * d * d;
                    }
                });
        }
        p
    }

    #[allow(clippy::too_many_arguments)]
    pub(crate) fn finish_report(
        &self,
        pop: Array2<f64>,
        energies: Array1<f64>,
        x: Array1<f64>,
        fun: f64,
        success: bool,
        message: String,
        nit: usize,
        nfev: usize,
    ) -> DEReport {
        DEReport {
            x,
            fun,
            success,
            message,
            nit,
            nfev,
            population: pop,
            population_energies: energies,
        }
    }

    pub(crate) fn polish(&self, x0: &Array1<f64>) -> (Array1<f64>, f64, usize) {
        let polish_cfg = match &self.config.polish {
            Some(cfg) if cfg.enabled => cfg,
            _ => {
                let f = self.energy(x0);
                return (x0.clone(), f, 1);
            }
        };

        let n = x0.len();
        let mut x = x0.clone();
        let mut best_f = self.energy(&x);
        let mut nfev = 1;

        let initial_step = 0.1;
        let min_step = 1e-8;
        let mut step = initial_step;

        let max_eval = polish_cfg.maxeval.min(200 * n);

        while nfev < max_eval && step > min_step {
            let mut improved = false;

            for i in 0..n {
                if nfev >= max_eval {
                    break;
                }

                let original = x[i];
                let bounds_span = self.upper[i] - self.lower[i];
                let dim_step = step * bounds_span.max(1.0);

                for delta in [dim_step, -dim_step] {
                    if nfev >= max_eval {
                        break;
                    }
                    x[i] = (original + delta).clamp(self.lower[i], self.upper[i]);
                    let f = self.energy(&x);
                    nfev += 1;

                    if f < best_f {
                        best_f = f;
                        improved = true;
                        break;
                    }
                    x[i] = original;
                }
            }

            if !improved {
                step *= 0.5;
            }
        }

        (x, best_f, nfev)
    }
}