math_audio_optimisation/
impl_helpers.rs

1use crate::{DEReport, DifferentialEvolution};
2use ndarray::{Array1, Array2, Zip};
3
4// ------------------------------ Internal helpers ------------------------------
5
6impl<'a, F> DifferentialEvolution<'a, F>
7where
8    F: Fn(&Array1<f64>) -> f64 + Sync,
9{
10    pub(crate) fn energy(&self, x: &Array1<f64>) -> f64 {
11        let base = (self.func)(x);
12        base + self.penalty(x)
13    }
14
15    pub(crate) fn penalty(&self, x: &Array1<f64>) -> f64 {
16        let mut p = 0.0;
17        // Nonlinear ineq: fc(x) <= 0 feasible
18        for (f, w) in &self.config.penalty_ineq {
19            let v = f(x);
20            let viol = v.max(0.0);
21            p += w * viol * viol;
22        }
23        // Nonlinear eq: h(x) = 0
24        for (h, w) in &self.config.penalty_eq {
25            let v = h(x);
26            p += w * v * v;
27        }
28        // Linear penalties: lb <= A x <= ub
29        if let Some(lp) = &self.config.linear_penalty {
30            let ax = lp.a.dot(&x.view());
31            Zip::from(&ax)
32                .and(&lp.lb)
33                .and(&lp.ub)
34                .for_each(|&v, &lo, &hi| {
35                    if v < lo {
36                        let d = lo - v;
37                        p += lp.weight * d * d;
38                    } else if v > hi {
39                        let d = v - hi;
40                        p += lp.weight * d * d;
41                    }
42                });
43        }
44        p
45    }
46
47    #[allow(clippy::too_many_arguments)]
48    pub(crate) fn finish_report(
49        &self,
50        pop: Array2<f64>,
51        energies: Array1<f64>,
52        x: Array1<f64>,
53        fun: f64,
54        success: bool,
55        message: String,
56        nit: usize,
57        nfev: usize,
58    ) -> DEReport {
59        DEReport {
60            x,
61            fun,
62            success,
63            message,
64            nit,
65            nfev,
66            population: pop,
67            population_energies: energies,
68        }
69    }
70
71    pub(crate) fn polish(&self, x0: &Array1<f64>) -> (Array1<f64>, f64, usize) {
72        let polish_cfg = match &self.config.polish {
73            Some(cfg) if cfg.enabled => cfg,
74            _ => {
75                let f = self.energy(x0);
76                return (x0.clone(), f, 1);
77            }
78        };
79
80        let n = x0.len();
81        let mut x = x0.clone();
82        let mut best_f = self.energy(&x);
83        let mut nfev = 1;
84
85        let initial_step = 0.1;
86        let min_step = 1e-8;
87        let mut step = initial_step;
88
89        let max_eval = polish_cfg.maxeval.min(200 * n);
90
91        while nfev < max_eval && step > min_step {
92            let mut improved = false;
93
94            for i in 0..n {
95                if nfev >= max_eval {
96                    break;
97                }
98
99                let original = x[i];
100                let bounds_span = self.upper[i] - self.lower[i];
101                let dim_step = step * bounds_span.max(1.0);
102
103                for delta in [dim_step, -dim_step] {
104                    if nfev >= max_eval {
105                        break;
106                    }
107                    x[i] = (original + delta).clamp(self.lower[i], self.upper[i]);
108                    let f = self.energy(&x);
109                    nfev += 1;
110
111                    if f < best_f {
112                        best_f = f;
113                        improved = true;
114                        break;
115                    }
116                    x[i] = original;
117                }
118            }
119
120            if !improved {
121                step *= 0.5;
122            }
123        }
124
125        (x, best_f, nfev)
126    }
127}
math_audio_optimisation/impl_helpers.rs

math_audio_optimisation/
impl_helpers.rs
