math_audio_optimisation/
impl_helpers.rs

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

math_audio_optimisation/
impl_helpers.rs
